diff options
author | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2021-01-05 18:11:21 +0100 |
---|---|---|
committer | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2021-01-05 18:11:21 +0100 |
commit | 5ee6815f0734b2089c5b4c068cc21f2983bdba24 (patch) | |
tree | 179e2ec28f57662ce91491d79799270f2727b144 /stereographic.hpp | |
parent | 7c3e71970b6f2d48530bc2ab4fc0f2932522b98b (diff) | |
download | code-5ee6815f0734b2089c5b4c068cc21f2983bdba24.tar.gz code-5ee6815f0734b2089c5b4c068cc21f2983bdba24.tar.bz2 code-5ee6815f0734b2089c5b4c068cc21f2983bdba24.zip |
A lot of work, and fixed a huge bug regarding the meaning of .dot in the Eigen library for complex vectors.
Diffstat (limited to 'stereographic.hpp')
-rw-r--r-- | stereographic.hpp | 67 |
1 files changed, 45 insertions, 22 deletions
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}; } |