diff options
| -rw-r--r-- | dynamics.hpp | 9 | ||||
| -rw-r--r-- | langevin.cpp | 77 | ||||
| -rw-r--r-- | p-spin.hpp | 5 | ||||
| -rw-r--r-- | 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 <class Scalar> -Vector<Scalar> normalize(const Vector<Scalar>& 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 <class Scalar, int p> std::tuple<double, Vector<Scalar>> gradientDescent(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, double ε, double γ0 = 1, double δγ = 2) { @@ -41,7 +36,7 @@ std::tuple<double, Vector<Scalar>> gradientDescent(const Tensor<Scalar, p>& 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 <getopt.h> #include <stereographic.hpp> +#include <unordered_map> +#include <list> #include "complex_normal.hpp" #include "p-spin.hpp" @@ -19,6 +21,20 @@ using ComplexTensor = Tensor<Complex, p>; using Rng = randutils::random_generator<pcg32>; +std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_map<uint64_t, ComplexVector>& saddles, double δ) { + std::list<std::array<ComplexVector, 2>> 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<Complex, PSPIN_P>(N, complex_normal_distribution<>(0, σ, κ), r.engine()); - ComplexVector z0 = normalize(randomVector<Complex>(N, d, r.engine())); - - ComplexVector zSaddle = findSaddle(J, z0, ε); - ComplexVector z = zSaddle; + ComplexVector z = normalize(randomVector<Complex>(N, d, r.engine())); std::function<double(const ComplexTensor&, const ComplexVector&)> 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<uint64_t, ComplexVector> saddles; + std::list<std::array<ComplexVector, 2>> 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<Complex> 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<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ = - [&dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) { + [&εJ, &dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) { Complex Ji = getJ<Complex, p>(JJ, ind); - setJ<Complex, p>(JJ, ind, Ji + dJ(r.engine())); + setJ<Complex, p>(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<Complex, p>(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<Complex> stokes(20, z1, z2); + Rope<Complex> 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<Complex> 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; @@ -21,6 +21,11 @@ std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>> hamGradHess(const Tensor<Scal } template <class Scalar> +Vector<Scalar> normalize(const Vector<Scalar>& z) { + return z * sqrt((double)z.size() / (Scalar)(z.transpose() * z)); +} + +template <class Scalar> Vector<Scalar> project(const Vector<Scalar>& z, const Vector<Scalar>& x) { Scalar xz = x.transpose() * z; return x - (xz / z.squaredNorm()) * z.conjugate(); @@ -1,6 +1,4 @@ #include "p-spin.hpp" -#include <iostream> -#include "dynamics.hpp" template <class Scalar> Vector<Scalar> zDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) { @@ -13,6 +11,12 @@ double segmentCost(const Vector<Scalar>& z, const Vector<Scalar>& 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 <class Scalar> Vector<Scalar> variation(const Vector<Scalar>& z, const Vector<Scalar>& z´, const Vector<Scalar>& z´´, const Vector<Scalar>& dH, const Matrix<Scalar>& ddH) { double z² = z.squaredNorm(); @@ -78,29 +82,58 @@ class Rope { } template <int p> - void perturb(const Tensor<Scalar, p>& J, double δ) { + std::vector<Vector<Scalar>> δz(const Tensor<Scalar, p>& J) const { + std::vector<Vector<Scalar>> δz(z.size()); + +#pragma omp parallel for for (unsigned i = 1; i < z.size() - 1; i++) { Vector<Scalar> dH; Matrix<Scalar> ddH; std::tie(std::ignore, dH, ddH) = hamGradHess(J, z[i]); - Vector<Scalar> δ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 <int p> + double perturb(const Tensor<Scalar, p>& J, double δ0) { + std::vector<Vector<Scalar>> Δ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<Vector<Scalar>> 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<Scalar> δ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 <int p> - void step(const Tensor<Scalar, p>& J, double δ) { - this->perturb(J, δ); + double step(const Tensor<Scalar, p>& J, double δ) { + double δNew = this->perturb(J, δ); this->spread(); + + return δNew; } template <int p> @@ -137,21 +174,25 @@ class Rope { } template <int p> - void relax(const Tensor<Scalar, p>& J, unsigned N, double δ) { - for (unsigned i = 0; i < N; i++) { - this->step(J, δ); + void relax(const Tensor<Scalar, p>& J, unsigned N, double δ0) { + double δ = δ0; + try { + for (unsigned i = 0; i < N; i++) { + δ = 2 * this->step(J, δ); + } + } catch (std::exception& e) { } } Rope<Scalar> interpolate() const { - Rope<Scalar> r(2 * z.size() - 1, z[0], z[z.size() - 1]); + Rope<Scalar> 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; |