A lot of uncompleted work.

1 files changed, 185 insertions, 57 deletions

diff --git a/langevin.cpp b/langevin.cpp
index 241abaf..b6f16d7 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -1,11 +1,15 @@
+#include <algorithm>
 #include <getopt.h>
 #include <cstdlib>
 #include <random>
 #include <complex>
 #include <array>
 #include <functional>
+#include <list>
 
 #include <eigen3/unsupported/Eigen/CXX11/Tensor>
+#include <eigen3/Eigen/Core>
+#include <eigen3/Eigen/Cholesky>
 
 #include "pcg-cpp/include/pcg_random.hpp"
 #include "randutils/randutils.hpp"
@@ -13,10 +17,14 @@
 #define P_SPIN_P 3
 const unsigned p = P_SPIN_P;
 
-using Tensor = Eigen::Tensor<std::complex<double>, p>;
-using Scalar = Eigen::Tensor<std::complex<double>, 0>;
-using Vector = Eigen::Tensor<std::complex<double>, 1>;
-using Matrix = Eigen::Tensor<std::complex<double>, 2>;
+using Scalar = std::complex<double>;
+using TensorP = Eigen::Tensor<Scalar, p>;
+using Tensor0 = Eigen::Tensor<Scalar, 0>;
+using Tensor1 = Eigen::Tensor<Scalar, 1>;
+using Tensor2 = Eigen::Tensor<Scalar, 2>;
+using Tensor3 = Eigen::Tensor<Scalar, 3>;
+using Matrix = Eigen::MatrixXcd;
+using Vector = Eigen::VectorXcd;
 
 const std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(0,0)};
 
@@ -49,11 +57,21 @@ long unsigned factorial(unsigned p) {
   }
 }
 
-Tensor generateCouplings(unsigned N, std::complex<double> κ, Rng& r) {
+/*
+void populateCouplings(TensorP& J, unsigned level, std::list<unsigned> indices, complex_normal_distribution<> dist, Rng& r) {
+  if (level == 0) {
+    Scalar x = dist(r.engine());
+    do {
+    } while (std::next_permutation(indices.begin(), indices.end()))
+  }
+}
+*/
+
+Tensor3 generateCouplings(unsigned N, std::complex<double> κ, Rng& r) {
   double σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1)));
   complex_normal_distribution<> dist(0, σ, κ);
 
