Better saddle-finding algorithm. Cristoffel symbols still need fixing.

1 files changed, 4 insertions, 3 deletions

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};
 }