1 files changed, 31 insertions, 13 deletions

diff --git a/walk.cpp b/walk.cpp
index 3d7020d..7c5d01c 100644
--- a/walk.cpp
+++ b/walk.cpp
@@ -1,5 +1,6 @@
 #include <getopt.h>
 #include <iomanip>
+#include <random>
 
 #include "pcg-cpp/include/pcg_random.hpp"
 #include "randutils/randutils.hpp"
@@ -16,31 +17,45 @@ Vector normalize(const Vector& x) {
   return x * sqrt(x.size() / x.squaredNorm());
 }
 
+Real wignerCDF(Real λ) {
+  return 0.5 + (λ * sqrt(4 - pow(λ, 2)) / 4 + atan(λ / sqrt(4 - pow(λ, 2)))) / M_PI;
+}
+
+Real wignerInverse(Real p, Real ε = 1e-14) {
+  Real a = -2;
+  Real b = 2;
+
+  while (b - a > ε) {
+    Real c = (a + b) / 2;
+    if ((wignerCDF(a) - p) * (wignerCDF(c) - p) > 0) {
+      a = c;
+    } else {
+      b = c;
+    }
+  }
+
+  return (a + b) / 2;
+}
+
 class QuadraticModel {
 private:
-  Matrix J;
+  Vector J;
 
 public:
   unsigned N;
 
-  QuadraticModel(unsigned N, Rng& r) : J(N, N), N(N) {
-
-    for (unsigned i = 0; i < N; i++) {
-      for (unsigned j = i; j < N; j++) {
-        J(i, j) = r.variate<Real, std::normal_distribution>(0, 1 / sqrt(N));
-        if (i != j) {
-          J(j, i) = J(i, j);
-        }
-      }
+  QuadraticModel(unsigned N, Rng& r) : J(N), N(N) {
+    for (Real& Jᵢ : J) {
+      Jᵢ = wignerInverse(r.uniform(0.0, 1.0));
     }
   }
 
   Real H(const Vector& x) const {
-    return 0.5 * (J * x).dot(x);
+    return 0.5 * (J.cwiseProduct(x)).dot(x);
   }
 
   Vector ∇H(const Vector& x) const {
-    Vector ∂H = J * x;
+    Vector ∂H = J.cwiseProduct(x);
     return ∂H - (∂H.dot(x) / x.squaredNorm()) * x;
   }
 };
@@ -112,7 +127,10 @@ int main(int argc, char* argv[]) {
   QuadraticModel ls(N, r);
 
   Vector x₀ = Vector::Zero(N);
-  x₀(0) = sqrt(N);
+  for (Real& xi : x₀) {
+    xi = r.variate<Real, std::normal_distribution>(0, 1);
+  }
+  x₀ = normalize(x₀);
   x₀ = gradientDescent(ls, x₀, E);
 
   std::cout << std::setprecision(15);