-  Tensor J(N, N, N);
+  TensorP J(N, N, N);
 
   for (unsigned i1 = 0; i1 < N; i1++) {
     for (unsigned i2 = i1; i2 < N; i2++) {
@@ -84,70 +102,186 @@ double norm(const Eigen::Tensor<std::complex<double>, r>& t) {
   return t2(0);
 }
 
-Vector initializeVector(unsigned N) {
+Vector initializeVector(unsigned N, Rng& r) {
   Vector z(N);
-  z.setConstant(1);
+  complex_normal_distribution<> dist(0, 1, 0);
+
+  for (unsigned i = 0; i < N; i++) {
+    z(i) = dist(r.engine());
+  }
+
+  z *= sqrt(N) / sqrt(z.dot(z));
+
   return z;
 }
 
-std::complex<double> hamiltonian(const Tensor& J, const Vector& z) {
-  Eigen::Tensor<std::complex<double>, 0> t = z.contract(z.contract(z.contract(J, ip), ip), ip);
-  return t(0) / (double)factorial(p);
+template <int r>
+Eigen::Tensor<Scalar, 3> contractDown(const Eigen::Tensor<Scalar, r>& J, const Eigen::Tensor<Scalar, 1>& z) {
+  return contractDown<r - 1>(J.contract(z, ip), z);
 }
 
-Matrix identity(unsigned N) {
-  Matrix I(N, N);
-  for (unsigned i = 0; i < N; i++) {
-    I(i, i) = 1;
+template <>
+Eigen::Tensor<Scalar, 3> contractDown(const Eigen::Tensor<Scalar, 3>& J, const Eigen::Tensor<Scalar, 1>& z) {
+  return J;
+}
+
+std::tuple<std::complex<double>, Vector, Matrix> hamGradHess(const Tensor3& J, const Vector& z) {
+  Tensor1 zT = Eigen::TensorMap<Eigen::Tensor<const std::complex<double>, 1>>(z.data(), {z.size()});
+
+  Tensor2 J1 = zT.contract(J, ip);
+  Tensor1 J2 = zT.contract(J1, ip);
+  Tensor0 J3 = zT.contract(J2, ip);
+
+  Matrix hess = ((p - 1) * (double)p / factorial(p)) * Eigen::Map<Matrix>(J1.data(), z.size(), z.size());
+  Vector grad = ((double)p / factorial(p)) * Eigen::Map<Vector>(J2.data(), z.size());
+  std::complex<double> hamiltonian = (1.0 / factorial(p)) * J3(0);
+
+  return {hamiltonian, grad, hess};
+}
+
+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);
+
+  for (unsigned i = 0; i < N - 1; i++) {
+    z(i) = b * ζ(i);
+  }
+
+  z(N - 1) = b * (a - 1.0) / 2.0;
+
+  return z;
+}
+
+Vector euclideanToStereographic(const Vector& z) {
+  unsigned N = z.size();
+  Vector ζ(N - 1);
+  
+  for (unsigned i = 0; i < N - 1; i++) {
+    ζ(i) = z(i) / (sqrt(N) - z(N - 1));
+  }
+
+  return ζ;
+}
+
+Matrix stereographicJacobian(const Vector& ζ) {
+  unsigned N = ζ.size() + 1;
+  Matrix J(N - 1, N);
+
+  Scalar b = ζ.dot(ζ);
+
+  for (unsigned i = 0; i < N - 1; i++) {
+    for (unsigned j = 0; j < N - 1; j++) {
+      J(i, j) = - 4.0 * ζ(i) * ζ(j);
+      if (i == j) {
+        J(i, j) += 2.0 * (1.0 + b);
+      }
+      J(i, j) *= sqrt(N) / pow(1.0 + b, 2);
+    }
+
+    J(i, N - 1) = 4.0 * sqrt(N) * ζ(i) / pow(1.0 + b, 2);
   }
-  return I;
+
+  return J;
 }
 
-std::tuple<Vector, Matrix> gradHess(const Tensor& J, const Vector& z, std::complex<double> ε) {
-  Matrix J1 = z.contract(J, ip);
-  Matrix hess = ((p - 1) * (double)p / factorial(p)) * J1 - ((double)p * ε) * identity(z.size());
+Matrix stereographicMetric(const Vector& ζ) {
+  unsigned N = ζ.size();
+  Matrix g(N, N);
 
-  Vector J2 = z.contract(J1, ip);
-  Vector grad = ((double)p / factorial(p)) * J2 - ((double)p * ε) * z;
+  double a = ζ.cwiseAbs2().sum();
+  Scalar b = ζ.dot(ζ);
 
-  return {grad, hess};
+  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)));
+      if (i == j) {
+        g(i, j) += 4.0 * std::abs(1.0 + b);
+      }
+      g(i, j) *= (N + 1) / std::norm(pow(b, 2));
+    }
+  }
+
+  return g;
 }
 
