diff options
author | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2021-01-15 14:45:19 +0100 |
---|---|---|
committer | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2021-01-15 14:45:19 +0100 |
commit | 136fcddcd38d0b8f3b40faf7c1cb7365d9b2a753 (patch) | |
tree | 38723818277fda2ed2573ee4479cc2b9a66c39a9 /dynamics.hpp | |
parent | c10b0a99ecb9a6aab144aeca9c7538d1a897fd49 (diff) | |
download | code-136fcddcd38d0b8f3b40faf7c1cb7365d9b2a753.tar.gz code-136fcddcd38d0b8f3b40faf7c1cb7365d9b2a753.tar.bz2 code-136fcddcd38d0b8f3b40faf7c1cb7365d9b2a753.zip |
Converted more of library to templates accepting generic Scalar types and p
Diffstat (limited to 'dynamics.hpp')
-rw-r--r-- | dynamics.hpp | 117 |
1 files changed, 117 insertions, 0 deletions
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 <exception> +#include <eigen3/Eigen/LU> +#include <random> + +#include "p-spin.hpp" +#include "stereographic.hpp" + +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) { + Vector<Scalar> z = z0; + double γ = γ0; + + auto [W, dW] = WdW(J, z); + + while (W > ε) { + Vector<Scalar> zNewTmp = z - γ * dW.conjugate(); + Vector<Scalar> 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 <class Scalar, int p> +Vector<Scalar> findSaddle(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, double ε, double δW = 2, double γ0 = 1, double δγ = 2) { + Vector<Scalar> z = z0; + Vector<Scalar> ζ = euclideanToStereographic(z); + + double W; + std::tie(W, std::ignore) = WdW(J, z); + + Vector<Scalar> dH; + Matrix<Scalar> 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<Scalar> dζ = ddH.partialPivLu().solve(dH); + Vector<Scalar> ζNew = ζ - dζ; + Vector<Scalar> 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 <class Scalar, class Distribution, class Generator> +Vector<Scalar> randomVector(unsigned N, Distribution d, Generator& r) { + Vector<Scalar> z(N); + + for (unsigned i = 0; i < N; i++) { + z(i) = d(r); + } + + return z; +} + +template <class Scalar, int p, class Distribution, class Generator> +std::tuple<double, Vector<Scalar>> langevin(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, double T, double γ, unsigned N, Distribution d, Generator& r) { + Vector<Scalar> z = z0; + + double W; + std::tie(W, std::ignore) = WdW(J, z); + + std::uniform_real_distribution<double> D(0, 1); + + for (unsigned i = 0; i < N; i++) { + Vector<Scalar> zNewTmp = z + randomVector<Scalar>(z.size(), d, r); + Vector<Scalar> 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}; +} |