From 92c319247fae9d8be3f3291d6fa64266224ef61c Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 2 Apr 2024 17:24:14 +0200 Subject: Added tensor support. --- factorial.hpp | 9 +++ least_squares.cpp | 174 ++++++++++++++++++++++++++++++------------------------ tensor.hpp | 108 +++++++++++++++++++++++++++++++++ types.hpp | 3 + 4 files changed, 216 insertions(+), 78 deletions(-) create mode 100644 factorial.hpp create mode 100644 tensor.hpp create mode 100644 types.hpp 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/least_squares.cpp b/least_squares.cpp index 093892d..28397ec 100644 --- a/least_squares.cpp +++ b/least_squares.cpp @@ -1,122 +1,136 @@ #include -#include #include #include "pcg-cpp/include/pcg_random.hpp" #include "randutils/randutils.hpp" +#include "tensor.hpp" -template -using Vector = Eigen::Matrix; - -template -using Matrix = Eigen::Matrix; +Vector normalize(const Vector& z) { + return z * sqrt((Real)z.size() / (Real)(z.transpose() * z)); +} -template +template class Model { private: - Matrix A; - Vector b; + std::tuple...> Js; public: - template - Model(Real σ, unsigned N, unsigned M, Generator& r) : A(M, N), b(M) { - std::normal_distribution aDistribution(0, 1 / sqrt(N)); - - for (unsigned i = 0; i < M; i++) { - for (unsigned j =0; j < N; j++) { - A(i, j) = aDistribution(r); - } - } - - std::normal_distribution bDistribution(0, σ); - - for (unsigned i = 0; i < M; i++) { - b(i) = bDistribution(r); - } + unsigned N; + unsigned M; + template + Model(unsigned N, unsigned M, Generator& r, T... μs) : N(N), M(M) { + Js = std::make_tuple(μs * generateRealPSpinCouplings(N, M, r)...); } - unsigned N() const { - return A.cols(); + unsigned numPs() const { + return std::tuple_size(Js); } - unsigned M() const { - return A.rows(); - } +private: + std::tuple> hamGradTensorHelper(const Vector& z, const Tensor<1>& J) const { + Tensor<3> J3(N, N, M);; + J3.setZero(); + Matrix Jz = Matrix::Zero(N, M); + Vector Jzz = Eigen::Map(J.data(), M); - Vector V(const Vector& x) const { - return A * x + b; + return {Jzz, Jz, J3}; } - Matrix dV(const Vector& x) const { - return A; + std::tuple> hamGradTensorHelper(const Vector& z, const Tensor<2>& J) const { + Tensor<3> J3(N, N, M);; + J3.setZero(); + Matrix Jz = Eigen::Map(J.data(), N, M); + Vector Jzz = z.transpose() * Jz; + + return {Jzz, Jz, J3}; } -// const Real ddV(const Vector& x) { -// return Matrix::Zero(; -// } + template + std::tuple> hamGradTensorHelper(const Vector z, const Tensor

