diff options
| -rw-r--r-- | langevin.cpp | 38 | ||||
| -rw-r--r-- | stereographic.hpp | 95 |
2 files changed, 37 insertions, 96 deletions
diff --git a/langevin.cpp b/langevin.cpp index 49efab4..e5ed0b8 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -69,41 +69,36 @@ Vector findSaddle(const Tensor& J, const Vector& z0, double ε, double δW, doub double W; std::tie(W, std::ignore) = WdW(J, z0); - Vector ζ = euclideanToStereographic(z0); + Vector z = z0; + Vector ζ = euclideanToStereographic(z); + Vector dH; Matrix ddH; - std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ); - - unsigned steps = 0; - unsigned gradSteps = 0; + std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ, z); while (W > ε) { - // ddH is complex symmetric, which is (almost always) invertible + // ddH is complex symmetric, which is (almost always) invertible, so a + // partial pivot LU decomposition can be used. Vector dζ = ddH.partialPivLu().solve(dH); Vector ζNew = ζ - dζ; + Vector zNew = stereographicToEuclidean(ζNew); double WNew; - std::tie(WNew, std::ignore) = WdW(J, stereographicToEuclidean(ζNew)); + std::tie(WNew, std::ignore) = WdW(J, zNew); if (WNew < W) { // If Newton's step lowered the objective, accept it! ζ = ζNew; + z = zNew; W = WNew; } else { // Otherwise, do gradient descent until W is a factor δW smaller. - Vector z; - std::tie(W, z) = gradientDescent(J, stereographicToEuclidean(ζ), γ0, W / δW); + std::tie(W, z) = gradientDescent(J, z, γ0, W / δW); ζ = euclideanToStereographic(z); - gradSteps++; } - std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ); - steps++; - - if (steps % 100 == 0) { - std::cerr << steps << " minimization steps, W is " << W << " " << gradSteps << "." << std::endl; - } + std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ, z); } - return stereographicToEuclidean(ζ); + return z; } std::tuple<double, Vector> langevin(const Tensor& J, const Vector& z0, double T, double γ, unsigned N, Rng& r) { @@ -189,13 +184,7 @@ int main(int argc, char* argv[]) { Rng r; Tensor J = generateCouplings<Scalar, PSPIN_P>(N, complex_normal_distribution<>(0, σ, κ), r.engine()); - Vector z0 = randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine()); - z0 *= sqrt(N) / sqrt((Scalar)(z0.transpose() * z0)); // Normalize. - - std::function<bool(double, unsigned)> f = [δ](double W, unsigned) { - std::cout << W << std::endl; - return W < δ; - }; + Vector z0 = normalize(randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine())); Vector zSaddle = findSaddle(J, z0, ε, δ, γ); @@ -209,6 +198,7 @@ int main(int argc, char* argv[]) { std::tie(H, std::ignore, ddH) = hamGradHess(J, zNewSaddle); Eigen::SelfAdjointEigenSolver<Matrix> es(ddH.adjoint() * ddH); std::cout << (zNewSaddle - zSaddle).norm() << " " << real(H) << " " << imag(H) << " " << es.eigenvalues().transpose() << std::endl; + std::cerr << M * (i+1) << " steps taken to move " << (zNewSaddle - zSaddle).norm() << ", saddle information saved." << std::endl; } return 0; diff --git a/stereographic.hpp b/stereographic.hpp index a3f2cc6..8313f25 100644 --- a/stereographic.hpp +++ b/stereographic.hpp @@ -1,22 +1,21 @@ #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 = ζ.transpose() * ζ; Scalar b = 2 * sqrt(N) / (1.0 + a); for (unsigned i = 0; i < N - 1; i++) { - z(i) = b * ζ(i); + z(i) = ζ(i); } - z(N - 1) = b * (a - 1.0) / 2.0; + z(N - 1) = (a - 1.0) / 2.0; - return z; + return b * z; } Vector euclideanToStereographic(const Vector& z) { @@ -24,94 +23,46 @@ Vector euclideanToStereographic(const Vector& z) { Vector ζ(N - 1); for (unsigned i = 0; i < N - 1; i++) { - ζ(i) = z(i) / (sqrt(N) - z(N - 1)); + ζ(i) = z(i); } - return ζ; + return ζ / (sqrt(N) - z(N - 1)); } Matrix stereographicJacobian(const Vector& ζ) { - unsigned N = ζ.size() + 1; - Matrix J(N - 1, N); - - Scalar b = ζ.transpose() * ζ; - - for (unsigned i = 0; i < N - 1; i++) { - for (unsigned j = 0; j < N - 1; j++) { - J(i, j) = - 4 * sqrt(N) * ζ(i) * ζ(j) / pow(1.0 + b, 2); - if (i == j) { - J(i, j) += 2 * sqrt(N) * (1.0 + b) / pow(1.0 + b, 2); - } - } - - J(i, N - 1) = 4.0 * sqrt(N) * ζ(i) / pow(1.0 + b, 2); - } - - return J; -} - -Matrix stereographicMetric(const Matrix& J) { - return J * J.adjoint(); -} - -// 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); + unsigned N = ζ.size(); + Matrix J(N, N + 1); Scalar b = 1.0 + (Scalar)(ζ.transpose() * ζ); - 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(b, 3); - if (i == j) { - dJ(i, j, k) -= 4 * sqrt(N) * ζ(k) / pow(b, 2); - } - if (i == k) { - dJ(i, j, k) -= 4 * sqrt(N) * ζ(j) / pow(b, 2); - } - if (j == k) { - dJ(i, j, k) -= 4 * sqrt(N) * ζ(i) / pow(b, 2); - } - } - dJ(i, j, N - 1) = - 16 * sqrt(N) * ζ(i) * ζ(j) / pow(b, 3); + for (unsigned i = 0; i < N; i++) { + for (unsigned j = 0; j < N; j++) { + J(i, j) = - ζ(i) * ζ(j); + if (i == j) { - dJ(i, j, N - 1) += 4 * sqrt(N) * ζ(i) / pow(b, 2); + J(i, j) += 0.5 * b; } } - } - std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(2, 1)}; + J(i, N) = ζ(i); + } - return dJ.contract(Eigen::TensorMap<Eigen::Tensor<const Scalar, 2>>(gInvJacConj.data(), N - 1, N), ip); + return 4 * sqrt(N + 1) * J / pow(b, 2); } -std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, const Vector& ζ) { - Vector grad; - Matrix hess; - - Matrix jacobian = stereographicJacobian(ζ); - Vector z = stereographicToEuclidean(ζ); - - std::complex<double> hamiltonian; +std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, const Vector& ζ, const Vector& z) { + Scalar hamiltonian; Vector gradZ; Matrix hessZ; std::tie(hamiltonian, gradZ, hessZ) = hamGradHess(J, z); - Matrix metric = stereographicMetric(jacobian); - - grad = metric.llt().solve(jacobian) * gradZ; + Matrix jacobian = stereographicJacobian(ζ); - /* 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>(2, 0)}; - Eigen::Tensor<Scalar, 2> H2 = Γ.contract(dfT, ip); - */ + Matrix metric = jacobian * jacobian.adjoint(); - hess = jacobian * hessZ * jacobian.transpose(); // - Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size()); + // The metric is Hermitian and positive definite, so a Cholesky decomposition can be used. + Vector grad = metric.llt().solve(jacobian) * gradZ; + Matrix hess = jacobian * hessZ * jacobian.transpose(); return {hamiltonian, grad, hess}; } |