From 8209ca60b99594f26f3e9b21ccdbc8695526eb93 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 5 Nov 2021 09:03:36 +0100 Subject: Work on Stokes lines, and new method involving parametric fits. --- dynamics.hpp | 17 +++ langevin.cpp | 129 +++++++-------------- p-spin.hpp | 15 +++ stokes.hpp | 372 +++++++++++++++++++++++++++++++++-------------------------- 4 files changed, 284 insertions(+), 249 deletions(-) diff --git a/dynamics.hpp b/dynamics.hpp index 10b1be2..561714e 100644 --- a/dynamics.hpp +++ b/dynamics.hpp @@ -113,3 +113,20 @@ std::tuple> metropolis(const Tensor& J, const Ve return {E, z}; } + +template +Vector randomSaddle(const Tensor& J, Distribution d, Generator& r, Real ε) { + Vector zSaddle; + bool foundSaddle = false; + + while (!foundSaddle) { + Vector z0 = normalize(randomVector(J.dimension(0), d, r.engine())); + + try { + zSaddle = findSaddle(J, z0, ε); + foundSaddle = true; + } catch (std::exception& e) {} + } + + return zSaddle; +} diff --git a/langevin.cpp b/langevin.cpp index dc7c5cf..fe26332 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -1,8 +1,10 @@ #include +#include #include #include #include +#include "Eigen/src/Eigenvalues/ComplexEigenSolver.h" #include "complex_normal.hpp" #include "p-spin.hpp" #include "dynamics.hpp" @@ -84,11 +86,10 @@ int main(int argc, char* argv[]) { Real T = 1; // temperature Real Rκ = 0; // real part of distribution parameter Real Iκ = 0; // imaginary part of distribution parameter - // simulation parameters Real ε = 1e-15; Real εJ = 1e-5; - Real δ = 1e-2; // threshold for determining saddle + Real δ = 1; // threshold for determining saddle Real Δ = 1e-3; Real γ = 1e-1; // step size unsigned t = 1000; // number of Langevin steps @@ -140,7 +141,7 @@ int main(int argc, char* argv[]) { complex_normal_distribution d(0, 1, 0); - ComplexTensor J = generateCouplings(N, complex_normal_distribution(0, σ, κ), r.engine()); + ComplexTensor J = generateCouplings(N, complex_normal_distribution(0, σ, κ), r.engine()); std::function energyNormGrad = [] (const ComplexTensor& J, const ComplexVector& z) { @@ -149,106 +150,58 @@ int main(int argc, char* argv[]) { return W; }; - - std::function energyThres = [N, κ] - (const ComplexTensor& J, const ComplexVector& z) { - Real a = z.squaredNorm() / (Real)N; - /* - Complex ε = getHamiltonian(J, z) / (Real)N; - Real t = norm(ε) / getThresholdEnergyDensity(p, κ, ε, a); - */ - - Real amin = 1.1547; - Real amax = 1.46806; - - /* - Real tmin = 1; - Real tmax = 1.17514; - */ - - Real E = 0; - - if (a > amax) { - E += pow(a - amax, 2); - } else if (a < amin) { - E += pow(a - amin, 2); - } - - /* - if (t > tmax) { - E += pow(t - tmax, 2); - } else if (t < tmin) { - E += pow(t - tmin, 2); - } - */ - - return E; - }; - - Real largestThreshold = 0; - ComplexVector saddlePastThreshold; - bool foundSaddle = false; - -#pragma omp parallel default(none) shared(largestThreshold, saddlePastThreshold, foundSaddle, M, J, energyThres, d, N, γ, ε, κ) - { - Rng rPrivate; - ComplexVector z = normalize(randomVector(N, d, rPrivate.engine())); - - while (largestThreshold < 1) { // Until we find a saddle past the threshold... - std::tie(std::ignore, z) = metropolis(J, z, energyThres, (Real)0.1, γ, M, d, rPrivate.engine()); - try { - ComplexVector latestSaddle = findSaddle(J, z, ε); - Real threshold = getProportionOfThreshold(κ, J, latestSaddle); - -#pragma omp critical - { - if (threshold > largestThreshold + 1e-6) { - largestThreshold = threshold; - saddlePastThreshold = latestSaddle; - std::cerr << "Found saddle with threshold porportion " << largestThreshold << std::endl;; - } - } - } catch (std::exception& e) {} - } - } - - std::cerr << "Found saddle with energy beyond threshold." << std::endl; + ComplexVector zSaddle = randomSaddle(J, d, r, ε); + std::cerr << "Found saddle." << std::endl; ComplexVector zSaddleNext; - Real closestSaddle = std::numeric_limits::infinity(); - ComplexVector z = saddlePastThreshold; - - while (closestSaddle > δ) { // Until we find two saddles sufficiently close... - std::tie(std::ignore, z) = metropolis(J, saddlePastThreshold, energyNormGrad, T, γ, 1e4, d, r.engine()); + bool foundSaddle = false; + while (!foundSaddle) { + ComplexVector z0 = normalize(zSaddle + δ * randomVector(N, d, r.engine())); try { - zSaddleNext = findSaddle(J, z, ε); - Real saddleDistance = (zSaddleNext - saddlePastThreshold).norm(); - if (saddleDistance < closestSaddle && saddleDistance > 1e-2) { - closestSaddle = saddleDistance; - std::cerr << "Nearby saddle: found saddle at distance " << saddleDistance << std::endl; + zSaddleNext = findSaddle(J, z0, ε); + Real saddleDistance = (zSaddleNext - zSaddle).norm(); + if (saddleDistance / N > 1e-2) { + foundSaddle = true; } } catch (std::exception& e) {} } - std::cerr << "Found sufficiently nearby saddles, perturbing J to equalize Im H." << std::endl; + auto [H1, dH1, ddH1] = hamGradHess(J, zSaddle); + auto [H2, dH2, ddH2] = hamGradHess(J, zSaddleNext); - ComplexVector z1 = saddlePastThreshold; - ComplexVector z2 = zSaddleNext; + Eigen::ComplexEigenSolver ces; + ces.compute(ddH1); - std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-14, ε, r); + Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real()); + std::cerr << (zSaddle - zSaddleNext).norm() / (Real)N << " " << φ << " " << H1 * exp(Complex(0, φ)) << " " << H2 * exp(Complex(0, φ)) << std::endl; + Real smallestNorm = std::numeric_limits::infinity(); + for (Complex z : ces.eigenvalues()) { + if (norm(z) < smallestNorm) { + smallestNorm = norm(z); + } + } + std::cerr << smallestNorm << std::endl; - std::cerr << "Im H is now sufficently close, starting to relax rope." << std::endl; + J = exp(Complex(0, φ)) * J; - std::cerr << "Threshold proportions of saddles are " << getProportionOfThreshold(κ, J, z1) << " and " << getProportionOfThreshold(κ, J, z2) << std::endl; + /* + if (stokesLineTest(J, zSaddle, zSaddleNext, 10, 4)) { + std::cerr << "Found a Stokes line" << std::endl; +// stokesLineTestNew(J, zSaddle, zSaddleNext, 10, 3); + } else { + std::cerr << "Didn't find a Stokes line" << std::endl; + } + */ - Rope stokes(10, z1, z2, J); - for (unsigned i = 0; i < 9; i++) { - stokes.relaxDiscreteGradient(J, 1e6, 0.1, 0); + Cord test(J, zSaddle, zSaddleNext, 5); + test.relax(J, 20, 1, 1e5); - std::cout << stokes.n() << " " << stokes.cost(J) << " " << stokes.length() << std::endl; + std::cout << test.z0.transpose() << std::endl; + std::cout << test.z1.transpose() << std::endl; - stokes = stokes.interpolate(); + for (Vector& g : test.gs) { + std::cout << g.transpose() << std::endl; } return 0; diff --git a/p-spin.hpp b/p-spin.hpp index b2bdf07..e5cd195 100644 --- a/p-spin.hpp +++ b/p-spin.hpp @@ -89,3 +89,18 @@ std::tuple> WdW(const Tensor& J, const Vector +Matrix dzDot(const Vector& z, const Vector& dH) { + Real z² = z.squaredNorm(); + return (dH.conjugate() - (dH.dot(z) / z²) * z.conjugate()) * z.adjoint() / z²; +} + +template +Matrix dzDotConjugate(const Vector& z, const Vector& dH, const Matrix& ddH) { + Real z² = z.squaredNorm(); + return -ddH + (ddH * z.conjugate()) * z.transpose() / z² + + (z.dot(dH) / z²) * ( + Matrix::Identity(ddH.rows(), ddH.cols()) - z.conjugate() * z.transpose() / z² + ); +} diff --git a/stokes.hpp b/stokes.hpp index cbe5437..b2694a9 100644 --- a/stokes.hpp +++ b/stokes.hpp @@ -4,55 +4,14 @@ #include "complex_normal.hpp" #include "dynamics.hpp" +#include + 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) { - Real z² = z.squaredNorm(); - Real z´² = z´.squaredNorm(); - - Vector ż = zDot(z, dH); - Real ż² = ż.squaredNorm(); - - Real Reż·z´ = real(ż.dot(z´)); - - Matrix dż = (dH.conjugate() - (dH.dot(z) / z²) * z.conjugate()) * z.adjoint() / z²; - Matrix dżc = -ddH + (ddH * z.conjugate()) * z.transpose() / z² - + (z.dot(dH) / z²) * ( - Matrix::Identity(ddH.rows(), ddH.cols()) - z.conjugate() * z.transpose() / z² - ); - - Vector dLdz = - ( - dżc * z´ + dż * z´.conjugate() - (dż * ż.conjugate() + dżc * ż) * Reż·z´ / ż² - ) / sqrt(ż² * z´²) / 2; - - Vector ż´ = -(ddH * z´).conjugate() + ((ddH * z´).dot(z) / z²) * z.conjugate() + ( - dH.dot(z) * z´.conjugate() + dH.dot(z´) * z.conjugate() - ( - dH.dot(z) * (z´.dot(z) + z.dot(z´)) / z² - ) * z.conjugate() - ) / z²; - - Real dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´)); - - Vector ddtdLdz´ = - ( - ( - ż´.conjugate() - ( - Reż·z´ * z´´.conjugate() + dReż·z´ * z´.conjugate() - - (Reż·z´ / z´²) * (z´´.dot(z´) + z´.dot(z´´)) * z´.conjugate() - ) / z´² - ) - - 0.5 * ( - (ż.dot(ż´) + ż´.dot(ż)) / ż² + (z´´.dot(z´) + z´.dot(z´´)) / z´² - ) * (ż.conjugate() - Reż·z´ / z´² * z´.conjugate()) - ) / sqrt(ż² * z´²) / 2; - - return dLdz - ddtdLdz´; -} - template class Rope { public: @@ -93,65 +52,10 @@ class Rope { return l; } - template - Real error(const Tensor& J) const { - Scalar H0, HN; - std::tie(H0, std::ignore, std::ignore) = hamGradHess(J, z.front()); - std::tie(HN, std::ignore, std::ignore) = hamGradHess(J, z.back()); - - Real ImH = imag((H0 + HN) / 2.0); - - Real err = 0; - - for (unsigned i = 1; i < z.size() - 1; i++) { - Scalar Hi; - std::tie(Hi, std::ignore, std::ignore) = hamGradHess(J, z[i]); - - err += pow(imag(Hi) - ImH, 2); - } - - return sqrt(err); - } - Vector dz(unsigned i) const { return z[i + 1] - z[i - 1]; } - Vector ddz(unsigned i) const { - return 4.0 * (z[i + 1] + z[i - 1] - 2.0 * z[i]); - } - - template - std::vector> generateGradientδ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]); - - δz[i] = variation(z[i], dz(i), ddz(i), dH, ddH); - } - - for (unsigned i = 1; i < z.size() - 1; i++) { - δz[i] = δz[i].conjugate() - (δz[i].dot(z[i]) / z[i].squaredNorm()) * z[i].conjugate(); -// δz[i] = δz[i] - ((δz[i].conjugate().dot(dz(i))) / dz(i).squaredNorm()) * dz(i).conjugate(); - } - - // We return a δz with average norm of one. - Real mag = 0; - for (unsigned i = 1; i < z.size() - 1; i++) { - mag += δz[i].norm(); - } - - for (unsigned i = 1; i < z.size() - 1; i++) { - δz[i] /= mag / n(); - } - - return δz; - } - template std::vector> generateDiscreteGradientδz(const Tensor& J, Real γ) const { std::vector> δz(z.size()); @@ -185,29 +89,38 @@ class Rope { dC += 0.5 * (ż[i + 1].conjugate() - dz(i + 1).conjugate() * real(ż[i + 1].dot(dz(i + 1))) / dz(i + 1).squaredNorm()) / (ż[i + 1].norm() * dz(i + 1).norm()); } - dC += - γ * (z[i - 1] + z[i + 1]).conjugate(); + dC += γ * (2 * z[i] - z[i - 1] - z[i + 1]).conjugate(); - δz[i] = dC; - } + δz[i] = dC.conjugate(); - Real size = 0; - for (unsigned i = 1; i < z.size() - 1; i++) { - δz[i] = δz[i].conjugate() - (δz[i].dot(z[i]) / z[i].squaredNorm()) * z[i].conjugate(); + δz[i] -= z[i].conjugate() * z[i].conjugate().dot(δz[i]) / z²; } return δz; } - template - std::vector> generateRandomδz(Gen& r) const { - std::vector> δz(z.size()); + void spread() { + Real l = length(); + + Real a = 0; + unsigned pos = 0; + + std::vector> zNew = z; - complex_normal_distribution<> d(0, 1, 0); for (unsigned i = 1; i < z.size() - 1; i++) { - δz[i] = randomVector(z[0].size(), d, r); + Real b = i * l / (z.size() - 1); + + while (b > a) { + pos++; + a += (z[pos] - z[pos - 1]).norm(); + } + + Vector δz = z[pos] - z[pos - 1]; + + zNew[i] = normalize(z[pos] - (a - b) / δz.norm() * δz); } - return δz; + z = zNew; } template @@ -222,8 +135,6 @@ class Rope { rNew.z[i] = normalize(z[i] - (δ * Δl) * δz[i]); } - rNew.spread(); - if (rNew.cost(J, γ) < cost(J, γ)) { break; } else { @@ -235,71 +146,33 @@ class Rope { } } +// rNew.spread(); + z = rNew.z; return δ; } - void spread() { - Real l = length(); - - Real a = 0; - unsigned pos = 0; - - std::vector> zNew = z; - - for (unsigned i = 1; i < z.size() - 1; i++) { - Real b = i * l / (z.size() - 1); - - while (b > a) { - pos++; - a += (z[pos] - z[pos - 1]).norm(); - } - - Vector δz = z[pos] - z[pos - 1]; - - zNew[i] = normalize(z[pos] - (a - b) / δz.norm() * δz); - } - - z = zNew; - } - - template - void relaxGradient(const Tensor& J, unsigned N, Real δ0) { - Real δ = δ0; - try { - for (unsigned i = 0; i < N; i++) { - std::vector> δz = generateGradientδz(J); - δ = 1.1 * perturb(J, δ, δz); - } - } catch (std::exception& e) { - } - } - template void relaxDiscreteGradient(const Tensor& J, unsigned N, Real δ0, Real γ) { Real δ = δ0; try { for (unsigned i = 0; i < N; i++) { std::vector> δz = generateDiscreteGradientδz(J, γ); - δ = 1.1 * perturb(J, δ, δz, γ); + double stepSize = 0; + for (const Vector& v : δz) { + stepSize += v.norm(); + } + if (stepSize / δz.size() < 1e-6) { + break; + } + std::cout << cost(J) << " " << stepSize / δz.size() << std::endl; + δ = 2 * perturb(J, δ, δz, γ); } } catch (std::exception& e) { } } - template - void relaxRandom(const Tensor& J, unsigned N, Real δ0, Gen& r) { - Real δ = δ0; - for (unsigned i = 0; i < N; i++) { - try { - std::vector> δz = generateRandomδz(r); - δ = 1.1 * perturb(J, δ, δz); - } catch (std::exception& e) { - } - } - } - template Real cost(const Tensor& J, Real γ = 0) const { Real c = 0; @@ -334,3 +207,180 @@ class Rope { return r; } }; + +template +bool stokesLineTest(const Tensor& J, const Vector& z1, const Vector& z2, unsigned n0, unsigned steps) { + Rope stokes(n0, z1, z2, J); + + Real oldCost = stokes.cost(J); + + for (unsigned i = 0; i < steps; i++) { + stokes.relaxDiscreteGradient(J, 1e6, 1, pow(2, steps)); + + Real newCost = stokes.cost(J); + + if (newCost > oldCost) { + return false; + } + + oldCost = newCost; + + stokes = stokes.interpolate(); + } + return true; +} + +template +class Cord { +public: + std::vector> gs; + Vector z0; + Vector z1; + + template + Cord(const Tensor& J, const Vector& z2, const Vector& z3, unsigned ng) : gs(ng, Vector::Zero(z2.size())) { + Scalar H2 = getHamiltonian(J, z2); + Scalar H3 = getHamiltonian(J, z3); + + if (real(H2) > real(H3)) { + z0 = z2; + z1 = z3; + } else { + z0 = z3; + z1 = z2; + } + } + + Real gCoeff(unsigned i, Real t) const { + return (1 - t) * t * pow(t, i); + } + + Real dgCoeff(unsigned i, Real t) const { + return (i + 1) * (1 - t) * pow(t, i) - pow(t, i + 1); + } + + Vector f(Real t) const { + Vector z = (1 - t) * z0 + t * z1; + + for (unsigned i = 0; i < gs.size(); i++) { + z += gCoeff(i, t) * gs[i]; + } + + return z; + } + + Vector df(Real t) const { + Vector z = z1 - z0; + + for (unsigned i = 0; i < gs.size(); i++) { + z += dgCoeff(i, t) * gs[i]; + } + + return z; + } + + template + Real cost(const Tensor& J, Real t) const { + Vector z = f(t); + Scalar H; + Vector dH; + std::tie(H, dH, std::ignore) = hamGradHess(J, z); + Vector ż = zDot(z, dH); + Vector dz = df(t); + + return 1 - real(ż.dot(dz)) / ż.norm() / dz.norm(); + } + + template + Real totalCost(const Tensor& J, unsigned nt) const { + Real tc = 0; + + for (unsigned i = 0; i < nt; i++) { + Real t = (i + 1.0) / (nt + 1.0); + tc += cost(J, t); + } + + return tc; + } + + template + std::vector> dgs(const Tensor& J, Real t) const { + Vector z = f(t); + auto [H, dH, ddH] = hamGradHess(J, z); + Vector ż = zDot(z, dH); + Vector dz = df(t); + Matrix dż = dzDot(z, dH); + Matrix dżc = dzDotConjugate(z, dH, ddH); + + std::vector> x; + x.reserve(gs.size()); + + for (unsigned i = 0; i < gs.size(); i++) { + Real fdg = gCoeff(i, t); + Real dfdg = dgCoeff(i, t); + Vector dC = - 0.5 / ż.norm() / dz.norm() * ( + dfdg * ż.conjugate() + fdg * dżc * dz + fdg * dż * dz.conjugate() + - real(dz.dot(ż)) * ( + dfdg * dz.conjugate() / dz.squaredNorm() + + fdg * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm() + ) + ); + + x.push_back(dC.conjugate()); + } + + return x; + } + + template + std::pair relaxStep(const Tensor& J, unsigned nt, Real δ₀) { + std::vector> dgsTot(gs.size(), Vector::Zero(z0.size())); + + for (unsigned i = 0; i < nt; i++) { + Real t = (i + 1.0) / (nt + 1.0); + std::vector> dgsI = dgs(J, t); + + for (unsigned j = 0; j < gs.size(); j++) { + dgsTot[j] += dgsI[j] / nt; + } + } + + Real stepSize = 0; + for (const Vector& dgi : dgsTot) { + stepSize += dgi.squaredNorm(); + } + stepSize = sqrt(stepSize); + + Cord cNew(*this); + + Real δ = δ₀; + Real oldCost = totalCost(J, nt); + Real newCost = std::numeric_limits::infinity(); + + while (newCost > oldCost) { + for (unsigned i = 0; i < gs.size(); i++) { + cNew.gs[i] = gs[i] - δ * dgsTot[i]; + } + + newCost = cNew.totalCost(J, nt); + + δ /= 2; + } + + gs = cNew.gs; + + return {2 * δ, stepSize}; + } + + template + void relax(const Tensor& J, unsigned nt, Real δ₀, unsigned maxSteps) { + Real δ = δ₀; + Real stepSize = std::numeric_limits::infinity(); + unsigned steps = 0; + while (stepSize > 1e-7 && steps < maxSteps) { + std::tie(δ, stepSize) = relaxStep(J, nt, δ); + std::cout << totalCost(J, nt) << " " << δ << " " << stepSize << std::endl; + steps++; + } + } +}; -- cgit v1.2.3-54-g00ecf