diff options
| -rw-r--r-- | complex_normal.hpp | 26 | ||||
| -rw-r--r-- | factorial.hpp | 9 | ||||
| -rw-r--r-- | langevin.cpp | 318 | ||||
| -rw-r--r-- | stereographic.hpp | 88 | ||||
| -rw-r--r-- | tensor.hpp | 65 |
5 files changed, 242 insertions, 264 deletions
diff --git a/complex_normal.hpp b/complex_normal.hpp new file mode 100644 index 0000000..0819758 --- /dev/null +++ b/complex_normal.hpp @@ -0,0 +1,26 @@ +#pragma once +#include <complex> +#include <random> + +template <class T = double> class complex_normal_distribution { +private: + std::normal_distribution<T> d; + double δx; + double δy; + std::complex<T> μ; + std::complex<T> eθ; + +public: + complex_normal_distribution(std::complex<T> μ, T σ, std::complex<T> κ) : μ(μ), d(0, σ / sqrt(2)) { + δx = sqrt(1 + std::abs(κ)); + δy = sqrt(1 - std::abs(κ)); + eθ = std::polar(1.0, std::arg(κ)); + } + + template <class Generator> std::complex<T> operator()(Generator& g) { + // First generate an axis-aligned complex number centered at the origin. + std::complex<T> z(d(g) * δx, d(g) * δy); + // Shift and rotate the number. + return μ + eθ * z; + } +}; diff --git a/factorial.hpp b/factorial.hpp new file mode 100644 index 0000000..ddc0bac --- /dev/null +++ b/factorial.hpp @@ -0,0 +1,9 @@ +#pragma once + +long unsigned factorial(unsigned p) { + if (p == 0) { + return 1; + } else { + return p * factorial(p - 1); + } +} diff --git a/langevin.cpp b/langevin.cpp index b6f16d7..7590d2b 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -14,94 +14,19 @@ #include "pcg-cpp/include/pcg_random.hpp" #include "randutils/randutils.hpp" +#include "complex_normal.hpp" +#include "tensor.hpp" + #define P_SPIN_P 3 const unsigned p = P_SPIN_P; using Scalar = std::complex<double>; -using TensorP = Eigen::Tensor<Scalar, p>; -using Tensor0 = Eigen::Tensor<Scalar, 0>; -using Tensor1 = Eigen::Tensor<Scalar, 1>; -using Tensor2 = Eigen::Tensor<Scalar, 2>; -using Tensor3 = Eigen::Tensor<Scalar, 3>; -using Matrix = Eigen::MatrixXcd; using Vector = Eigen::VectorXcd; - -const std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(0,0)}; +using Matrix = Eigen::MatrixXcd; +using Tensor = Eigen::Tensor<Scalar, p>; using Rng = randutils::random_generator<pcg32>; -template <class T = double> -class complex_normal_distribution { - private: - std::complex<T> μ; - T σ; - std::complex<T> κ; - std::normal_distribution<T> dist; - - public: - complex_normal_distribution(std::complex<T> μ, T σ, std::complex<T> κ) : - μ(μ), σ(σ), κ(κ), dist(0, σ / sqrt(2)) {} - - template <class Generator> - std::complex<T> operator()(Generator& g) { - std::complex<T> x(dist(g) * sqrt(1 + std::abs(κ)), dist(g) * sqrt(1 - std::abs(κ))); - return μ + std::polar((T)1, std::arg(κ)) * x; - } -}; - -long unsigned factorial(unsigned p) { - if (p == 0) { - return 1; - } else { - return p * factorial(p - 1); - } -} - -/* -void populateCouplings(TensorP& J, unsigned level, std::list<unsigned> indices, complex_normal_distribution<> dist, Rng& r) { - if (level == 0) { - Scalar x = dist(r.engine()); - do { - } while (std::next_permutation(indices.begin(), indices.end())) - } -} -*/ - -Tensor3 generateCouplings(unsigned N, std::complex<double> κ, Rng& r) { - double σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1))); - complex_normal_distribution<> dist(0, σ, κ); - - TensorP J(N, N, N); - - for (unsigned i1 = 0; i1 < N; i1++) { - for (unsigned i2 = i1; i2 < N; i2++) { - for (unsigned i3 = i2; i3 < N; i3++) { - std::complex<double> x = dist(r.engine()); - - J(i1, i2, i3) = x; - J(i1, i3, i2) = x; - J(i2, i1, i3) = x; - J(i2, i3, i1) = x; - J(i3, i1, i2) = x; - J(i3, i2, i1) = x; - } - } - } - - return J; -} - -template <int r> -Eigen::Tensor<std::complex<double>, r> conj(const Eigen::Tensor<std::complex<double>, r>& t) { - return t.unaryExpr(static_cast<std::complex<double> (*)(const std::complex<double>&)>(&std::conj)); -} - -template <int r> -double norm(const Eigen::Tensor<std::complex<double>, r>& t) { - Eigen::Tensor<double, 0> t2 = t.unaryExpr(static_cast<double (*)(const std::complex<double>&)>(&std::norm)).sum(); - return t2(0); -} - Vector initializeVector(unsigned N, Rng& r) { Vector z(N); complex_normal_distribution<> dist(0, 1, 0); @@ -115,207 +40,64 @@ Vector initializeVector(unsigned N, Rng& r) { return z; } -template <int r> -Eigen::Tensor<Scalar, 3> contractDown(const Eigen::Tensor<Scalar, r>& J, const Eigen::Tensor<Scalar, 1>& z) { - return contractDown<r - 1>(J.contract(z, ip), z); -} +std::tuple<std::complex<double>, Vector, Matrix> hamGradHess(const Tensor& J, const Vector& z) { + Matrix Jz = contractDown(J, z); + Vector Jzz = Jz * z; -template <> -Eigen::Tensor<Scalar, 3> contractDown(const Eigen::Tensor<Scalar, 3>& J, const Eigen::Tensor<Scalar, 1>& z) { - return J; -} - -std::tuple<std::complex<double>, Vector, Matrix> hamGradHess(const Tensor3& J, const Vector& z) { - Tensor1 zT = Eigen::TensorMap<Eigen::Tensor<const std::complex<double>, 1>>(z.data(), {z.size()}); - - Tensor2 J1 = zT.contract(J, ip); - Tensor1 J2 = zT.contract(J1, ip); - Tensor0 J3 = zT.contract(J2, ip); + Matrix hessian = ((p - 1) * (double)p / factorial(p)) * Jz; + Vector gradient = ((double)p / factorial(p)) * Jzz; + Scalar hamiltonian = (1.0 / factorial(p)) * Jzz.dot(z); - Matrix hess = ((p - 1) * (double)p / factorial(p)) * Eigen::Map<Matrix>(J1.data(), z.size(), z.size()); - Vector grad = ((double)p / factorial(p)) * Eigen::Map<Vector>(J2.data(), z.size()); - std::complex<double> hamiltonian = (1.0 / factorial(p)) * J3(0); - - return {hamiltonian, grad, hess}; + return {hamiltonian, gradient, hessian}; } -Vector stereographicToEuclidean(const Vector& ζ) { - unsigned N = ζ.size() + 1; - Vector z(N); - Scalar a = ζ.dot(ζ); - std::complex<double> b = 2.0 * sqrt(N) / (1.0 + a); - - for (unsigned i = 0; i < N - 1; i++) { - z(i) = b * ζ(i); - } - - z(N - 1) = b * (a - 1.0) / 2.0; +std::tuple<double, Vector> WdW(const Tensor& J, const Vector& z) { + Vector gradient; + Matrix hessian; + std::tie(std::ignore, gradient, hessian) = hamGradHess(J, z); - return z; -} + Vector projectedGradient = gradient - (gradient.dot(z) / (double)z.size()) * z; -Vector euclideanToStereographic(const Vector& z) { - unsigned N = z.size(); - Vector ζ(N - 1); - - for (unsigned i = 0; i < N - 1; i++) { - ζ(i) = z(i) / (sqrt(N) - z(N - 1)); - } + double W = projectedGradient.cwiseAbs2().sum(); + Vector dW = hessian.conjugate() * projectedGradient; - return ζ; + return {W, dW}; } -Matrix stereographicJacobian(const Vector& ζ) { - unsigned N = ζ.size() + 1; - Matrix J(N - 1, N); - - Scalar b = ζ.dot(ζ); - - for (unsigned i = 0; i < N - 1; i++) { - for (unsigned j = 0; j < N - 1; j++) { - J(i, j) = - 4.0 * ζ(i) * ζ(j); - if (i == j) { - J(i, j) += 2.0 * (1.0 + b); - } - J(i, j) *= sqrt(N) / pow(1.0 + b, 2); - } - - J(i, N - 1) = 4.0 * sqrt(N) * ζ(i) / pow(1.0 + b, 2); - } - - return J; -} +Vector langevin(const Tensor& J, const Vector& z0, double T, double γ0, std::function<bool(double, unsigned)> quit, Rng& r) { + Vector z = z0; -Matrix stereographicMetric(const Vector& ζ) { - unsigned N = ζ.size(); - Matrix g(N, N); + double W; + Vector dW; + std::tie(W, dW) = WdW(J, z); - double a = ζ.cwiseAbs2().sum(); - Scalar b = ζ.dot(ζ); + unsigned steps = 0; + complex_normal_distribution<> d(0, T, 0); - for (unsigned i = 0; i < N; i++) { - for (unsigned j = 0; j < N; j++) { - g(i, j) = 16.0 * (std::conj(ζ(i)) * ζ(j) * (1.0 + a) - std::real(std::conj(ζ(i) * ζ(j)) * (1.0 + b))); - if (i == j) { - g(i, j) += 4.0 * std::abs(1.0 + b); - } - g(i, j) *= (N + 1) / std::norm(pow(b, 2)); + while (!quit(W, steps)) { + double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2); + Vector η(z.size()); + for (unsigned i = 0; i < z.size(); i++) { + η(i) = d(r.engine()); } - } - - return g; -} - -std::tuple<std::complex<double>, Vector, Matrix> stereographicHamGradHess(const Tensor3& J, const Vector& ζ) { - Vector grad; - Matrix hess; - - Matrix jacobian = stereographicJacobian(ζ); - Vector z = stereographicToEuclidean(ζ); - - std::complex<double> hamiltonian; - Vector gradZ; - Matrix hessZ; - std::tie(hamiltonian, gradZ, hessZ) = hamGradHess(J, z); - - Matrix g = stereographicMetric(ζ); - Matrix gj = g.llt().solve(jacobian); - - grad = gj * gradZ; - hess = gj * hessZ * gj.transpose(); - return {hamiltonian, grad, hess}; -} - -std::tuple<double, Vector, Matrix> WdW(const Tensor3& J, const Vector& z) { - Vector grad; - Matrix hess; - - std::tie(std::ignore, grad, hess) = hamGradHess(J, z); - - Tensor1 gradT = Eigen::TensorMap<Eigen::Tensor<const Scalar, 1>>(grad.data(), {z.size()}); - - Vector gtmp = grad - (grad.dot(z) / (double)z.size()) * z; - double W = gtmp.cwiseAbs2().sum(); - Vector dW = (hess - (z.dot(grad) * Matrix::Identity(z.size(), z.size()) + grad * z.transpose() + z * (hess * z).transpose()) / (double)z.size()).conjugate() * gtmp; - Tensor2 ddWT = ((p - 2) * (p - 1) * (double)p / factorial(p)) * conj(J).contract(gradT, ip); - Matrix ddW = Eigen::Map<Matrix>(ddWT.data(), z.size(), z.size()); + Vector zNext = z - γ * dW + η; + zNext *= sqrt(zNext.size()) / sqrt(zNext.dot(zNext)); - return {W, dW, ddW}; -} + double WNext; + Vector dWNext; + std::tie(WNext, dWNext) = WdW(J, zNext); -Vector findSaddle(const Tensor3& J, const Vector& z0, double δ, double γ) { - Vector z = z0; - double W; - Vector dW; - Matrix ddW; - - std::tie(W, dW, ddW) = WdW(J, z); - - while (W > δ * z.size()) { - // Vector dz = ddW.ldlt().solve(dW); - Vector dz = dW; - while (true) { - Vector newz = z - γ * dz; - newz *= sqrt(newz.size()) / sqrt(newz.dot(newz)); - - double newW; - Vector newdW; - Matrix newddW; - std::tie(newW, newdW, newddW) = WdW(J, newz); - - if (newW < W) { - γ *= 1.01; - z = newz; - W = newW; - dW= newdW; - ddW = newddW; - break; - } else { - γ /= 0.9; - } + if (exp((W - WNext) / T) > r.uniform(0.0, 1.0)) { + z = zNext; + W = WNext; + dW = dWNext; } - std::cout << W << std::endl; } return z; } -/* -Vector generateKick(const Vector& z, Rng& r) { - Vector k(z.size()); - - for (unsigned i = 0; i < k.size(); i++) { - k(i) = std::complex<double>(r.variate<double>(0, 1), r.variate<double>(0, 1)) / sqrt(2); - } - - Scalar kz = z.contract(k, ip); - k = k - kz(0) * z; - - return k; -} - -void langevin(const Tensor& J, Vector& z, std::complex<double>& ε, double T, unsigned t, double γ, Rng& r) { - double W; - Vector dz; - - std::tie(W, dz) = WdW(J, z); - - for (unsigned i = 0; i < t; i++) { - std::cout << W << std::endl; - - Vector η = generateKick(z, r); - Vector δz = (γ / z.size()) * (-dz + T * η); - Scalar δε = dεdz.contract(δz, ip); - - z = z + δz; - ε = ε + δε(0); - - std::tie(W, dz, std::ignore, dεdz) = WdW(J, z, ε); - } -} -*/ - int main(int argc, char* argv[]) { // model parameters unsigned N = 10; // number of spins @@ -358,20 +140,28 @@ int main(int argc, char* argv[]) { } std::complex<double> κ(Rκ, Iκ); + double σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1))); Rng r; + complex_normal_distribution<> d(0, σ, κ); - Tensor3 J = generateCouplings(N, κ, r); + Tensor J = generateCouplings<std::complex<double>, p>(N, d, r.engine()); Vector z = initializeVector(N, r); - Vector zm = findSaddle(J, z, δ, γ); + std::function<bool(double, unsigned)> findSaddle = [δ](double W, unsigned) { + std::cout << W << std::endl; + return W < δ; + }; + + Vector zm = langevin(J, z, T, γ, findSaddle, r); std::complex<double> H; Vector dH; std::tie(H, dH, std::ignore) = hamGradHess(J, zm); - std::cout << zm << std::endl; - std::cout << dH - (H / (double)N) * zm << std::endl; -} + Vector constraint = dH - ((double)p * H / (double)N) * zm; + std::cout << std::endl << zm.dot(zm) << std::endl; + std::cout << constraint.cwiseAbs2().sum() << std::endl; +} diff --git a/stereographic.hpp b/stereographic.hpp new file mode 100644 index 0000000..7a6b411 --- /dev/null +++ b/stereographic.hpp @@ -0,0 +1,88 @@ + +Vector stereographicToEuclidean(const Vector& ζ) { + unsigned N = ζ.size() + 1; + Vector z(N); + Scalar a = ζ.dot(ζ); + std::complex<double> b = 2.0 * sqrt(N) / (1.0 + a); + + for (unsigned i = 0; i < N - 1; i++) { + z(i) = b * ζ(i); + } + + z(N - 1) = b * (a - 1.0) / 2.0; + + return z; +} + +Vector euclideanToStereographic(const Vector& z) { + unsigned N = z.size(); + Vector ζ(N - 1); + + for (unsigned i = 0; i < N - 1; i++) { + ζ(i) = z(i) / (sqrt(N) - z(N - 1)); + } + + return ζ; +} + +Matrix stereographicJacobian(const Vector& ζ) { + unsigned N = ζ.size() + 1; + Matrix J(N - 1, N); + + Scalar b = ζ.dot(ζ); + + for (unsigned i = 0; i < N - 1; i++) { + for (unsigned j = 0; j < N - 1; j++) { + J(i, j) = - 4.0 * ζ(i) * ζ(j); + if (i == j) { + J(i, j) += 2.0 * (1.0 + b); + } + J(i, j) *= sqrt(N) / pow(1.0 + b, 2); + } + + J(i, N - 1) = 4.0 * sqrt(N) * ζ(i) / pow(1.0 + b, 2); + } + + return J; +} + +Matrix stereographicMetric(const Vector& ζ) { + unsigned N = ζ.size(); + Matrix g(N, N); + + double a = ζ.cwiseAbs2().sum(); + Scalar b = ζ.dot(ζ); + + for (unsigned i = 0; i < N; i++) { + for (unsigned j = 0; j < N; j++) { + g(i, j) = 16.0 * (std::conj(ζ(i)) * ζ(j) * (1.0 + a) - std::real(std::conj(ζ(i) * ζ(j)) * (1.0 + b))); + if (i == j) { + g(i, j) += 4.0 * std::abs(1.0 + b); + } + g(i, j) *= (N + 1) / std::norm(pow(b, 2)); + } + } + + return g; +} + +std::tuple<std::complex<double>, Vector, Matrix> stereographicHamGradHess(const Tensor3& J, const Vector& ζ) { + Vector grad; + Matrix hess; + + Matrix jacobian = stereographicJacobian(ζ); + Vector z = stereographicToEuclidean(ζ); + + std::complex<double> hamiltonian; + Vector gradZ; + Matrix hessZ; + std::tie(hamiltonian, gradZ, hessZ) = hamGradHess(J, z); + + Matrix g = stereographicMetric(ζ); + Matrix gj = g.llt().solve(jacobian); + + grad = gj * gradZ; + hess = gj * hessZ * gj.transpose(); + + return {hamiltonian, grad, hess}; +} diff --git a/tensor.hpp b/tensor.hpp new file mode 100644 index 0000000..8345d83 --- /dev/null +++ b/tensor.hpp @@ -0,0 +1,65 @@ +#pragma once + +#include <algorithm> +#include <array> +#include <utility> +#include <eigen3/unsupported/Eigen/CXX11/Tensor> + +#include "factorial.hpp" + +template <class Scalar, unsigned p, std::size_t... Indices> +void setJ(Scalar z, Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& is, std::index_sequence<Indices...>) { + J(std::get<Indices>(is)...) = z; +} + +template <class Scalar, unsigned p, std::size_t... Indices> +Eigen::Tensor<Scalar, p> initializeJ(unsigned N, std::index_sequence<Indices...>) { + std::array<unsigned, p> Ns; + std::fill_n(Ns.begin(), p, N); + return Eigen::Tensor<Scalar, p>(std::get<Indices>(Ns)...); +} + +template <class Scalar, unsigned p, class Distribution, class Generator> +void populateCouplings(Eigen::Tensor<Scalar, p>& J, unsigned N, unsigned l, std::array<unsigned, p> is, Distribution d, Generator& r) { + if (l == 0) { + Scalar z = d(r); + do { + setJ<Scalar, p>(z, J, is, std::make_index_sequence<p>()); + } while (std::next_permutation(is.begin(), is.end())); + } else { + unsigned iMin; + if (l == p) { + iMin = 0; + } else { + iMin = is[p - l - 1]; + } + for (unsigned i = iMin; i < N; i++) { + std::array<unsigned, p> js = is; + js[p - l] = i; + populateCouplings<Scalar, p, Distribution, Generator>(J, N, l - 1, js, d, r); + } + } +} + +template <class Scalar, unsigned p, class Distribution, class Generator> +Eigen::Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator& r) { + Eigen::Tensor<Scalar, p> J = initializeJ<Scalar, p>(N, std::make_index_sequence<p>()); + + std::array<unsigned, p> is; + populateCouplings<Scalar, p, Distribution, Generator>(J, N, p, is, d, r); + + return J; +} + +template <class Scalar> +Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> contractDown(const Eigen::Tensor<Scalar, 2>& J, const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>& z) { + return Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>(J.data(), z.size(), z.size()); +} + +template <class Scalar, int r> +Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> contractDown(const Eigen::Tensor<Scalar, r>& J, const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>& z) { + Eigen::Tensor<Scalar, 1> zT = Eigen::TensorMap<Eigen::Tensor<const Scalar, 1>>(z.data(), {z.size()}); + std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(0,0)}; + Eigen::Tensor<Scalar, r - 1> Jz = J.contract(zT, ip); + return contractDown(Jz, z); +} |