A lot of work, and fixed a huge bug regarding the meaning of .dot in the Eigen library for complex vectors.

3 files changed, 114 insertions, 41 deletions

diff --git a/langevin.cpp b/langevin.cpp
index 9da4ef6..3bc12ed 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -9,6 +9,7 @@
 
 #include "complex_normal.hpp"
 #include "p-spin.hpp"
+#include "stereographic.hpp"
 
 using Rng = randutils::random_generator<pcg32>;
 
@@ -23,6 +24,57 @@ Vector randomVector(unsigned N, Distribution d, Generator& r) {
   return z;
 }
 
+Vector findSaddle(const Tensor& J, const Vector& z0, double γ0, double δ, double ε, Rng& r) {
+  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);
+
+    Vector zNew = z - γ * dW.conjugate();
+    zNew *= sqrt(z.size()) / sqrt((Scalar)(zNew.transpose() * zNew));
+
+    double WNew;
+    Vector dWNew;
+    std::tie(WNew, dWNew) = WdW(J, zNew);
+
+    if (WNew < W) {
+      z = zNew;
+      W = WNew;
+      dW = dWNew;
+      std::cout << W << std::endl;
+    } 
+  }
+
+  Vector ζ = euclideanToStereographic(z);
+  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);
+
+    Vector ζNew = ζ - ddH.ldlt().solve(dH);
+
+    double WNew;
+    Vector dWNew;
+    std::tie(WNew, dWNew) = WdW(J, stereographicToEuclidean(ζNew));
+
+    if (WNew < W) {
+      ζ = ζNew;
+      W = WNew;
+      dW = dWNew;
+      std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ);
+      std::cout << WNew << std::endl;
+    } 
+  }
+
+  return stereographicToEuclidean(ζ);
+}
+
 Vector langevin(const Tensor& J, const Vector& z0, double T, double γ0,
                 std::function<bool(double, unsigned)> quit, Rng& r) {
   Vector z = z0;
@@ -32,24 +84,16 @@ Vector langevin(const Tensor& J, const Vector& z0, double T, double γ0,
   std::tie(W, dW) = WdW(J, z);
 
   unsigned steps = 0;
-  complex_normal_distribution<> d(0, T, 0);
+  complex_normal_distribution<> d(0, sqrt(T), 0);
 
   while (!quit(W, steps)) {
     double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2);
     Vector η = randomVector(z.size(), d, r.engine());
 
-    Vector zNext = z - γ * dW + η;
-    zNext *= sqrt(zNext.size()) / sqrt(zNext.dot(zNext));
+    z -= γ * (dW.conjugate() / pow((double)z.size(), 2) + η);
+    z *= sqrt(z.size()) / sqrt((Scalar)(z.transpose() * z));
 
-    double WNext;
-    Vector dWNext;
-    std::tie(WNext, dWNext) = WdW(J, zNext);
-
-    if (exp((W - WNext) / T) > r.uniform(0.0, 1.0)) {
-      z = zNext;
-      W = WNext;
-      dW = dWNext;
-    }
+    std::tie(W, dW) = WdW(J, z);
   }
 
   return z;
@@ -63,13 +107,14 @@ int main(int argc, char* argv[]) {
   double Iκ = 0;   // imaginary part of distribution parameter
 
   // simulation parameters
+  double ε = 1e-4;
   double δ = 1e-2;   // threshold for determining saddle
   double γ = 1e-2;   // step size
   unsigned t = 1000; // number of Langevin steps
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "N:T:e:r:i:g:t:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:T:e:r:i:g:t:E:")) != -1) {
     switch (opt) {
     case 'N':
       N = (unsigned)atof(optarg);
@@ -83,6 +128,9 @@ int main(int argc, char* argv[]) {
     case 'e':
       δ = atof(optarg);
       break;
+    case 'E':
+      ε = atof(optarg);
+      break;
     case 'g':
       γ = atof(optarg);
     case 'r':
@@ -103,14 +151,15 @@ int main(int argc, char* argv[]) {
 
   Tensor J = generateCouplings<Scalar, PSPIN_P>(N, complex_normal_distribution<>(0, σ, κ), r.engine());
   Vector z = randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine());
-  z *= sqrt(N) / sqrt(z.dot(z)); // Normalize.
+  z *= sqrt(N) / sqrt((Scalar)(z.transpose() * z)); // Normalize.
 
-  std::function<bool(double, unsigned)> findSaddle = [δ](double W, unsigned) {
+  std::function<bool(double, unsigned)> f = [δ](double W, unsigned) {
     std::cout << W << std::endl;
     return W < δ;
   };
 
-  Vector zm = langevin(J, z, T, γ, findSaddle, r);
+  //Vector zm = langevin(J, z, T, γ, f, r);
+  Vector zm = findSaddle(J, z, γ, δ, ε, r);
 
   Scalar H;
   Vector dH;
@@ -119,5 +168,6 @@ int main(int argc, char* argv[]) {
 
   Vector constraint = dH - ((double)p * H / (double)N) * zm;
 
+  std::cout << constraint << std::endl;
   std::cout << H / (double)N << std::endl;
 }
diff --git a/p-spin.hpp b/p-spin.hpp
index 1ed4d8e..a522aee 100644
--- a/p-spin.hpp
+++ b/p-spin.hpp
@@ -24,7 +24,7 @@ std::tuple<Scalar, Vector, Matrix> hamGradHess(const Tensor& J, const Vector& z)
 
   Matrix hessian = ((p - 1) * p / f) * Jz;
   Vector gradient = (p / f) * Jzz;
-  Scalar hamiltonian = (1 / f) * Jzz.dot(z);
+  Scalar hamiltonian = (1 / f) * Jzz.transpose() * z;
 
   return {hamiltonian, gradient, hessian};
 }
@@ -34,10 +34,10 @@ std::tuple<double, Vector> WdW(const Tensor& J, const Vector& z) {
   Matrix hessian;
   std::tie(std::ignore, gradient, hessian) = hamGradHess(J, z);
 
-  Vector projectedGradient = gradient - (gradient.dot(z) / (double)z.size()) * z;
+  Vector projectedGradient = (gradient - ((Scalar)(gradient.transpose() * z) / (double)z.size()) * z).conjugate();
 
   double W = projectedGradient.cwiseAbs2().sum();
-  Vector dW = hessian.conjugate() * projectedGradient;
+  Vector dW = hessian * projectedGradient - ((z.transpose() * gradient) * projectedGradient + (z.transpose() * projectedGradient) * (gradient + hessian * z)) / (double)z.size();
 
   return {W, dW};
 }
diff --git a/stereographic.hpp b/stereographic.hpp
index 4109bc7..59432b7 100644
--- a/stereographic.hpp
+++ b/stereographic.hpp
@@ -1,12 +1,14 @@
 #include <eigen3/Eigen/Cholesky>
 
+#include "Eigen/src/Core/util/Meta.h"
 #include "p-spin.hpp"
+#include "unsupported/Eigen/CXX11/src/Tensor/TensorMeta.h"
 
 Vector stereographicToEuclidean(const Vector& ζ) {
   unsigned N = ζ.size() + 1;
   Vector z(N);
-  Scalar a = ζ.dot(ζ);
-  std::complex<double> b = 2.0 * sqrt(N) / (1.0 + a);
+  Scalar a = ζ.transpose() * ζ;
+  Scalar b = 2 * sqrt(N) / (1.0 + a);
 
   for (unsigned i = 0; i < N - 1; i++) {
     z(i) = b * ζ(i);
@@ -32,15 +34,14 @@ Matrix stereographicJacobian(const Vector& ζ) {
   unsigned N = ζ.size() + 1;
   Matrix J(N - 1, N);
 
-  Scalar b = ζ.dot(ζ);
+  Scalar b = ζ.transpose() * ζ;
 
   for (unsigned i = 0; i < N - 1; i++) {
     for (unsigned j = 0; j < N - 1; j++) {
-      J(i, j) = -4.0 * ζ(i) * ζ(j);
+      J(i, j) = - 4 * sqrt(N) * ζ(i) * ζ(j) / pow(1.0 + b, 2);
       if (i == j) {
-        J(i, j) += 2.0 * (1.0 + b);
+        J(i, j) += 2 * sqrt(N) * (1.0 + b) / pow(1.0 + b, 2);
       }
-      J(i, j) *= sqrt(N) / pow(1.0 + b, 2);
     }
 
     J(i, N - 1) = 4.0 * sqrt(N) * ζ(i) / pow(1.0 + b, 2);
@@ -49,25 +50,40 @@ Matrix stereographicJacobian(const Vector& ζ) {
   return J;
 }
 
-Matrix stereographicMetric(const Vector& ζ) {
-  unsigned N = ζ.size();
-  Matrix g(N, N);
+Matrix stereographicMetric(const Matrix& J) {
+  return J * J.adjoint();
+}
+
+Eigen::Tensor<Scalar, 3> stereographicChristoffel(const Vector& ζ, const Matrix& gInvJac) {
+  unsigned N = ζ.size() + 1;
+  Eigen::Tensor<Scalar, 3> dJ(N - 1, N - 1, N);
 
-  double a = ζ.cwiseAbs2().sum();
-  Scalar b = ζ.dot(ζ);
+  Scalar b = ζ.transpose() * ζ;
 
-  for (unsigned i = 0; i < N; i++) {
-    for (unsigned j = 0; j < N; j++) {
-      g(i, j) = 16.0 * (std::conj(ζ(i)) * ζ(j) * (1.0 + a) -
-                        std::real(std::conj(ζ(i) * ζ(j)) * (1.0 + b)));
+  for (unsigned i = 0; i < N - 1; i++) {
+    for (unsigned j = 0; j < N - 1; j++) {
+      for (unsigned k = 0; k < N - 1; k++) {
+        dJ(i, j, k) = 16 * sqrt(N) * ζ(i) * ζ(j) * ζ(k) / pow(1.0 + b, 3);
+        if (i == j) {
+          dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(k) / pow(1.0 + b, 3);
+        }
+        if (i == k) {
+          dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(j) / pow(1.0 + b, 3);
+        }
+        if (j == k) {
+          dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(i) / pow(1.0 + b, 3);
+        }
+      }
+      dJ(i, j, N - 1) = - 16 * sqrt(N) * ζ(i) * ζ(j) / pow(1.0 + b, 3);
       if (i == j) {
-        g(i, j) += 4.0 * std::abs(1.0 + b);
+        dJ(i, j, N - 1) += 4 * sqrt(N) * (1.0 + b) * ζ(i) / pow(1.0 + b, 3);
       }
-      g(i, j) *= (N + 1) / std::norm(pow(b, 2));
     }
   }
 
-  return g;
+  std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(2, 1)};
+
+  return dJ.contract(Eigen::TensorMap<Eigen::Tensor<const Scalar, 2>>(gInvJac.data(), N - 1, N), ip);
 }
 
 std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, const Vector& ζ) {
@@ -82,11 +98,18 @@ std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, con
   Matrix hessZ;
   std::tie(hamiltonian, gradZ, hessZ) = hamGradHess(J, z);
 
-  Matrix g = stereographicMetric(ζ);
-  Matrix gj = g.llt().solve(jacobian);
+  Matrix g = stereographicMetric(jacobian);
+  Matrix gInvJac = g.ldlt().solve(jacobian);
+
+  grad = gInvJac * gradZ;
+
+  Eigen::Tensor<Scalar, 3> Γ = stereographicChristoffel(ζ, gInvJac.conjugate());
+  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>(0, 0)};
+  Eigen::Tensor<Scalar, 2> H2 = Γ.contract(dfT, ip);
 
-  grad = gj * gradZ;
-  hess = gj * hessZ * gj.transpose();
+  hess = jacobian * hessZ * jacobian.transpose() + (Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size())).transpose();
 
   return {hamiltonian, grad, hess};
 }