From 95df02c90a455c2e539e795acb1921b688e8bc66 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 17 Feb 2021 15:46:00 +0100 Subject: Algorithmic tweaks to rope relaxation now suceeds in finding Stokes lines relatively quickly. --- dynamics.hpp | 9 ++----- langevin.cpp | 77 +++++++++++++++++++++++++++++++++++++----------------------- p-spin.hpp | 5 ++++ stokes.hpp | 73 +++++++++++++++++++++++++++++++++++++++++++------------- 4 files changed, 112 insertions(+), 52 deletions(-) diff --git a/dynamics.hpp b/dynamics.hpp index a27ccb4..22d590a 100644 --- a/dynamics.hpp +++ b/dynamics.hpp @@ -7,16 +7,11 @@ #include "p-spin.hpp" #include "stereographic.hpp" -template -Vector normalize(const Vector& z) { - return z * sqrt((double)z.size() / (Scalar)(z.transpose() * z)); -} - class gradientDescentStallException: public std::exception { virtual const char* what() const throw() { return "Gradient descent stalled."; } -} gradientDescentStall; +}; template std::tuple> gradientDescent(const Tensor& J, const Vector& z0, double ε, double γ0 = 1, double δγ = 2) { @@ -41,7 +36,7 @@ std::tuple> gradientDescent(const Tensor& J, c } if (γ < 1e-50) { - throw gradientDescentStall; + throw gradientDescentStallException(); } } diff --git a/langevin.cpp b/langevin.cpp index 5d9acc9..e1464ec 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -1,5 +1,7 @@ #include #include +#include +#include #include "complex_normal.hpp" #include "p-spin.hpp" @@ -19,6 +21,20 @@ using ComplexTensor = Tensor; using Rng = randutils::random_generator; +std::list> saddlesCloserThan(const std::unordered_map& saddles, double δ) { + std::list> pairs; + + for (auto it1 = saddles.begin(); it1 != saddles.end(); it1++) { + for (auto it2 = std::next(it1); it2 != saddles.end(); it2++) { + if ((it1->second - it2->second).norm() < δ) { + pairs.push_back({it1->second, it2->second}); + } + } + } + + return pairs; +} + int main(int argc, char* argv[]) { // model parameters unsigned N = 10; // number of spins @@ -82,10 +98,7 @@ int main(int argc, char* argv[]) { complex_normal_distribution<> d(0, 1, 0); ComplexTensor J = generateCouplings(N, complex_normal_distribution<>(0, σ, κ), r.engine()); - ComplexVector z0 = normalize(randomVector(N, d, r.engine())); - - ComplexVector zSaddle = findSaddle(J, z0, ε); - ComplexVector z = zSaddle; + ComplexVector z = normalize(randomVector(N, d, r.engine())); std::function energyNormGrad = [] (const ComplexTensor& J, const ComplexVector& z) { @@ -94,49 +107,55 @@ int main(int argc, char* argv[]) { return W; }; - ComplexVector zSaddlePrev = ComplexVector::Zero(N); + std::unordered_map saddles; + std::list> nearbySaddles; - while (δ < (zSaddle - zSaddlePrev).norm()) { // Until we find two saddles sufficiently close... + while (true) { // Until we find two saddles sufficiently close... std::tie(std::ignore, z) = metropolis(J, z, energyNormGrad, T, γ, M, d, r.engine()); try { ComplexVector zSaddleNext = findSaddle(J, z, ε); - if (Δ < (zSaddleNext - zSaddle).norm()) { // Ensure we are finding distinct saddles. - zSaddlePrev = zSaddle; - zSaddle = zSaddleNext; + uint64_t saddleKey = round(1e2 * real(zSaddleNext(0))); + auto got = saddles.find(saddleKey); + + if (got == saddles.end()) { + saddles[saddleKey] = zSaddleNext; + nearbySaddles = saddlesCloserThan(saddles, δ); + if (nearbySaddles.size() > 0) { + break; + } + std::cerr << "Found " << saddles.size() << " distinct saddles." << std::endl; } } catch (std::exception& e) { - std::cerr << "Could not find a saddle: " << e.what() << std::endl; +// std::cerr << "Could not find a saddle: " << e.what() << std::endl; } - std::cerr << "Current saddles are " << (zSaddle - zSaddlePrev).norm() << " apart." << std::endl; } - Rope stokest(20, zSaddle, zSaddlePrev); - - stokest.relax(J, 10000, 0.0001); - std::cerr << "Found sufficiently nearby saddles, perturbing J to equalize Im H." << std::endl; complex_normal_distribution<> dJ(0, εJ * σ, 0); std::function)> perturbJ = - [&dJ, &r] (ComplexTensor& JJ, std::array ind) { + [&εJ, &dJ, &r] (ComplexTensor& JJ, std::array ind) { Complex Ji = getJ(JJ, ind); - setJ(JJ, ind, Ji + dJ(r.engine())); + setJ(JJ, ind, Ji + εJ * dJ(r.engine())); }; ComplexTensor Jp = J; - ComplexVector z1, z2; Complex H1, H2; - double prevdist = 100; + ComplexVector z1, z2; + std::tie(H1, std::ignore, std::ignore) = hamGradHess(Jp, nearbySaddles.front()[0]); + std::tie(H2, std::ignore, std::ignore) = hamGradHess(Jp, nearbySaddles.front()[1]); + double prevdist = abs(imag(H1-H2)); + εJ = 1e4 * prevdist; while (true) { ComplexTensor Jpp = Jp; iterateOver(Jpp, perturbJ); try { - z1 = findSaddle(Jpp, zSaddle, ε); - z2 = findSaddle(Jpp, zSaddlePrev, ε); + z1 = findSaddle(Jpp, nearbySaddles.front()[0], ε); + z2 = findSaddle(Jpp, nearbySaddles.front()[1], ε); double dist = (z1 - z2).norm(); @@ -146,10 +165,10 @@ int main(int argc, char* argv[]) { if (abs(imag(H1 - H2)) < prevdist && dist > 1e-2) { Jp = Jpp; prevdist = abs(imag(H1 - H2)); + εJ = 1e4 * prevdist; - std::cout << abs(imag(H1 - H2)) << " " << dist << std::endl; - if (abs(imag(H1 - H2)) < 1e-8 && dist > 1e-2) { - std::cout << "Found distinct critical points with sufficiently similar Im H." << std::endl; + if (abs(imag(H1 - H2)) < 1e-10 && dist > 1e-2) { + std::cerr << "Found distinct critical points with sufficiently similar Im H." << std::endl; break; } } @@ -161,29 +180,29 @@ int main(int argc, char* argv[]) { } } - Rope stokes(20, z1, z2); + Rope stokes(10, z1, z2); std::cout << stokes.cost(Jp) << std::endl; - stokes.relax(Jp, 10000, 0.0001); + stokes.relax(Jp, 10000, 1e-4); std::cout << stokes.cost(Jp) << std::endl; Rope stokes2 = stokes.interpolate(); - stokes2.relax(Jp, 10000, 0.001); + stokes2.relax(Jp, 10000, 1e-4); std::cout << stokes2.cost(Jp) << std::endl; stokes2 = stokes2.interpolate(); - stokes2.relax(Jp, 10000, 0.001); + stokes2.relax(Jp, 10000, 1e-4); std::cout << stokes2.cost(Jp) << std::endl; stokes2 = stokes2.interpolate(); - stokes2.relax(Jp, 10000, 0.001); + stokes2.relax(Jp, 10000, 1e-4); std::cout << stokes2.cost(Jp) << std::endl; diff --git a/p-spin.hpp b/p-spin.hpp index 13fbd98..480b3ca 100644 --- a/p-spin.hpp +++ b/p-spin.hpp @@ -20,6 +20,11 @@ std::tuple, Matrix> hamGradHess(const Tensor +Vector normalize(const Vector& z) { + return z * sqrt((double)z.size() / (Scalar)(z.transpose() * z)); +} + template Vector project(const Vector& z, const Vector& x) { Scalar xz = x.transpose() * z; diff --git a/stokes.hpp b/stokes.hpp index 6a00637..95bb9c4 100644 --- a/stokes.hpp +++ b/stokes.hpp @@ -1,6 +1,4 @@ #include "p-spin.hpp" -#include -#include "dynamics.hpp" template Vector zDot(const Vector& z, const Vector& dH) { @@ -13,6 +11,12 @@ double segmentCost(const Vector& z, const Vector& dz, const Vect return 1.0 - pow(real(zD.dot(dz)), 2) / zD.squaredNorm() / dz.squaredNorm(); } +class ropeRelaxationStallException: public std::exception { + virtual const char* what() const throw() { + return "Gradient descent stalled."; + } +}; + template Vector variation(const Vector& z, const Vector& z´, const Vector& z´´, const Vector& dH, const Matrix& ddH) { double z² = z.squaredNorm(); @@ -78,29 +82,58 @@ class Rope { } template - void perturb(const Tensor& J, double δ) { + std::vector> δz(const Tensor& J) const { + std::vector> δz(z.size()); + +#pragma omp parallel for for (unsigned i = 1; i < z.size() - 1; i++) { Vector dH; Matrix ddH; std::tie(std::ignore, dH, ddH) = hamGradHess(J, z[i]); - Vector δz = variation(z[i], this->dz(i), this->ddz(i), dH, ddH); + δz[i] = variation(z[i], this->dz(i), this->ddz(i), dH, ddH); + } - z[i] -= δ * δz.conjugate(); - z[i] = normalize(z[i]); + return δz; + } + + template + double perturb(const Tensor& J, double δ0) { + std::vector> Δz = this->δz(J); + + Rope rNew = *this; + double δ = δ0; + + while (rNew.cost(J) >= this->cost(J)) { + for (unsigned i = 1; i < z.size() - 1; i++) { + rNew.z[i] = z[i] - δ * Δz[i].conjugate(); + rNew.z[i] = normalize(rNew.z[i]); + } + + δ /= 2; + + if (δ < 1e-30) { + throw ropeRelaxationStallException(); + } } + + z = rNew.z; + + return δ; } void spread() { double l = 0; - for (unsigned i= 0; i < z.size() - 1; i++) { + for (unsigned i = 0; i < z.size() - 1; i++) { l += (z[i + 1] - z[i]).norm(); } double a = 0; unsigned pos = 0; + std::vector> zNew = z; + for (unsigned i = 1; i < z.size() - 1; i++) { double b = i * l / (z.size() - 1); @@ -111,15 +144,19 @@ class Rope { Vector δz = z[pos] - z[pos - 1]; - z[i] = z[pos] - (a - b) / δz.norm() * δz; - z[i] = normalize(z[i]); + zNew[i] = z[pos] - (a - b) / δz.norm() * δz; + zNew[i] = normalize(zNew[i]); } + + z = zNew; } template - void step(const Tensor& J, double δ) { - this->perturb(J, δ); + double step(const Tensor& J, double δ) { + double δNew = this->perturb(J, δ); this->spread(); + + return δNew; } template @@ -137,21 +174,25 @@ class Rope { } template - void relax(const Tensor& J, unsigned N, double δ) { - for (unsigned i = 0; i < N; i++) { - this->step(J, δ); + void relax(const Tensor& J, unsigned N, double δ0) { + double δ = δ0; + try { + for (unsigned i = 0; i < N; i++) { + δ = 2 * this->step(J, δ); + } + } catch (std::exception& e) { } } Rope interpolate() const { - Rope r(2 * z.size() - 1, z[0], z[z.size() - 1]); + Rope r(2 * (z.size() - 2) + 1, z[0], z[z.size() - 1]); for (unsigned i = 0; i < z.size(); i++) { r.z[2 * i] = z[i]; } for (unsigned i = 0; i < z.size() - 1; i++) { - r.z[2 * i + 1] = (z[i] + z[i + 1]) / 2.0; + r.z[2 * i + 1] = normalize(((z[i] + z[i + 1]) / 2.0).eval()); } return r; -- cgit v1.2.3-70-g09d2