diff options
| -rw-r--r-- | dynamics.hpp | 17 | ||||
| -rw-r--r-- | langevin.cpp | 129 | ||||
| -rw-r--r-- | p-spin.hpp | 15 | ||||
| -rw-r--r-- | 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<Real, Vector<Scalar>> metropolis(const Tensor<Scalar, p>& J, const Ve return {E, z}; } + +template <class Real, class Scalar, int p, class Distribution, class Generator> +Vector<Scalar> randomSaddle(const Tensor<Scalar, p>& J, Distribution d, Generator& r, Real ε) { + Vector<Scalar> zSaddle; + bool foundSaddle = false; + + while (!foundSaddle) { + Vector<Scalar> z0 = normalize(randomVector<Scalar>(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 <getopt.h> +#include <limits> #include <stereographic.hpp> #include <unordered_map> #include <list> +#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<Real> d(0, 1, 0); - ComplexTensor J = generateCouplings<Complex, PSPIN_P>(N, complex_normal_distribution<Real>(0, σ, κ), r.engine()); + ComplexTensor J = generateCouplings<Complex, p>(N, complex_normal_distribution<Real>(0, σ, κ), r.engine()); std::function<Real(const ComplexTensor&, const ComplexVector&)> energyNormGrad = [] (const ComplexTensor& J, const ComplexVector& z) { @@ -149,106 +150,58 @@ int main(int argc, char* argv[]) { return W; }; - - std::function<Real(const ComplexTensor&, const ComplexVector&)> 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<Complex>(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<Real>::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<Complex>(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<ComplexMatrix> 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<Real>::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<Real>(J, zSaddle, zSaddleNext, 10, 4)) { + std::cerr << "Found a Stokes line" << std::endl; +// stokesLineTestNew<Real>(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<Complex>& g : test.gs) { + std::cout << g.transpose() << std::endl; } return 0; @@ -89,3 +89,18 @@ std::tuple<Real, Vector<Scalar>> WdW(const Tensor<Scalar, p>& J, const Vector<Sc return {W, dW}; } + +template <class Scalar> +Matrix<Scalar> dzDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) { + Real z² = z.squaredNorm(); + return (dH.conjugate() - (dH.dot(z) / z²) * z.conjugate()) * z.adjoint() / z²; +} + +template <class Scalar> +Matrix<Scalar> dzDotConjugate(const Vector<Scalar>& z, const Vector<Scalar>& dH, const Matrix<Scalar>& ddH) { + Real z² = z.squaredNorm(); + return -ddH + (ddH * z.conjugate()) * z.transpose() / z² + + (z.dot(dH) / z²) * ( + Matrix<Scalar>::Identity(ddH.rows(), ddH.cols()) - z.conjugate() * z.transpose() / z² + ); +} @@ -4,6 +4,8 @@ #include "complex_normal.hpp" #include "dynamics.hpp" +#include <iostream> + class ropeRelaxationStallException: public std::exception { virtual const char* what() const throw() { return "Gradient descent stalled."; @@ -11,49 +13,6 @@ class ropeRelaxationStallException: public std::exception { }; 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) { - Real z² = z.squaredNorm(); - Real z´² = z´.squaredNorm(); - - Vector<Scalar> ż = zDot(z, dH); - Real ż² = ż.squaredNorm(); - - Real Reż·z´ = real(ż.dot(z´)); - - Matrix<Scalar> dż = (dH.conjugate() - (dH.dot(z) / z²) * z.conjugate()) * z.adjoint() / z²; - Matrix<Scalar> dżc = -ddH + (ddH * z.conjugate()) * z.transpose() / z² - + (z.dot(dH) / z²) * ( - Matrix<Scalar>::Identity(ddH.rows(), ddH.cols()) - z.conjugate() * z.transpose() / z² - ); - - Vector<Scalar> dLdz = - ( - dżc * z´ + dż * z´.conjugate() - (dż * ż.conjugate() + dżc * ż) * Reż·z´ / ż² - ) / sqrt(ż² * z´²) / 2; - - Vector<Scalar> ż´ = -(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<Scalar> 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 Scalar> class Rope { public: std::vector<Vector<Scalar>> z; @@ -93,65 +52,10 @@ class Rope { return l; } - template <int p> - Real error(const Tensor<Scalar, p>& 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<Scalar> dz(unsigned i) const { return z[i + 1] - z[i - 1]; } - Vector<Scalar> ddz(unsigned i) const { - return 4.0 * (z[i + 1] + z[i - 1] - 2.0 * z[i]); - } - - template <int p> - std::vector<Vector<Scalar>> generateGradientδ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]); - - δ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 <int p> std::vector<Vector<Scalar>> generateDiscreteGradientδz(const Tensor<Scalar, p>& J, Real γ) const { std::vector<Vector<Scalar>> δ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<class Gen> - std::vector<Vector<Scalar>> generateRandomδz(Gen& r) const { - std::vector<Vector<Scalar>> δz(z.size()); + void spread() { + Real l = length(); + + Real a = 0; + unsigned pos = 0; + + std::vector<Vector<Scalar>> zNew = z; - complex_normal_distribution<> d(0, 1, 0); for (unsigned i = 1; i < z.size() - 1; i++) { - δz[i] = randomVector<Scalar>(z[0].size(), d, r); + Real b = i * l / (z.size() - 1); + + while (b > a) { + pos++; + a += (z[pos] - z[pos - 1]).norm(); + } + + Vector<Scalar> δz = z[pos] - z[pos - 1]; + + zNew[i] = normalize(z[pos] - (a - b) / δz.norm() * δz); } - return δz; + z = zNew; } template<int p> @@ -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<Vector<Scalar>> 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<Scalar> δz = z[pos] - z[pos - 1]; - - zNew[i] = normalize(z[pos] - (a - b) / δz.norm() * δz); - } - - z = zNew; - } - - template <int p> - void relaxGradient(const Tensor<Scalar, p>& J, unsigned N, Real δ0) { - Real δ = δ0; - try { - for (unsigned i = 0; i < N; i++) { - std::vector<Vector<Scalar>> δz = generateGradientδz(J); - δ = 1.1 * perturb(J, δ, δz); - } - } catch (std::exception& e) { - } - } - template <int p> void relaxDiscreteGradient(const Tensor<Scalar, p>& J, unsigned N, Real δ0, Real γ) { Real δ = δ0; try { for (unsigned i = 0; i < N; i++) { std::vector<Vector<Scalar>> δz = generateDiscreteGradientδz(J, γ); - δ = 1.1 * perturb(J, δ, δz, γ); + double stepSize = 0; + for (const Vector<Scalar>& 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 <int p, class Gen> - void relaxRandom(const Tensor<Scalar, p>& J, unsigned N, Real δ0, Gen& r) { - Real δ = δ0; - for (unsigned i = 0; i < N; i++) { - try { - std::vector<Vector<Scalar>> δz = generateRandomδz(r); - δ = 1.1 * perturb(J, δ, δz); - } catch (std::exception& e) { - } - } - } - template <int p> Real cost(const Tensor<Scalar, p>& J, Real γ = 0) const { Real c = 0; @@ -334,3 +207,180 @@ class Rope { return r; } }; + +template <class Real, class Scalar, int p> +bool stokesLineTest(const Tensor<Scalar, p>& J, const Vector<Scalar>& z1, const Vector<Scalar>& 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 Scalar> +class Cord { +public: + std::vector<Vector<Scalar>> gs; + Vector<Scalar> z0; + Vector<Scalar> z1; + + template <int p> + Cord(const Tensor<Scalar, p>& J, const Vector<Scalar>& z2, const Vector<Scalar>& z3, unsigned ng) : gs(ng, Vector<Scalar>::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<Scalar> f(Real t) const { + Vector<Scalar> z = (1 - t) * z0 + t * z1; + + for (unsigned i = 0; i < gs.size(); i++) { + z += gCoeff(i, t) * gs[i]; + } + + return z; + } + + Vector<Scalar> df(Real t) const { + Vector<Scalar> z = z1 - z0; + + for (unsigned i = 0; i < gs.size(); i++) { + z += dgCoeff(i, t) * gs[i]; + } + + return z; + } + + template <int p> + Real cost(const Tensor<Scalar, p>& J, Real t) const { + Vector<Scalar> z = f(t); + Scalar H; + Vector<Scalar> dH; + std::tie(H, dH, std::ignore) = hamGradHess(J, z); + Vector<Scalar> ż = zDot(z, dH); + Vector<Scalar> dz = df(t); + + return 1 - real(ż.dot(dz)) / ż.norm() / dz.norm(); + } + + template <int p> + Real totalCost(const Tensor<Scalar, p>& 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 <int p> + std::vector<Vector<Scalar>> dgs(const Tensor<Scalar, p>& J, Real t) const { + Vector<Scalar> z = f(t); + auto [H, dH, ddH] = hamGradHess(J, z); + Vector<Scalar> ż = zDot(z, dH); + Vector<Scalar> dz = df(t); + Matrix<Scalar> dż = dzDot(z, dH); + Matrix<Scalar> dżc = dzDotConjugate(z, dH, ddH); + + std::vector<Vector<Scalar>> x; + x.reserve(gs.size()); + + for (unsigned i = 0; i < gs.size(); i++) { + Real fdg = gCoeff(i, t); + Real dfdg = dgCoeff(i, t); + Vector<Scalar> 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 <int p> + std::pair<Real, Real> relaxStep(const Tensor<Scalar, p>& J, unsigned nt, Real δ₀) { + std::vector<Vector<Scalar>> dgsTot(gs.size(), Vector<Scalar>::Zero(z0.size())); + + for (unsigned i = 0; i < nt; i++) { + Real t = (i + 1.0) / (nt + 1.0); + std::vector<Vector<Scalar>> dgsI = dgs(J, t); + + for (unsigned j = 0; j < gs.size(); j++) { + dgsTot[j] += dgsI[j] / nt; + } + } + + Real stepSize = 0; + for (const Vector<Scalar>& dgi : dgsTot) { + stepSize += dgi.squaredNorm(); + } + stepSize = sqrt(stepSize); + + Cord cNew(*this); + + Real δ = δ₀; + Real oldCost = totalCost(J, nt); + Real newCost = std::numeric_limits<Real>::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 <int p> + void relax(const Tensor<Scalar, p>& J, unsigned nt, Real δ₀, unsigned maxSteps) { + Real δ = δ₀; + Real stepSize = std::numeric_limits<Real>::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++; + } + } +}; |