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.

---
 stereographic.hpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

(limited to 'stereographic.hpp')

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-54-g00ecf