diff options
| -rw-r--r-- | langevin.cpp | 82 | ||||
| -rw-r--r-- | p-spin.hpp | 6 | ||||
| -rw-r--r-- | stereographic.hpp | 67 |
3 files changed, 114 insertions, 41 deletions
diff --git a/langevin.cpp b/langevin.cpp index 9da4ef6..3bc12ed 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -9,6 +9,7 @@ #include "complex_normal.hpp" #include "p-spin.hpp" +#include "stereographic.hpp" using Rng = randutils::random_generator<pcg32>; @@ -23,6 +24,57 @@ Vector randomVector(unsigned N, Distribution d, Generator& r) { return z; } +Vector findSaddle(const Tensor& J, const Vector& z0, double γ0, double δ, double ε, Rng& r) { + Vector z = z0; + + double W; + Vector dW; + std::tie(W, dW) = WdW(J, z); + + while (W > δ) { + double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2); + + Vector zNew = z - γ * dW.conjugate(); + zNew *= sqrt(z.size()) / sqrt((Scalar)(zNew.transpose() * zNew)); + + double WNew; + Vector dWNew; + std::tie(WNew, dWNew) = WdW(J, zNew); + + if (WNew < W) { + z = zNew; + W = WNew; + dW = dWNew; + std::cout << W << std::endl; + } + } + + Vector ζ = euclideanToStereographic(z); + Vector dH; + Matrix ddH; + std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ); + + while (W > ε) { + double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2); + + Vector ζNew = ζ - ddH.ldlt().solve(dH); + + double WNew; + Vector dWNew; + std::tie(WNew, dWNew) = WdW(J, stereographicToEuclidean(ζNew)); + + if (WNew < W) { + ζ = ζNew; + W = WNew; + dW = dWNew; + std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ); + std::cout << WNew << std::endl; + } + } + + return stereographicToEuclidean(ζ); +} + Vector langevin(const Tensor& J, const Vector& z0, double T, double γ0, std::function<bool(double, unsigned)> quit, Rng& r) { Vector z = z0; @@ -32,24 +84,16 @@ Vector langevin(const Tensor& J, const Vector& z0, double T, double γ0, std::tie(W, dW) = WdW(J, z); unsigned steps = 0; - complex_normal_distribution<> d(0, T, 0); + complex_normal_distribution<> d(0, sqrt(T), 0); while (!quit(W, steps)) { double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2); Vector η = randomVector(z.size(), d, r.engine()); - Vector zNext = z - γ * dW + η; - zNext *= sqrt(zNext.size()) / sqrt(zNext.dot(zNext)); + z -= γ * (dW.conjugate() / pow((double)z.size(), 2) + η); + z *= sqrt(z.size()) / sqrt((Scalar)(z.transpose() * z)); - double WNext; - Vector dWNext; - std::tie(WNext, dWNext) = WdW(J, zNext); - - if (exp((W - WNext) / T) > r.uniform(0.0, 1.0)) { - z = zNext; - W = WNext; - dW = dWNext; - } + std::tie(W, dW) = WdW(J, z); } return z; @@ -63,13 +107,14 @@ int main(int argc, char* argv[]) { double Iκ = 0; // imaginary part of distribution parameter // simulation parameters + double ε = 1e-4; double δ = 1e-2; // threshold for determining saddle double γ = 1e-2; // step size unsigned t = 1000; // number of Langevin steps int opt; - while ((opt = getopt(argc, argv, "N:T:e:r:i:g:t:")) != -1) { + while ((opt = getopt(argc, argv, "N:T:e:r:i:g:t:E:")) != -1) { switch (opt) { case 'N': N = (unsigned)atof(optarg); @@ -83,6 +128,9 @@ int main(int argc, char* argv[]) { case 'e': δ = atof(optarg); break; + case 'E': + ε = atof(optarg); + break; case 'g': γ = atof(optarg); case 'r': @@ -103,14 +151,15 @@ int main(int argc, char* argv[]) { Tensor J = generateCouplings<Scalar, PSPIN_P>(N, complex_normal_distribution<>(0, σ, κ), r.engine()); Vector z = randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine()); - z *= sqrt(N) / sqrt(z.dot(z)); // Normalize. + z *= sqrt(N) / sqrt((Scalar)(z.transpose() * z)); // Normalize. - std::function<bool(double, unsigned)> findSaddle = [δ](double W, unsigned) { + std::function<bool(double, unsigned)> f = [δ](double W, unsigned) { std::cout << W << std::endl; return W < δ; }; - Vector zm = langevin(J, z, T, γ, findSaddle, r); + //Vector zm = langevin(J, z, T, γ, f, r); + Vector zm = findSaddle(J, z, γ, δ, ε, r); Scalar H; Vector dH; @@ -119,5 +168,6 @@ int main(int argc, char* argv[]) { Vector constraint = dH - ((double)p * H / (double)N) * zm; + std::cout << constraint << std::endl; std::cout << H / (double)N << std::endl; } @@ -24,7 +24,7 @@ std::tuple<Scalar, Vector, Matrix> hamGradHess(const Tensor& J, const Vector& z) Matrix hessian = ((p - 1) * p / f) * Jz; Vector gradient = (p / f) * Jzz; - Scalar hamiltonian = (1 / f) * Jzz.dot(z); + Scalar hamiltonian = (1 / f) * Jzz.transpose() * z; return {hamiltonian, gradient, hessian}; } @@ -34,10 +34,10 @@ std::tuple<double, Vector> WdW(const Tensor& J, const Vector& z) { Matrix hessian; std::tie(std::ignore, gradient, hessian) = hamGradHess(J, z); - Vector projectedGradient = gradient - (gradient.dot(z) / (double)z.size()) * z; + Vector projectedGradient = (gradient - ((Scalar)(gradient.transpose() * z) / (double)z.size()) * z).conjugate(); double W = projectedGradient.cwiseAbs2().sum(); - Vector dW = hessian.conjugate() * projectedGradient; + Vector dW = hessian * projectedGradient - ((z.transpose() * gradient) * projectedGradient + (z.transpose() * projectedGradient) * (gradient + hessian * z)) / (double)z.size(); return {W, dW}; } 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}; } |