From 0c6efee6b0f099d71ab4daa9f79715e9dc2dc902 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 2 Apr 2024 20:20:48 +0200 Subject: Simplied implementation, specific for the quadractic models. --- least_squares.cpp | 161 +++++++++++++++++++++++++++--------------------------- tensor.hpp | 108 ------------------------------------ 2 files changed, 80 insertions(+), 189 deletions(-) delete mode 100644 tensor.hpp diff --git a/least_squares.cpp b/least_squares.cpp index f667e28..3aa9948 100644 --- a/least_squares.cpp +++ b/least_squares.cpp @@ -1,95 +1,91 @@ #include #include +#include + #include "pcg-cpp/include/pcg_random.hpp" #include "randutils/randutils.hpp" -#include "tensor.hpp" +using Rng = randutils::random_generator; -Vector normalize(const Vector& z) { - return z * sqrt((Real)z.size() / (Real)(z.transpose() * z)); -} +using Real = double; +using Vector = Eigen::Matrix; +using Matrix = Eigen::Matrix; -template -class Model { -private: - std::tuple...> Js; +/* Eigen tensor manipulations are quite annoying, especially with the need to convert other types + * into tensors beforehand. Here I overload multiplication operators to allow contraction between + * vectors and the first or last index of a tensor. + */ +class Tensor : public Eigen::Tensor { + using Eigen::Tensor::Tensor; public: - 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 numPs() const { - return std::tuple_size(Js); + Matrix operator*(const Vector& x) const { + const std::array, 1> ip20 = {Eigen::IndexPair(2, 0)}; + const Eigen::Tensor xT = Eigen::TensorMap>(x.data(), x.size()); + const Eigen::Tensor JxT = contract(xT, ip20); + return Eigen::Map(JxT.data(), dimension(0), dimension(1)); } +}; -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); - - return {Jzz, Jz, J3}; - } - - 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}; - } +Matrix operator*(const Eigen::Matrix& x, const Tensor& J) { + const std::array, 1> ip00 = {Eigen::IndexPair(0, 0)}; + const Eigen::Tensor xT = Eigen::TensorMap>(x.data(), x.size()); + const Eigen::Tensor JxT = J.contract(xT, ip00); + return Eigen::Map(JxT.data(), J.dimension(1), J.dimension(2)); +} - 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; +Vector normalize(const Vector& x) { + return x * sqrt((Real)x.size() / x.squaredNorm()); +} - return {Jzz, Jz, J3}; - } +class QuadraticModel { +private: + Tensor J; + Matrix A; + Vector b; - template - std::tuple> hamGradHessHelper(const Vector& z, const Tensor

& J, const Tensor& ...Js) const { - auto [Jzz, Jz, J3] = hamGradTensorHelper(z, J); +public: + unsigned N; + unsigned M; - Real pBang = factorial(p-1); + template + QuadraticModel(unsigned N, unsigned M, Generator& r, double μ1, double μ2, double μ3) : N(N), M(M), J(M, N, N), A(M, N), b(M) { + std::normal_distribution distribution(0, 1); - Tensor<3> ddH = ((p - 1) * p / pBang) * J3; - Matrix dH = (p / pBang) * Jz; - Vector H = Jzz / pBang; + for (unsigned i = 0; i < M; i++) { + for (unsigned j = 0; j < N; j++) { + for (unsigned k = 0; k < N; k++) { + J(i, j, k) = (2 * μ3 / N) * distribution(r); + } + } + } - if constexpr (sizeof...(Js) > 0) { - auto [Hs, dHs, ddHs] = hamGradHessHelper(z, Js...); - H += Hs; - dH += dHs; - ddH += ddHs; + for (unsigned i = 0; i < M; i++) { + for (unsigned j = 0; j < N; j++) { + A(i, j) = (μ2 / sqrt(N)) * distribution(r); + } } - return {H, dH, ddH}; + for (unsigned i = 0; i < M; i++) { + b(i) = μ1 * distribution(r); + } } -public: - std::tuple> VdVddV(const Vector& z) const { - return std::apply([&z, this](const Tensor& ...Ks) -> std::tuple> { return hamGradHessHelper(z, Ks...); }, Js); + std::tuple VdVddV(const Vector& x) const { + Matrix Jx = J * x; + Vector V1 = (A + 0.5 * Jx) * x; + Matrix dV = A + Jx; + + return {b + V1, dV, J}; } - std::tuple HdHddH(const Vector& z) const { - auto [V, dV, ddV] = VdVddV(z); + std::tuple HdHddH(const Vector& x) const { + auto [V, dV, ddV] = VdVddV(x); 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(); + Vector dH = V.transpose() * dV; + Matrix ddH = V.transpose() * ddV + dV.transpose() * dV; return {H, dH, ddH}; } @@ -103,7 +99,7 @@ public: return {H, gradH, hessH}; } - Vector HessSpectrum(const Vector& x) const { + Vector spectrum(const Vector& x) const { Matrix hessH; std::tie(std::ignore, std::ignore, hessH) = hamGradHess(x); Eigen::EigenSolver eigenS(hessH); @@ -111,24 +107,22 @@ public: } }; -template -Vector findMinimum(const Model& M, const Vector& x0, Real ε) { +Vector findMinimum(const QuadraticModel& M, const Vector& x0, Real ε) { Vector x = x0; Real λ = 100; 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 / M.N; - Vector zNew = normalize(x - dz); + Vector dx = (m + λ * (Matrix)abs(m.diagonal().array()).matrix().asDiagonal()).partialPivLu().solve(g); + dx -= x.dot(dx) * x / M.N; + Vector xNew = normalize(x - dx); - auto [HNew, gNew, mNew] = M.hamGradHess(zNew); + auto [HNew, gNew, mNew] = M.hamGradHess(xNew); if (HNew * 1.0001 <= H) { - x = zNew; + x = xNew; H = HNew; - g = gNew; m = mNew; @@ -141,17 +135,16 @@ Vector findMinimum(const Model& M, const Vector& x0, Real ε) { return x; } -using Rng = randutils::random_generator; -using Real = double; - int main(int argc, char* argv[]) { unsigned N = 10; Real α = 1; Real σ = 1; + Real A = 1; + Real J = 1; int opt; - while ((opt = getopt(argc, argv, "N:a:s:")) != -1) { + while ((opt = getopt(argc, argv, "N:a:s:A:J:")) != -1) { switch (opt) { case 'N': N = (unsigned)atof(optarg); @@ -162,6 +155,12 @@ int main(int argc, char* argv[]) { case 's': σ = atof(optarg); break; + case 'A': + A = atof(optarg); + break; + case 'J': + J = atof(optarg); + break; default: exit(1); } @@ -171,7 +170,7 @@ int main(int argc, char* argv[]) { Rng r; - Model<1, 2> leastSquares(N, M, r.engine(), σ, 1); + QuadraticModel leastSquares(N, M, r.engine(), σ, A, J); Vector x = Vector::Zero(N); x(0) = sqrt(N); @@ -185,7 +184,7 @@ int main(int argc, char* argv[]) { std::tie(energy, std::ignore, std::ignore) = leastSquares.hamGradHess(xMin); std::cout << energy / N << std::endl; - std::cout << leastSquares.HessSpectrum(xMin)(1) / N << std::endl; + std::cout << leastSquares.spectrum(xMin)(1) / N << std::endl; return 0; } diff --git a/tensor.hpp b/tensor.hpp deleted file mode 100644 index 7f98d02..0000000 --- a/tensor.hpp +++ /dev/null @@ -1,108 +0,0 @@ -#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)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); -} -- cgit v1.2.3-70-g09d2