Got Newton's method working well.

1 files changed, 10 insertions, 9 deletions

diff --git a/stereographic.hpp b/stereographic.hpp
index 515890e..a3f2cc6 100644
--- a/stereographic.hpp
+++ b/stereographic.hpp
@@ -54,7 +54,8 @@ Matrix stereographicMetric(const Matrix& J) {
   return J * J.adjoint();
 }
 
-Eigen::Tensor<Scalar, 3> stereographicChristoffel(const Vector& ζ, const Matrix& gInvJac) {
+// Gives the Christoffel symbol, with Γ_(i1, i2)^(i3).
+Eigen::Tensor<Scalar, 3> stereographicChristoffel(const Vector& ζ, const Matrix& gInvJacConj) {
   unsigned N = ζ.size() + 1;
   Eigen::Tensor<Scalar, 3> dJ(N - 1, N - 1, N);
 
@@ -83,7 +84,7 @@ Eigen::Tensor<Scalar, 3> stereographicChristoffel(const Vector& ζ, const Matrix
 
   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);
+  return dJ.contract(Eigen::TensorMap<Eigen::Tensor<const Scalar, 2>>(gInvJacConj.data(), N - 1, N), ip);
 }
 
 std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, const Vector& ζ) {
@@ -98,19 +99,19 @@ std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, con
   Matrix hessZ;
   std::tie(hamiltonian, gradZ, hessZ) = hamGradHess(J, z);
 
-  Matrix g = stereographicMetric(jacobian);
-  // g is Hermitian, which is positive definite.
-  Matrix gInvJac = g.llt().solve(jacobian);
+  Matrix metric = stereographicMetric(jacobian);
 
-  grad = gInvJac * gradZ;
+  grad = metric.llt().solve(jacobian) * gradZ;
 
-  Eigen::Tensor<Scalar, 3> Γ = stereographicChristoffel(ζ, gInvJac);
+  /* This is much slower to calculate than the marginal speedup it offers...
+  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)};
+  std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(2, 0)};
   Eigen::Tensor<Scalar, 2> H2 = Γ.contract(dfT, ip);
+  */
 
-  hess = jacobian * hessZ * jacobian.transpose() - Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size()).transpose();
+  hess = jacobian * hessZ * jacobian.transpose(); // - Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size());
 
   return {hamiltonian, grad, hess};
 }