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 +++++++++++++++++++++++++++--------------------------- 1 file changed, 80 insertions(+), 81 deletions(-) (limited to 'least_squares.cpp') 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; } -- cgit v1.2.3-70-g09d2