-std::tuple<double, Vector, std::complex<double>, Vector> WdW(const Tensor& J, const Vector& z, std::complex<double> ε) {
+std::tuple<std::complex<double>, Vector, Matrix> stereographicHamGradHess(const Tensor3& J, const Vector& ζ) {
   Vector grad;
   Matrix hess;
 
-  std::tie(grad, hess) = gradHess(J, z, ε);
+  Matrix jacobian = stereographicJacobian(ζ);
+  Vector z = stereographicToEuclidean(ζ);
 
-  double W = norm(grad);
-  Matrix conjHess = conj(hess);
-  Scalar zz = z.pow(2).sum();
-  Vector zc = conj(z);
+  std::complex<double> hamiltonian;
+  Vector gradZ;
+  Matrix hessZ;
+  std::tie(hamiltonian, gradZ, hessZ) = hamGradHess(J, z);
 
-  Vector dWdz = grad.contract(conjHess, ip) - (pow(p, 2) / 2.0 * ((double)z.size() - zz(0))) * zc;
-  Scalar dWdε = (-(double)p) * grad.contract(zc, ip);
-  Vector dεdz = (1 / (double)z.size()) * ((1 - 1/(double)p) * grad + ((double)p * ε) * z - (1 /(double)p) * z.contract(hess, ip));
+  Matrix g = stereographicMetric(ζ);
+  Matrix gj = g.llt().solve(jacobian);
 
-  return {W, dWdz, dWdε(0), dεdz};
+  grad = gj * gradZ;
+  hess = gj * hessZ * gj.transpose();
 
+  return {hamiltonian, grad, hess};
 }
 
-void gradientDescent(const Tensor& J, Vector& z, std::complex<double>& ε, double δ, double γ) {
+std::tuple<double, Vector, Matrix> WdW(const Tensor3& J, const Vector& z) {
+  Vector grad;
+  Matrix hess;
+
+  std::tie(std::ignore, grad, hess) = hamGradHess(J, z);
+
+  Tensor1 gradT = Eigen::TensorMap<Eigen::Tensor<const Scalar, 1>>(grad.data(), {z.size()});
+
+  Vector gtmp = grad - (grad.dot(z) / (double)z.size()) * z;
+  double W = gtmp.cwiseAbs2().sum();
+  Vector dW = (hess - (z.dot(grad) * Matrix::Identity(z.size(), z.size()) + grad * z.transpose() + z * (hess * z).transpose()) / (double)z.size()).conjugate() * gtmp;
+  Tensor2 ddWT = ((p - 2) * (p - 1) * (double)p / factorial(p)) * conj(J).contract(gradT, ip);
+  Matrix ddW = Eigen::Map<Matrix>(ddWT.data(), z.size(), z.size());
+
+  return {W, dW, ddW};
+}
+
+Vector findSaddle(const Tensor3& J, const Vector& z0, double δ, double γ) {
+  Vector z = z0;
   double W;
-  Vector dz;
-  std::complex<double> dε;
+  Vector dW;
+  Matrix ddW;
 
-  std::tie(W, dz, dε, std::ignore) = WdW(J, z, ε);
+  std::tie(W, dW, ddW) = WdW(J, z);
 
   while (W > δ * z.size()) {
-    std::cout << W << std::endl;
-    z = z - (γ / z.size()) * dz;
-    ε = ε - (γ / z.size()) * dε;
-
-    std::tie(W, dz, dε, std::ignore) = WdW(J, z, ε);
+ //   Vector dz = ddW.ldlt().solve(dW);
+    Vector dz = dW;
+    while (true) {
+      Vector newz = z - γ * dz;
+      newz *= sqrt(newz.size()) / sqrt(newz.dot(newz));
+
+      double newW;
+      Vector newdW;
+      Matrix newddW;
+      std::tie(newW, newdW, newddW) = WdW(J, newz);
+
+      if (newW < W) {
+        γ *= 1.01;
+        z = newz;
+        W = newW;
+        dW= newdW;
+        ddW = newddW;
+        break;
+      } else {
+        γ /= 0.9;
+      }
+    }
+  std::cout << W << std::endl;
   }
+
+  return z;
 }
 
+/*
 Vector generateKick(const Vector& z, Rng& r) {
   Vector k(z.size());
 
@@ -164,9 +298,8 @@ Vector generateKick(const Vector& z, Rng& r) {
 void langevin(const Tensor& J, Vector& z, std::complex<double>& ε, double T, unsigned t, double γ, Rng& r) {
   double W;
   Vector dz;
-  Vector dεdz;
 
-  std::tie(W, dz, std::ignore, dεdz) = WdW(J, z, ε);
+  std::tie(W, dz) = WdW(J, z);
 
   for (unsigned i = 0; i < t; i++) {
     std::cout << W << std::endl;
@@ -181,6 +314,7 @@ void langevin(const Tensor& J, Vector& z, std::complex<double>& ε, double T, un
     std::tie(W, dz, std::ignore, dεdz) = WdW(J, z, ε);
   }
 }
+*/
 
 int main(int argc, char* argv[]) {
   // model parameters
@@ -227,23 +361,17 @@ int main(int argc, char* argv[]) {
 
   Rng r;
 
-  Tensor J = generateCouplings(N, κ, r);
-  Vector z = initializeVector(N);
-  std::complex<double> ε = hamiltonian(J, z) / (double)N;
-
-  Scalar zz = z.pow(2).sum();
+  Tensor3 J = generateCouplings(N, κ, r);
+  Vector z = initializeVector(N, r);
 
-  std::cout << zz(0) << std::endl;
+  Vector zm = findSaddle(J, z, δ, γ);
 
-  gradientDescent(J, z, ε, δ, γ);
-  langevin(J, z, ε, T, t, γ, r);
+  std::complex<double> H;
+  Vector dH;
 
-  zz = z.pow(2).sum();
-  double a = norm(z);
+  std::tie(H, dH, std::ignore) = hamGradHess(J, zm);
 
-  std::cout << zz(0) / (double)N << std::endl;
-  std::cout << ε << std::endl;
-  std::cout << hamiltonian(J, z) / (double)N << std::endl;
-  std::cout << a / N << std::endl;
+  std::cout << zm << std::endl;
+  std::cout << dH - (H / (double)N) * zm << std::endl;
 }