From 136fcddcd38d0b8f3b40faf7c1cb7365d9b2a753 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 15 Jan 2021 14:45:19 +0100 Subject: Converted more of library to templates accepting generic Scalar types and p --- dynamics.hpp | 117 ++++++++++++++++++++++++++++++++++++++++ langevin.cpp | 156 +++++++++--------------------------------------------- p-spin.hpp | 44 +++++++++------ stereographic.hpp | 26 +++++---- tensor.hpp | 20 ++++--- 5 files changed, 194 insertions(+), 169 deletions(-) create mode 100644 dynamics.hpp diff --git a/dynamics.hpp b/dynamics.hpp new file mode 100644 index 0000000..85eef71 --- /dev/null +++ b/dynamics.hpp @@ -0,0 +1,117 @@ +#pragma once + +#include +#include +#include + +#include "p-spin.hpp" +#include "stereographic.hpp" + +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) { + Vector z = z0; + double γ = γ0; + + auto [W, dW] = WdW(J, z); + + while (W > ε) { + Vector zNewTmp = z - γ * dW.conjugate(); + Vector zNew = normalize(zNewTmp); + + auto [WNew, dWNew] = WdW(J, zNew); + + if (WNew < W) { // If the step lowered the objective, accept it! + z = zNew; + W = WNew; + dW = dWNew; + γ = γ0; + } else { // Otherwise, shrink the step and try again. + γ /= δγ; + } + + if (γ < 1e-50) { + throw gradientDescentStall; + } + } + + return {W, z}; +} + +template +Vector findSaddle(const Tensor& J, const Vector& z0, double ε, double δW = 2, double γ0 = 1, double δγ = 2) { + Vector z = z0; + Vector ζ = euclideanToStereographic(z); + + double W; + std::tie(W, std::ignore) = WdW(J, z); + + Vector dH; + Matrix ddH; + std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ, z); + + while (W > ε) { + // ddH is complex symmetric, which is (almost always) invertible, so a + // partial pivot LU decomposition can be used. + Vector dζ = ddH.partialPivLu().solve(dH); + Vector ζNew = ζ - dζ; + Vector zNew = stereographicToEuclidean(ζNew); + + double WNew; + std::tie(WNew, std::ignore) = WdW(J, zNew); + + if (WNew < W) { // If Newton's step lowered the objective, accept it! + ζ = ζNew; + z = zNew; + W = WNew; + } else { // Otherwise, do gradient descent until W is a factor δW smaller. + std::tie(W, z) = gradientDescent(J, z, W / δW, γ0, δγ); + ζ = euclideanToStereographic(z); + } + + std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ, z); + } + + return z; +} + +template +Vector randomVector(unsigned N, Distribution d, Generator& r) { + Vector z(N); + + for (unsigned i = 0; i < N; i++) { + z(i) = d(r); + } + + return z; +} + +template +std::tuple> langevin(const Tensor& J, const Vector& z0, double T, double γ, unsigned N, Distribution d, Generator& r) { + Vector z = z0; + + double W; + std::tie(W, std::ignore) = WdW(J, z); + + std::uniform_real_distribution D(0, 1); + + for (unsigned i = 0; i < N; i++) { + Vector zNewTmp = z + randomVector(z.size(), d, r); + Vector zNew = normalize(zNewTmp); + + double WNew; + std::tie(WNew, std::ignore) = WdW(J, zNew); + + if (exp((W - WNew) / T) > D(r)) { + z = zNew; + W = WNew; + } + } + + return {W, z}; +} diff --git a/langevin.cpp b/langevin.cpp index a4c9cf3..870879b 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -1,127 +1,21 @@ -#include #include -#include -#include - #include "complex_normal.hpp" #include "p-spin.hpp" -#include "stereographic.hpp" +#include "dynamics.hpp" #include "pcg-cpp/include/pcg_random.hpp" #include "randutils/randutils.hpp" #include "unsupported/Eigen/CXX11/src/Tensor/TensorFFT.h" -using Rng = randutils::random_generator; - -Vector normalize(const Vector& z) { - return z * sqrt((double)z.size() / (Scalar)(z.transpose() * z)); -} - -template -Vector randomVector(unsigned N, Distribution d, Generator& r) { - Vector z(N); - - for (unsigned i = 0; i < N; i++) { - z(i) = d(r); - } - - return z; -} - -class gradientDescentStallException: public std::exception { - virtual const char* what() const throw() { - return "Gradient descent stalled."; - } -} gradientDescentStall; - -std::tuple gradientDescent(const Tensor& J, const Vector& z0, double ε, double γ0 = 1, double δγ = 2) { - Vector z = z0; - double γ = γ0; - - auto [W, dW] = WdW(J, z); - - while (W > ε) { - Vector zNew = normalize(z - γ * dW.conjugate()); - - auto [WNew, dWNew] = WdW(J, zNew); - - if (WNew < W) { // If the step lowered the objective, accept it! - z = zNew; - W = WNew; - dW = dWNew; - γ = γ0; - } else { // Otherwise, shrink the step and try again. - γ /= δγ; - } - - if (γ < 1e-50) { - throw gradientDescentStall; - } - } +#define PSPIN_P 3 +const int p = PSPIN_P; // polynomial degree of Hamiltonian +using Complex = std::complex; +using ComplexVector = Vector; +using ComplexMatrix = Matrix; +using ComplexTensor = Tensor; - return {W, z}; -} - -Vector findSaddle(const Tensor& J, const Vector& z0, double ε, double δW = 2, double γ0 = 1, double δγ = 2) { - Vector z = z0; - Vector ζ = euclideanToStereographic(z); - - double W; - std::tie(W, std::ignore) = WdW(J, z); - - Vector dH; - Matrix ddH; - std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ, z); - - while (W > ε) { - // ddH is complex symmetric, which is (almost always) invertible, so a - // partial pivot LU decomposition can be used. - Vector dζ = ddH.partialPivLu().solve(dH); - Vector ζNew = ζ - dζ; - Vector zNew = stereographicToEuclidean(ζNew); - - double WNew; - std::tie(WNew, std::ignore) = WdW(J, zNew); - - if (WNew < W) { // If Newton's step lowered the objective, accept it! - ζ = ζNew; - z = zNew; - W = WNew; - } else { // Otherwise, do gradient descent until W is a factor δW smaller. - std::tie(W, z) = gradientDescent(J, z, W / δW, γ0, δγ); - ζ = euclideanToStereographic(z); - } - - std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ, z); - } - - return z; -} - -std::tuple langevin(const Tensor& J, const Vector& z0, double T, double γ, unsigned N, Rng& r) { - Vector z = z0; - - double W; - std::tie(W, std::ignore) = WdW(J, z); - - complex_normal_distribution<> d(0, γ, 0); - - for (unsigned i = 0; i < N; i++) { - Vector dz = randomVector(z.size(), d, r.engine()); - Vector zNew = normalize(z + dz); - - double WNew; - std::tie(WNew, std::ignore) = WdW(J, zNew); - - if (exp((W - WNew) / T) > r.uniform(0.0, 1.0)) { - z = zNew; - W = WNew; - } - } - - return {W, z}; -} +using Rng = randutils::random_generator; int main(int argc, char* argv[]) { // model parameters @@ -178,22 +72,24 @@ int main(int argc, char* argv[]) { } } - Scalar κ(Rκ, Iκ); + Complex κ(Rκ, Iκ); double σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1))); Rng r; - Tensor J = generateCouplings(N, complex_normal_distribution<>(0, σ, κ), r.engine()); - Vector z0 = normalize(randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine())); + 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())); - Vector zSaddle = findSaddle(J, z0, ε); - Vector zSaddlePrev = Vector::Zero(N); - Vector z = zSaddle; + ComplexVector zSaddle = findSaddle(J, z0, ε); + ComplexVector zSaddlePrev = ComplexVector::Zero(N); + ComplexVector z = zSaddle; while (δ < (zSaddle - zSaddlePrev).norm()) { // Until we find two saddles sufficiently close... - std::tie(std::ignore, z) = langevin(J, z, T, γ, M, r); + std::tie(std::ignore, z) = langevin(J, z, T, γ, M, d, r.engine()); try { - Vector zSaddleNext = findSaddle(J, z, ε); + ComplexVector zSaddleNext = findSaddle(J, z, ε); if (Δ < (zSaddleNext - zSaddle).norm()) { // Ensure we are finding distinct saddles. zSaddlePrev = zSaddle; zSaddle = zSaddleNext; @@ -209,20 +105,20 @@ int main(int argc, char* argv[]) { complex_normal_distribution<> dJ(0, εJ * σ, 0); - std::function)> perturbJ = - [&dJ, &r] (Tensor& JJ, std::array ind) { - Scalar Ji = getJ(JJ, ind); - setJ(JJ, ind, Ji + dJ(r.engine())); + std::function)> perturbJ = + [&dJ, &r] (ComplexTensor& JJ, std::array ind) { + Complex Ji = getJ(JJ, ind); + setJ(JJ, ind, Ji + dJ(r.engine())); }; for (unsigned i = 0; i < n; i++) { - Tensor Jp = J; + ComplexTensor Jp = J; - iterateOver(Jp, perturbJ); + iterateOver(Jp, perturbJ); try { - Vector zSaddleNew = findSaddle(Jp, zSaddle, ε); - Vector zSaddlePrevNew = findSaddle(Jp, zSaddlePrev, ε); + ComplexVector zSaddleNew = findSaddle(Jp, zSaddle, ε); + ComplexVector zSaddlePrevNew = findSaddle(Jp, zSaddlePrev, ε); std::cout << zSaddleNew.transpose() << " " << zSaddlePrevNew.transpose() << std::endl; } catch (std::exception& e) { diff --git a/p-spin.hpp b/p-spin.hpp index 3ef10e7..4621db6 100644 --- a/p-spin.hpp +++ b/p-spin.hpp @@ -3,49 +3,59 @@ #include #include "tensor.hpp" +#include "factorial.hpp" -#define PSPIN_P 3 -const unsigned p = PSPIN_P; // polynomial degree of Hamiltonian +template +using Vector = Eigen::Matrix; -using Scalar = std::complex; -using Vector = Eigen::VectorXcd; -using Matrix = Eigen::MatrixXcd; -using Tensor = Eigen::Tensor; +template +using Matrix = Eigen::Matrix; -std::tuple hamGradHess(const Tensor& J, const Vector& z) { - Matrix Jz = contractDown(J, z); // Contracts J into p - 2 copies of z. - Vector Jzz = Jz * z; +template +using Tensor = Eigen::Tensor; + +template +std::tuple, Matrix> hamGradHess(const Tensor& J, const Vector& z) { + Matrix Jz = contractDown(J, z); // Contracts J into p - 2 copies of z. + Vector Jzz = Jz * z; Scalar Jzzz = Jzz.transpose() * z; double pBang = factorial(p); - Matrix hessian = ((p - 1) * p / pBang) * Jz; - Vector gradient = (p / pBang) * Jzz; + Matrix hessian = ((p - 1) * p / pBang) * Jz; + Vector gradient = (p / pBang) * Jzz; Scalar hamiltonian = Jzzz / pBang; return {hamiltonian, gradient, hessian}; } -Vector project(const Vector& z, const Vector& x) { +template +Vector project(const Vector& z, const Vector& x) { Scalar xz = x.transpose() * z; return x - (xz / z.squaredNorm()) * z.conjugate(); } -std::tuple WdW(const Tensor& J, const Vector& z) { - Vector dH; - Matrix ddH; +template +std::tuple> WdW(const Tensor& J, const Vector& z) { + Vector dH; + Matrix ddH; std::tie(std::ignore, dH, ddH) = hamGradHess(J, z); - Vector dzdt = project(z, dH.conjugate()); + Vector dzdt = project(z, dH.conjugate().eval()); double a = z.squaredNorm(); Scalar A = (Scalar)(z.transpose() * dzdt) / a; Scalar B = dH.dot(z) / a; double W = dzdt.squaredNorm(); - Vector dW = ddH * (dzdt - A * z.conjugate()) + Vector dW = ddH * (dzdt - A * z.conjugate()) + 2 * (conj(A) * B * z).real() - conj(B) * dzdt - conj(A) * dH.conjugate(); return {W, dW}; } + +template +Vector normalize(const Vector& z) { + return z * sqrt((double)z.size() / (Scalar)(z.transpose() * z)); +} diff --git a/stereographic.hpp b/stereographic.hpp index 638923f..2dfa735 100644 --- a/stereographic.hpp +++ b/stereographic.hpp @@ -4,9 +4,10 @@ #include "p-spin.hpp" -Vector stereographicToEuclidean(const Vector& ζ) { +template +Vector stereographicToEuclidean(const Vector& ζ) { unsigned N = ζ.size() + 1; - Vector z(N); + Vector z(N); Scalar a = ζ.transpose() * ζ; Scalar b = 2 * sqrt(N) / (1.0 + a); @@ -20,9 +21,10 @@ Vector stereographicToEuclidean(const Vector& ζ) { return b * z; } -Vector euclideanToStereographic(const Vector& z) { +template +Vector euclideanToStereographic(const Vector& z) { unsigned N = z.size(); - Vector ζ(N - 1); + Vector ζ(N - 1); for (unsigned i = 0; i < N - 1; i++) { ζ(i) = z(i); @@ -31,9 +33,10 @@ Vector euclideanToStereographic(const Vector& z) { return ζ / (sqrt(N) - z(N - 1)); } -Matrix stereographicJacobian(const Vector& ζ) { +template +Matrix stereographicJacobian(const Vector& ζ) { unsigned N = ζ.size(); - Matrix J(N, N + 1); + Matrix J(N, N + 1); Scalar b = 1.0 + (Scalar)(ζ.transpose() * ζ); @@ -52,15 +55,16 @@ Matrix stereographicJacobian(const Vector& ζ) { return 4 * sqrt(N + 1) * J / pow(b, 2); } -std::tuple stereographicHamGradHess(const Tensor& J, const Vector& ζ, const Vector& z) { +template +std::tuple, Matrix> stereographicHamGradHess(const Tensor& J, const Vector& ζ, const Vector& z) { auto [hamiltonian, gradZ, hessZ] = hamGradHess(J, z); - Matrix jacobian = stereographicJacobian(ζ); + Matrix jacobian = stereographicJacobian(ζ); - Matrix metric = jacobian * jacobian.adjoint(); + Matrix metric = jacobian * jacobian.adjoint(); // The metric is Hermitian and positive definite, so a Cholesky decomposition can be used. - Vector grad = metric.llt().solve(jacobian) * gradZ; - Matrix hess = jacobian * hessZ * jacobian.transpose(); + Vector grad = metric.llt().solve(jacobian) * gradZ; + Matrix hess = jacobian * hessZ * jacobian.transpose(); return {hamiltonian, grad, hess}; } diff --git a/tensor.hpp b/tensor.hpp index 66f6d59..21f2a89 100644 --- a/tensor.hpp +++ b/tensor.hpp @@ -5,43 +5,41 @@ #include -#include "factorial.hpp" - -template +template Eigen::Tensor initializeJHelper(unsigned N, std::index_sequence) { std::array Ns; std::fill_n(Ns.begin(), p, N); return Eigen::Tensor(std::get(Ns)...); } -template +template Eigen::Tensor initializeJ(unsigned N) { return initializeJHelper(N, std::make_index_sequence

()); } -template +template void setJHelper(Eigen::Tensor& J, const std::array& ind, Scalar val, std::index_sequence) { J(std::get(ind)...) = val; } -template +template void setJ(Eigen::Tensor& J, std::array ind, Scalar val) { do { setJHelper(J, ind, val, std::make_index_sequence

()); } while (std::next_permutation(ind.begin(), ind.end())); } -template +template Scalar getJHelper(const Eigen::Tensor& J, const std::array& ind, std::index_sequence) { return J(std::get(ind)...); } -template +template Scalar getJ(const Eigen::Tensor& J, const std::array& ind) { return getJHelper(J, ind, std::make_index_sequence

()); } -template +template void iterateOverHelper(Eigen::Tensor& J, std::function&, std::array)>& f, unsigned l, std::array is) { @@ -62,13 +60,13 @@ void iterateOverHelper(Eigen::Tensor& J, } } -template +template void iterateOver(Eigen::Tensor& J, std::function&, std::array)>& f) { std::array is; iterateOverHelper(J, f, p, is); } -template +template Eigen::Tensor generateCouplings(unsigned N, Distribution d, Generator& r) { Eigen::Tensor J = initializeJ(N); -- cgit v1.2.3-70-g09d2