From 752d0b61595166aafa93b48dd6c209381ac00f02 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 6 Jan 2021 17:36:28 +0100 Subject: Got Newton's method working well. --- langevin.cpp | 85 ++++++++++++++++++++++++++++++++++++++---------------------- 1 file changed, 54 insertions(+), 31 deletions(-) (limited to 'langevin.cpp') diff --git a/langevin.cpp b/langevin.cpp index acea609..49efab4 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -8,6 +8,7 @@ #include "randutils/randutils.hpp" #include +#include #include "complex_normal.hpp" #include "p-spin.hpp" @@ -15,6 +16,10 @@ using Rng = randutils::random_generator; +Vector normalize(const Vector& z) { + return z * sqrt(z.size()) / sqrt((Scalar)(z.transpose() * z)); +} + template Vector randomVector(unsigned N, Distribution d, Generator& r) { Vector z(N); @@ -36,19 +41,18 @@ std::tuple gradientDescent(const Tensor& J, const Vector& z0, do double γ = γ0; while (W > ε) { - Vector zNew = z - γ * dW.conjugate(); - zNew *= sqrt(z.size()) / sqrt((Scalar)(zNew.transpose() * zNew)); + Vector zNew = normalize(z - γ * dW.conjugate()); double WNew; Vector dWNew; std::tie(WNew, dWNew) = WdW(J, zNew); - if (WNew < W) { + if (WNew < W) { // If the step lowered the objective, accept it! z = zNew; W = WNew; dW = dWNew; γ = γ0; - } else { + } else { // Otherwise, shrink the step and try again. γ *= 0.5; } @@ -70,6 +74,9 @@ Vector findSaddle(const Tensor& J, const Vector& z0, double ε, double δW, doub Matrix ddH; std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ); + unsigned steps = 0; + unsigned gradSteps = 0; + while (W > ε) { // ddH is complex symmetric, which is (almost always) invertible Vector dζ = ddH.partialPivLu().solve(dH); @@ -85,37 +92,42 @@ Vector findSaddle(const Tensor& J, const Vector& z0, double ε, double δW, doub Vector z; std::tie(W, z) = gradientDescent(J, stereographicToEuclidean(ζ), γ0, W / δW); ζ = euclideanToStereographic(z); + gradSteps++; } std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ); - std::cout << W << std::endl; + steps++; + + if (steps % 100 == 0) { + std::cerr << steps << " minimization steps, W is " << W << " " << gradSteps << "." << std::endl; + } } return stereographicToEuclidean(ζ); } -Vector langevin(const Tensor& J, const Vector& z0, double T, double γ0, - std::function quit, Rng& r) { +std::tuple langevin(const Tensor& J, const Vector& z0, double T, double γ, unsigned N, Rng& r) { Vector z = z0; double W; - Vector dW; - std::tie(W, dW) = WdW(J, z); + std::tie(W, std::ignore) = WdW(J, z); - unsigned steps = 0; - complex_normal_distribution<> d(0, sqrt(T), 0); + complex_normal_distribution<> d(0, γ, 0); - while (!quit(W, steps)) { - double γ = pow(r.variate(0, γ0), 2); - Vector η = randomVector(z.size(), d, r.engine()); + for (unsigned i = 0; i < N; i++) { + Vector dz = randomVector(z.size(), d, r.engine()); + Vector zNew = normalize(z + dz); - z -= γ * (dW.conjugate() / pow((double)z.size(), 2) + η); - z *= sqrt(z.size()) / sqrt((Scalar)(z.transpose() * z)); + double WNew; + std::tie(WNew, std::ignore) = WdW(J, zNew); - std::tie(W, dW) = WdW(J, z); + if (exp((W - WNew) / T) > r.uniform(0.0, 1.0)) { + z = zNew; + W = WNew; + } } - return z; + return {W, z}; } int main(int argc, char* argv[]) { @@ -130,10 +142,12 @@ int main(int argc, char* argv[]) { double δ = 1e-2; // threshold for determining saddle double γ = 1e-2; // step size unsigned t = 1000; // number of Langevin steps + unsigned M = 100; + unsigned n = 100; int opt; - while ((opt = getopt(argc, argv, "N:T:e:r:i:g:t:E:")) != -1) { + while ((opt = getopt(argc, argv, "N:M:n:T:e:r:i:g:t:E:")) != -1) { switch (opt) { case 'N': N = (unsigned)atof(optarg); @@ -158,6 +172,12 @@ int main(int argc, char* argv[]) { case 'i': Iκ = atof(optarg); break; + case 'n': + n = (unsigned)atof(optarg); + break; + case 'M': + M = (unsigned)atof(optarg); + break; default: exit(1); } @@ -169,24 +189,27 @@ int main(int argc, char* argv[]) { Rng r; Tensor J = generateCouplings(N, complex_normal_distribution<>(0, σ, κ), r.engine()); - Vector z = randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine()); - z *= sqrt(N) / sqrt((Scalar)(z.transpose() * z)); // Normalize. + Vector z0 = randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine()); + z0 *= sqrt(N) / sqrt((Scalar)(z0.transpose() * z0)); // Normalize. std::function f = [δ](double W, unsigned) { std::cout << W << std::endl; return W < δ; }; - //Vector zm = langevin(J, z, T, γ, f, r); - Vector zm = findSaddle(J, z, ε, δ, γ); - - Scalar H; - Vector dH; - - std::tie(H, dH, std::ignore) = hamGradHess(J, zm); + Vector zSaddle = findSaddle(J, z0, ε, δ, γ); - Vector constraint = dH - ((double)p * H / (double)N) * zm; + double W; + Vector z = zSaddle; + for (unsigned i = 0; i < n; i++) { + std::tie(W, z) = langevin(J, z, T, γ, M, r); + Vector zNewSaddle = findSaddle(J, z, ε, δ, γ); + Scalar H; + Matrix ddH; + std::tie(H, std::ignore, ddH) = hamGradHess(J, zNewSaddle); + Eigen::SelfAdjointEigenSolver es(ddH.adjoint() * ddH); + std::cout << (zNewSaddle - zSaddle).norm() << " " << real(H) << " " << imag(H) << " " << es.eigenvalues().transpose() << std::endl; + } - std::cout << constraint << std::endl; - std::cout << H / (double)N << std::endl; + return 0; } -- cgit v1.2.3-70-g09d2