& J) const { + Tensor<3> J3 = contractDown(J, z); + Tensor<1> zT = Eigen::TensorMap>(z.data(), N); + Tensor<2> J3zT = J3.contract(zT, ip00); + Matrix Jz = Eigen::Map(J3zT.data(), N, M); + Vector Jzz = z.transpose() * Jz; - Real H(const Vector& x) const { - return V(x).squaredNorm(); + return {Jzz, Jz, J3}; } - Vector dH(const Vector& x) const { - return dV(x).transpose() * V(x); + template + std::tuple> hamGradHessHelper(const Vector& z, const Tensor

& J, const Tensor& ...Js) const { + auto [Jzz, Jz, J3] = hamGradTensorHelper(z, J); + + Real pBang = factorial(p-1); + + Tensor<3> ddH = ((p - 1) * p / pBang) * J3; + Matrix dH = (p / pBang) * Jz; + Vector H = Jzz / pBang; + + if constexpr (sizeof...(Js) > 0) { + auto [Hs, dHs, ddHs] = hamGradHessHelper(z, Js...); + H += Hs; + dH += dHs; + ddH += ddHs; + } + + return {H, dH, ddH}; } - Matrix ddH(const Vector& x) const { - return dV(x).transpose() * dV(x); +public: + std::tuple> VdVddV(const Vector& z) const { + return std::apply([&z, this](const Tensor& ...Ks) -> std::tuple> { return hamGradHessHelper(z, Ks...); }, Js); } - Vector ∇H(const Vector& x) const { - return dH(x) - dH(x).dot(x) * x / x.squaredNorm(); + std::tuple HdHddH(const Vector& z) const { + auto [V, dV, ddV] = VdVddV(z); + + Real H = 0.5 * V.squaredNorm(); + Vector dH = dV * V; + Tensor<1> VT = Eigen::TensorMap>(V.data(), M); + Tensor<2> ddVzT = ddV.contract(VT, ip20); + Matrix ddH = Eigen::Map(ddVzT.data(), N, N) + dV * dV.transpose(); + + return {H, dH, ddH}; } - Matrix HessH(const Vector& x) const { - Matrix hess = ddH(x) - x.dot(dH(x)) * Matrix::Identity(N(), N()); - return hess - (hess * x) * x.transpose() / x.squaredNorm(); + std::tuple hamGradHess(const Vector& x) const { + auto [H, dH, ddH] = HdHddH(x); + + Vector gradH = dH - dH.dot(x) * x / N; + Matrix hessH = ddH - (dH * x.transpose() + x.dot(dH) * Matrix::Identity(N, N) + (ddH * x) * x.transpose()) / (Real)N + 2.0 * x * x.transpose(); + + return {H, gradH, hessH}; } - Vector HessSpectrum(const Vector& x) const { - Eigen::EigenSolver> eigenS(HessH(x)); + Vector HessSpectrum(const Vector& x) const { + Matrix hessH; + std::tie(std::ignore, std::ignore, hessH) = hamGradHess(x); + Eigen::EigenSolver eigenS(hessH); return eigenS.eigenvalues().real(); } }; -template -Vector normalize(const Eigen::MatrixBase& z) { - return z * sqrt((double)z.size() / (typename Derived::Scalar)(z.transpose() * z)); -} - -template -Vector findMinimum(const Model& M, const Vector& x0, Real ε) { - Vector x = x0; +template +Vector findMinimum(const Model& M, const Vector& x0, Real ε) { + Vector x = x0; Real λ = 100; - Real H = M.H(x); - Vector dH = M.dH(x); - Matrix ddH = M.ddH(x); - - Vector g = dH - x.dot(dH) * x / x.squaredNorm(); - Matrix m = ddH - (dH * x.transpose() + x.dot(dH) * Matrix::Identity(M.N(), M.N()) + (ddH * x) * x.transpose()) / x.squaredNorm() + 2.0 * x * x.transpose(); + auto [H, g, m] = M.hamGradHess(x0); while (g.norm() / x.size() > ε && λ < 1e8) { - Vector dz = (m + λ * (Matrix)abs(m.diagonal().array()).matrix().asDiagonal()).partialPivLu().solve(g); - dz -= x.dot(dz) * x / x.squaredNorm(); - Vector zNew = normalize(x - dz); + Vector dz = (m + λ * (Matrix)abs(m.diagonal().array()).matrix().asDiagonal()).partialPivLu().solve(g); + dz -= x.dot(dz) * x / M.N; + Vector zNew = normalize(x - dz); - Real HNew = M.H(zNew); - Vector dHNew = M.dH(zNew); - Matrix ddHNew = M.ddH(zNew); + auto [HNew, gNew, mNew] = M.hamGradHess(zNew); if (HNew * 1.0001 <= H) { x = zNew; H = HNew; - dH = dHNew; - ddH = ddHNew; - g = dH - x.dot(dH) * x / (Real)x.size(); - m = ddH - (dH * x.transpose() + x.dot(dH) * Matrix::Identity(x.size(), x.size()) + (ddH * x) * x.transpose()) / (Real)x.size() + 2.0 * x * x.transpose(); + g = gNew; + m = mNew; λ /= 2; } else { @@ -157,16 +171,20 @@ int main(int argc, char* argv[]) { Rng r; - Model leastSquares(σ, N, M, r.engine()); + Model<1, 2> leastSquares(N, M, r.engine(), sqrt(2) * pow(σ, 2), sqrt(2)); - Vector x = Vector::Zero(N); + Vector x = Vector::Zero(N); x(0) = sqrt(N); - std::cout << leastSquares.H(x) / N << std::endl; + double energy; + std::tie(energy, std::ignore, std::ignore) = leastSquares.hamGradHess(x); + + std::cout << energy / N << std::endl; - Vector xMin = findMinimum(leastSquares, x, 1e-12); + Vector xMin = findMinimum(leastSquares, x, 1e-12); + std::tie(energy, std::ignore, std::ignore) = leastSquares.hamGradHess(xMin); - std::cout << leastSquares.H(xMin) / N << std::endl; + std::cout << energy / N << std::endl; std::cout << leastSquares.HessSpectrum(xMin)(1) / N << std::endl; return 0; diff --git a/tensor.hpp b/tensor.hpp new file mode 100644 index 0000000..58c2434 --- /dev/null +++ b/tensor.hpp @@ -0,0 +1,108 @@ +#pragma once + +#include +#include + +#include + +#include "types.hpp" +#include "factorial.hpp" + +using Vector = Eigen::Matrix; + +using Matrix = Eigen::Matrix; + +template +using Tensor = Eigen::Tensor; + +template +using constTensor = Eigen::Tensor; + +template +Tensor

initializeJHelper(unsigned N, unsigned M, std::index_sequence) { + std::array Ns; + std::fill_n(Ns.begin(), p, N); + Ns[p-1] = M; + return Tensor

(std::get(Ns)...); +} + +template +Tensor

initializeJ(unsigned N, unsigned M) { + return initializeJHelper

(N, M, std::make_index_sequence

()); +} + +template +void setJHelper(Tensor

& J, const std::array& ind, Real val, std::index_sequence) { + J(std::get(ind)...) = val; +} + +template +void setJ(Tensor

& J, std::array ind, Real val) { + setJHelper

(J, ind, val, std::make_index_sequence

()); +} + +template +Real getJHelper(const Tensor

& J, const std::array& ind, std::index_sequence) { + return J(std::get(ind)...); +} + +template +Real getJ(const Tensor

& J, const std::array& ind) { + return getJHelper

(J, ind, std::make_index_sequence

()); +} + +template +void iterateOverHelper(Tensor

& J, + std::function&, std::array)>& f, + unsigned l, std::array is) { + if (l == 0) { + f(J, is); + } else { + for (unsigned i = 0; i < J.dimension(p - l); i++) { + std::array js = is; + js[p - l] = i; + iterateOverHelper

(J, f, l - 1, js); + } + } +} + +template +void iterateOver(Tensor

& J, std::function&, std::array)>& f) { + std::array is; + iterateOverHelper

(J, f, p, is); +} + +template +Tensor

generateCouplings(unsigned N, unsigned M, Distribution d, Generator& r) { + Tensor

J = initializeJ

(N, M); + + std::function&, std::array)> setRandom = + [&d, &r] (Tensor

& JJ, std::array ind) { + setJ

(JJ, ind, d(r)); + }; + + iterateOver

(J, setRandom); + + return J; +} + +template +Tensor

generateRealPSpinCouplings(unsigned N, unsigned M, Generator& r) { + Real σp = sqrt(factorial(p-1) / ((Real)2 * pow(N, p - 1))); + + return generateCouplings

(N, M, std::normal_distribution(0, σp), r); +} + +Tensor<3> contractDown(const Tensor<3>& J, const Vector& z) { + return J; +} + +const std::array, 1> ip00 = {Eigen::IndexPair(0, 0)}; +const std::array, 1> ip20 = {Eigen::IndexPair(2, 0)}; + +template +Tensor<3> contractDown(const Tensor& J, const Vector& z) { + Tensor<1> zT = Eigen::TensorMap>(z.data(), {z.size()}); + Tensor Jz = J.contract(zT, ip00); + return contractDown(Jz, z); +} diff --git a/types.hpp b/types.hpp new file mode 100644 index 0000000..86ccb8f --- /dev/null +++ b/types.hpp @@ -0,0 +1,3 @@ +#pragma once + +using Real = double; -- cgit v1.2.3-70-g09d2