From 1c79cff86e5cfae48e00184842101a9678b89903 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 5 Jan 2021 10:59:51 +0100
Subject: Lots of work, refactoring.
---
complex_normal.hpp | 26 +++++
factorial.hpp | 9 ++
langevin.cpp | 318 +++++++++--------------------------------------------
stereographic.hpp | 88 +++++++++++++++
tensor.hpp | 65 +++++++++++
5 files changed, 242 insertions(+), 264 deletions(-)
create mode 100644 complex_normal.hpp
create mode 100644 factorial.hpp
create mode 100644 stereographic.hpp
create mode 100644 tensor.hpp
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);
+}
--
cgit v1.2.3-70-g09d2