From bfa9e59782d9c814568cca2e9ed56094d04c7bc7 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 6 Jan 2021 11:46:41 +0100
Subject: Better saddle-finding algorithm. Cristoffel symbols still need
 fixing.

---
 langevin.cpp      | 48 ++++++++++++++++++++++++++++++++----------------
 stereographic.hpp |  7 ++++---
 2 files changed, 36 insertions(+), 19 deletions(-)

diff --git a/langevin.cpp b/langevin.cpp
index 05b436e..acea609 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -26,16 +26,16 @@ Vector randomVector(unsigned N, Distribution d, Generator& r) {
   return z;
 }
 
-Vector findSaddle(const Tensor& J, const Vector& z0, double γ0, double δ, double ε, Rng& r) {
+std::tuple<double, Vector> gradientDescent(const Tensor& J, const Vector& z0, double γ0, double ε) {
   Vector z = z0;
 
   double W;
   Vector dW;
   std::tie(W, dW) = WdW(J, z);
 
-  while (W > δ) {
-    double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2);
+  double γ = γ0;
 
+  while (W > ε) {
     Vector zNew = z - γ * dW.conjugate();
     zNew *= sqrt(z.size()) / sqrt((Scalar)(zNew.transpose() * zNew));
 
@@ -47,32 +47,48 @@ Vector findSaddle(const Tensor& J, const Vector& z0, double γ0, double δ, doub
       z = zNew;
       W = WNew;
       dW = dWNew;
-      std::cout << W << std::endl;
-    } 
+      γ = γ0;
+    } else {
+      γ *= 0.5;
+    }
+
+    if (γ < 1e-15) {
+      std::cout << "Gradient descent stalled." << std::endl;
+      exit(1);
+    }
   }
 
-  Vector ζ = euclideanToStereographic(z);
+  return {W, z};
+}
+
+Vector findSaddle(const Tensor& J, const Vector& z0, double ε, double δW, double γ0) {
+  double W;
+  std::tie(W, std::ignore) = WdW(J, z0);
+
+  Vector ζ = euclideanToStereographic(z0);
   Vector dH;
   Matrix ddH;
   std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ);
 
   while (W > ε) {
-    double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2);
-
+    // ddH is complex symmetric, which is (almost always) invertible
     Vector dζ = ddH.partialPivLu().solve(dH);
     Vector ζNew = ζ - dζ;
 
     double WNew;
-    Vector dWNew;
-    std::tie(WNew, dWNew) = WdW(J, stereographicToEuclidean(ζNew));
+    std::tie(WNew, std::ignore) = WdW(J, stereographicToEuclidean(ζNew));
 
-    if (WNew < W) {
+    if (WNew < W) { // If Newton's step lowered the objective, accept it!
       ζ = ζNew;
       W = WNew;
-      dW = dWNew;
-      std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ);
-      std::cout << WNew << std::endl;
-    } 
+    } else { // Otherwise, do gradient descent until W is a factor δW smaller.
+      Vector z;
+      std::tie(W, z) = gradientDescent(J, stereographicToEuclidean(ζ), γ0, W / δW);
+      ζ = euclideanToStereographic(z);
+    }
+
+    std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ);
+    std::cout << W << std::endl;
   }
 
   return stereographicToEuclidean(ζ);
@@ -162,7 +178,7 @@ int main(int argc, char* argv[]) {
   };
 
   //Vector zm = langevin(J, z, T, γ, f, r);
-  Vector zm = findSaddle(J, z, γ, δ, ε, r);
+  Vector zm = findSaddle(J, z, ε, δ, γ);
 
   Scalar H;
   Vector dH;
diff --git a/stereographic.hpp b/stereographic.hpp
index 4c3c831..515890e 100644
--- a/stereographic.hpp
+++ b/stereographic.hpp
@@ -99,17 +99,18 @@ std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, con
   std::tie(hamiltonian, gradZ, hessZ) = hamGradHess(J, z);
 
   Matrix g = stereographicMetric(jacobian);
-  Matrix gInvJac = g.ldlt().solve(jacobian);
+  // g is Hermitian, which is positive definite.
+  Matrix gInvJac = g.llt().solve(jacobian);
 
   grad = gInvJac * gradZ;
 
   Eigen::Tensor<Scalar, 3> Γ = stereographicChristoffel(ζ, gInvJac);
   Vector df = jacobian * gradZ;
   Eigen::Tensor<Scalar, 1> dfT = Eigen::TensorMap<Eigen::Tensor<Scalar, 1>>(df.data(), {df.size()});
-  std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(1, 0)};
+  std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(0, 0)};
   Eigen::Tensor<Scalar, 2> H2 = Γ.contract(dfT, ip);
 
-  hess = jacobian * hessZ * jacobian.transpose() + Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size());
+  hess = jacobian * hessZ * jacobian.transpose() - Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size()).transpose();
 
   return {hamiltonian, grad, hess};
 }
-- 
cgit v1.2.3-70-g09d2