1 files changed, 20 insertions, 20 deletions

diff --git a/walk.cpp b/walk.cpp
index ab4dbd3..f064bb5 100644
--- a/walk.cpp
+++ b/walk.cpp
@@ -11,12 +11,19 @@ using Rng = randutils::random_generator<pcg32>;
 
 using Real = double;
 using Vector = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
-using Matrix = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
 
-Vector normalize(const Vector& x) {
+Vector normalizeVector(const Vector& x) {
   return x * sqrt(x.size() / x.squaredNorm());
 }
 
+Vector randomVector(unsigned N, Rng& r, Real σ = 1) {
+  Vector v(N);
+  for (Real& vᵢ : v) {
+    vᵢ = r.variate<Real, std::normal_distribution>(0, σ);
+  }
+  return v;
+}
+
 Real wignerCDF(Real λ) {
   return 0.5 + (λ * sqrt(4 - pow(λ, 2)) / 4 + atan(λ / sqrt(4 - pow(λ, 2)))) / M_PI;
 }
@@ -63,7 +70,7 @@ public:
 Vector gradientDescent(const QuadraticModel& M, const Vector& x₀, Real E, Real ε = 1e-14) {
   Vector xₜ = x₀;
   Real Hₜ = M.H(x₀);
-  Real α = 1;
+  Real α = 1.0 / M.N;
   Real m;
   Vector ∇H;
 
@@ -75,7 +82,7 @@ Vector gradientDescent(const QuadraticModel& M, const Vector& x₀, Real E, Real
     Real Hₜ₊₁;
 
     while (
-      xₜ₊₁ = normalize(xₜ - α * ∇H), Hₜ₊₁ = M.H(xₜ₊₁),
+      xₜ₊₁ = normalizeVector(xₜ - α * ∇H), Hₜ₊₁ = M.H(xₜ₊₁),
       pow(Hₜ₊₁ / M.N - E, 2) > pow(Hₜ / M.N - E, 2) - α * m
     ) {
       α /= 2;
@@ -90,14 +97,11 @@ Vector gradientDescent(const QuadraticModel& M, const Vector& x₀, Real E, Real
 }
 
 Vector randomStep(const QuadraticModel& M, const Vector& x₀, Real E, Rng& r, Real Δt = 1e-4) {
-  Vector η(M.N);
-  for (Real& ηᵢ : η) {
-    ηᵢ = r.variate<Real, std::normal_distribution>(0, sqrt(2 * Δt));
-  }
-  η -= η.dot(x₀) * x₀ / x₀.squaredNorm();
+  Vector η = randomVector(M.N, r, sqrt(2 * Δt));
+  η -= (η.dot(x₀) / x₀.squaredNorm()) * x₀;
   Vector ∇H₀ = M.∇H(x₀);
-  η -= η.dot(∇H₀) * ∇H₀ / ∇H₀.squaredNorm();
-  return gradientDescent(M, normalize(x₀ + η), E);
+  η -= (η.dot(∇H₀) / ∇H₀.squaredNorm()) * ∇H₀;
+  return gradientDescent(M, normalizeVector(x₀ + η), E);
 }
 
 int main(int argc, char* argv[]) {
@@ -132,26 +136,22 @@ int main(int argc, char* argv[]) {
   }
 
   Rng r;
-  QuadraticModel ls(N, r);
+  QuadraticModel model(N, r);
 
-  Vector x₀ = Vector::Zero(N);
-  for (Real& xi : x₀) {
-    xi = r.variate<Real, std::normal_distribution>(0, 1);
-  }
-  x₀ = normalize(x₀);
-  x₀ = gradientDescent(ls, x₀, E);
+  Vector x₀ = normalizeVector(randomVector(N, r));
+  x₀ = gradientDescent(model, x₀, E);
 
   std::cout << std::setprecision(15);
 
   for (Real t = 0; t < T₀; t += Δt) {
-    x₀ = randomStep(ls, x₀, E, r, Δt);
+    x₀ = randomStep(model, x₀, E, r, Δt);
   }
 
   Vector x = x₀;
 
   for (Real t = 0; t < T; t += Δt) {
     std::cout << t << " " << x.dot(x₀) / N << std::endl;
-    x = randomStep(ls, x, E, r, Δt);
+    x = randomStep(model, x, E, r, Δt);
   }
 
   return 0;