From bba01ec107ea1e0e6eccd79e4a8ebb30f46bfea7 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 25 Dec 2024 12:28:36 +0100
Subject: Wrote the program

---
 compile_flags.txt |   3 ++
 walk.cpp          | 128 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 131 insertions(+)
 create mode 100644 compile_flags.txt
 create mode 100644 walk.cpp

diff --git a/compile_flags.txt b/compile_flags.txt
new file mode 100644
index 0000000..4cb2489
--- /dev/null
+++ b/compile_flags.txt
@@ -0,0 +1,3 @@
+-std=c++17
+-Wno-mathematical-notation-identifier-extension
+-Wall
diff --git a/walk.cpp b/walk.cpp
new file mode 100644
index 0000000..3d7020d
--- /dev/null
+++ b/walk.cpp
@@ -0,0 +1,128 @@
+#include <getopt.h>
+#include <iomanip>
+
+#include "pcg-cpp/include/pcg_random.hpp"
+#include "randutils/randutils.hpp"
+
+#include "eigen/Eigen/Dense"
+
+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) {
+  return x * sqrt(x.size() / x.squaredNorm());
+}
+
+class QuadraticModel {
+private:
+  Matrix 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);
+        }
+      }
+    }
+  }
+
+  Real H(const Vector& x) const {
+    return 0.5 * (J * x).dot(x);
+  }
+
+  Vector ∇H(const Vector& x) const {
+    Vector ∂H = J * x;
+    return ∂H - (∂H.dot(x) / x.squaredNorm()) * x;
+  }
+};
+
+Vector gradientDescent(const QuadraticModel& M, const Vector& x₀, Real E, Real ε = 1e-14) {
+  Vector xₜ = x₀;
+  Real Hₜ = M.H(x₀);
+  Real α = 1;
+  Real m;
+  Vector ∇H;
+
+  while (
+    ∇H = (Hₜ / M.N - E) * M.∇H(xₜ) / M.N, m = ∇H.squaredNorm(),
+    m > ε
+  ) {
+    Vector xₜ₊₁;
+    Real Hₜ₊₁;
+
+    while (
+      xₜ₊₁ = normalize(xₜ - α * ∇H), Hₜ₊₁ = M.H(xₜ₊₁),
+      pow(Hₜ₊₁ / M.N - E, 2) > pow(Hₜ / M.N - E, 2) - α * m
+    ) {
+      α /= 2;
+    }
+
+    xₜ = xₜ₊₁;
+    Hₜ = Hₜ₊₁;
+    α *= 1.25;
+  }
+
+  return xₜ;
+}
+
+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 ∇H₀ = M.∇H(x₀);
+  η -= η.dot(∇H₀) * ∇H₀ / ∇H₀.squaredNorm();
+  return gradientDescent(M, normalize(x₀ + η), E);
+}
+
+int main(int argc, char* argv[]) {
+  unsigned N = 10;
+  unsigned steps = 10;
+  Real E = 0;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "N:E:n:")) != -1) {
+    switch (opt) {
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    case 'E':
+      E = atof(optarg);
+      break;
+    case 'n':
+      steps = (unsigned)atof(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  Rng r;
+  QuadraticModel ls(N, r);
+
+  Vector x₀ = Vector::Zero(N);
+  x₀(0) = sqrt(N);
+  x₀ = gradientDescent(ls, x₀, E);
+
+  std::cout << std::setprecision(15);
+
+  Vector x = x₀;
+
+  for (unsigned step = 0; step < steps; step++) {
+    x = randomStep(ls, x, E, r);
+    std::cout << ls.H(x) / N << " " << x.dot(x₀) / N << std::endl;
+  }
+
+  return 0;
+}
-- 
cgit v1.2.3-70-g09d2