diff options
| -rw-r--r-- | Makefile | 7 | ||||
| -rw-r--r-- | walk.cpp | 27 |
2 files changed, 30 insertions, 4 deletions
diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..7e1ed7d --- /dev/null +++ b/Makefile @@ -0,0 +1,7 @@ +all: walk + +CC := clang++ -std=c++17 -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native + +walk: walk.cpp + $(CC) walk.cpp -o walk + @@ -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); |