1 files changed, 23 insertions, 4 deletions

diff --git a/walk.cpp b/walk.cpp
index 2e22395..cd07253 100644
--- a/walk.cpp
+++ b/walk.cpp
@@ -5,11 +5,13 @@
 #include "randutils/randutils.hpp"
 
 #include "eigen/Eigen/Dense"
+#include "eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h"
 
 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 normalizeVector(const Vector& x) {
   return x * sqrt(x.size() / x.squaredNorm());
@@ -54,9 +56,26 @@ private:
 public:
   unsigned N;
 
-  QuadraticModel(unsigned N, Rng& r) : J(N), N(N) {
-    for (Real& Jᵢ : J) {
-      Jᵢ = wignerInverse(r.uniform(0.0, 1.0));
+  QuadraticModel(unsigned N, Rng& r, bool diag = false) : J(N), N(N) {
+    if (diag) {
+      Matrix Jtmp(N, N);
+
+      for (unsigned j = 0; j < N; j++) {
+        for (unsigned i = j; i < N; i++) {
+          Jtmp(i, j) = r.variate<Real, std::normal_distribution>(0, 1 / sqrt(N));
+          Jtmp(j, i) = Jtmp(i, j);
+        }
+      }
+
+      std::cerr << "Beginning diagonalization" << std::endl;
+      Eigen::SelfAdjointEigenSolver<Matrix> es;
+      es.compute(Jtmp);
+      J = es.eigenvalues();
+      std::cerr << "Finished diagonalization" << std::endl;
+    } else {
+      for (Real& Jᵢ : J) {
+        Jᵢ = wignerInverse(r.uniform(0.0, 1.0));
+      }
     }
   }
 
@@ -142,7 +161,7 @@ int main(int argc, char* argv[]) {
   }
 
   Rng r;
-  QuadraticModel model(N, r);
+  QuadraticModel model(N, r, true);
 
   Vector x₀ = normalizeVector(randomVector(N, r));
   x₀ = gradientDescent(model, x₀, E);