Added tensor support.

4 files changed, 216 insertions, 78 deletions

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 <eigen3/Eigen/Dense>
-#include <random>
 #include <getopt.h>
 
 #include "pcg-cpp/include/pcg_random.hpp"
 #include "randutils/randutils.hpp"
 
+#include "tensor.hpp"
 
-template <class Real>
-using Vector = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
-
-template <class Real>
-using Matrix = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
+Vector normalize(const Vector& z) {
+  return z * sqrt((Real)z.size() / (Real)(z.transpose() * z));
+}
 
-template <class Real>
+template <int... ps>
 class Model {
 private:
-  Matrix<Real> A;
-  Vector<Real> b;
+  std::tuple<Tensor<ps>...> Js;
 
 public:
-  template <class Generator>
-  Model(Real σ, unsigned N, unsigned M, Generator& r) : A(M, N), b(M) {
-    std::normal_distribution<Real> 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<Real> bDistribution(0, σ);
-
-    for (unsigned i = 0; i < M; i++) {
-      b(i) = bDistribution(r);
-    }
+  unsigned N;
+  unsigned M;
+  template <class Generator, typename... T>
+  Model(unsigned N, unsigned M, Generator& r, T... μs) : N(N), M(M) {
+    Js = std::make_tuple(μs * generateRealPSpinCouplings<ps>(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<Vector, Matrix, Tensor<3>> 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<const Vector>(J.data(), M);
 
-  Vector<Real> V(const Vector<Real>& x) const {
-    return A * x + b;
+    return {Jzz, Jz, J3};
   }
 
-  Matrix<Real> dV(const Vector<Real>& x) const {
-    return A;
+  std::tuple<Vector, Matrix, Tensor<3>> hamGradTensorHelper(const Vector& z, const Tensor<2>& J) const {
+    Tensor<3> J3(N, N, M);;
+    J3.setZero();
+    Matrix Jz = Eigen::Map<const Matrix>(J.data(), N, M);
+    Vector Jzz = z.transpose() * Jz;
+
+    return {Jzz, Jz, J3};
   }
 
-//  const Real ddV(const Vector<Real>& x) {
-//    return Matrix::Zero(;
-//  }
+  template <int p>
+  std::tuple<Vector, Matrix, Tensor<3>> hamGradTensorHelper(const Vector z, const Tensor<p>& J) const {
+    Tensor<3> J3 = contractDown(J, z);
+    Tensor<1> zT = Eigen::TensorMap<constTensor<1>>(z.data(), N);
+    Tensor<2> J3zT = J3.contract(zT, ip00);
+    Matrix Jz = Eigen::Map<const Matrix>(J3zT.data(), N, M);
+    Vector Jzz = z.transpose() * Jz;
 
-  Real H(const Vector<Real>& x) const {
-    return V(x).squaredNorm();
+    return {Jzz, Jz, J3};
   }
 
-  Vector<Real> dH(const Vector<Real>& x) const {
-    return dV(x).transpose() * V(x);
+  template <int p, int... qs>
+  std::tuple<Vector, Matrix, Tensor<3>> hamGradHessHelper(const Vector& z, const Tensor<p>& J, const Tensor<qs>& ...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<Real> ddH(const Vector<Real>& x) const {
-    return dV(x).transpose() * dV(x);
+public:
+  std::tuple<Vector, Matrix, Tensor<3>> VdVddV(const Vector& z) const {
+    return std::apply([&z, this](const Tensor<ps>& ...Ks) -> std::tuple<Vector, Matrix, Tensor<3>> { return hamGradHessHelper(z, Ks...); }, Js);
   }
 
-  Vector<Real> ∇H(const Vector<Real>& x) const {
-    return dH(x) - dH(x).dot(x) * x / x.squaredNorm();
+  std::tuple<Real, Vector, Matrix> 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<constTensor<1>>(V.data(), M);
+    Tensor<2> ddVzT = ddV.contract(VT, ip20);
+    Matrix ddH = Eigen::Map<const Matrix>(ddVzT.data(), N, N) + dV * dV.transpose();
+
+    return {H, dH, ddH};
   }
 
-  Matrix<Real> HessH(const Vector<Real>& x) const {
-    Matrix<Real> hess = ddH(x) - x.dot(dH(x)) * Matrix<Real>::Identity(N(), N());
-    return hess - (hess * x) * x.transpose() / x.squaredNorm();
+  std::tuple<Real, Vector, Matrix> 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<Real> HessSpectrum(const Vector<Real>& x) const {
-    Eigen::EigenSolver<Matrix<Real>> eigenS(HessH(x));
+  Vector HessSpectrum(const Vector& x) const {
+    Matrix hessH;
+    std::tie(std::ignore, std::ignore, hessH) = hamGradHess(x);
+    Eigen::EigenSolver<Matrix> eigenS(hessH);
     return eigenS.eigenvalues().real();
   }
 };
 
-template <typename Derived>
-Vector<typename Derived::Scalar> normalize(const Eigen::MatrixBase<Derived>& z) {
-  return z * sqrt((double)z.size() / (typename Derived::Scalar)(z.transpose() * z));
-}
-
-template <class Real>
-Vector<Real> findMinimum(const Model<Real>& M, const Vector<Real>& x0, Real ε) {
-  Vector<Real> x = x0;
+template <int ...ps>
+Vector findMinimum(const Model<ps...>& M, const Vector& x0, Real ε) {
+  Vector x = x0;
   Real λ = 100;
 
-  Real H = M.H(x);
-  Vector<Real> dH = M.dH(x);
-  Matrix<Real> ddH = M.ddH(x);
-
-  Vector<Real> g = dH - x.dot(dH) * x / x.squaredNorm();
-  Matrix<Real> m = ddH - (dH * x.transpose() + x.dot(dH) * Matrix<Real>::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<Real> dz = (m + λ * (Matrix<Real>)abs(m.diagonal().array()).matrix().asDiagonal()).partialPivLu().solve(g);
-    dz -= x.dot(dz) * x / x.squaredNorm();
-    Vector<Real> 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<Real> dHNew = M.dH(zNew);
-    Matrix<Real> 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<Real>::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<Real> leastSquares(σ, N, M, r.engine());
+  Model<1, 2> leastSquares(N, M, r.engine(), sqrt(2) * pow(σ, 2), sqrt(2));
 
-  Vector<Real> x = Vector<Real>::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<Real> 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 <array>
+#include <functional>
+
+#include <eigen3/unsupported/Eigen/CXX11/Tensor>
+
+#include "types.hpp"
+#include "factorial.hpp"
+
+using Vector = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
+
+using Matrix = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
+
+template <int p>
+using Tensor = Eigen::Tensor<Real, p>;
+
+template <int p>
+using constTensor = Eigen::Tensor<const Real, p>;
+
+template <int p, std::size_t... Indices>
+Tensor<p> initializeJHelper(unsigned N, unsigned M, std::index_sequence<Indices...>) {
+  std::array<unsigned, p> Ns;
+  std::fill_n(Ns.begin(), p, N);
+  Ns[p-1] = M;
+  return Tensor<p>(std::get<Indices>(Ns)...);
+}
+
+template <int p>
+Tensor<p> initializeJ(unsigned N, unsigned M) {
+  return initializeJHelper<p>(N, M, std::make_index_sequence<p>());
+}
+
+template <int p, std::size_t... Indices>
+void setJHelper(Tensor<p>& J, const std::array<unsigned, p>& ind, Real val, std::index_sequence<Indices...>) {
+  J(std::get<Indices>(ind)...) = val;
+}
+
+template <int p>
+void setJ(Tensor<p>& J, std::array<unsigned, p> ind, Real val) {
+  setJHelper<p>(J, ind, val, std::make_index_sequence<p>());
+}
+
+template <int p, std::size_t... Indices>
+Real getJHelper(const Tensor<p>& J, const std::array<unsigned, p>& ind, std::index_sequence<Indices...>) {
+  return J(std::get<Indices>(ind)...);
+}
+
+template <int p>
+Real getJ(const Tensor<p>& J, const std::array<unsigned, p>& ind) {
+  return getJHelper<p>(J, ind, std::make_index_sequence<p>());
+}
+
+template <int p>
+void iterateOverHelper(Tensor<p>& J,
+    std::function<void(Tensor<p>&, std::array<unsigned, p>)>& f,
+    unsigned l, std::array<unsigned, p> is) {
+  if (l == 0) {
+    f(J, is);
+  } else {
+    for (unsigned i = 0; i < J.dimension(p - l); i++) {
+      std::array<unsigned, p> js = is;
+      js[p - l] = i;
+      iterateOverHelper<p>(J, f, l - 1, js);
+    }
+  }
+}
+
+template <int p>
+void iterateOver(Tensor<p>& J, std::function<void(Tensor<p>&, std::array<unsigned, p>)>& f) {
+  std::array<unsigned, p> is;
+  iterateOverHelper<p>(J, f, p, is);
+}
+
+template <int p, class Distribution, class Generator>
+Tensor<p> generateCouplings(unsigned N, unsigned M, Distribution d, Generator& r) {
+  Tensor<p> J = initializeJ<p>(N, M);
+
+  std::function<void(Tensor<p>&, std::array<unsigned, p>)> setRandom =
+    [&d, &r] (Tensor<p>& JJ, std::array<unsigned, p> ind) {
+      setJ<p>(JJ, ind, d(r));
+    };
+
+  iterateOver<p>(J, setRandom);
+
+  return J;
+}
+
+template <int p, class Generator>
+Tensor<p> generateRealPSpinCouplings(unsigned N, unsigned M, Generator& r) {
+  Real σp = sqrt(factorial(p-1) / ((Real)2 * pow(N, p - 1)));
+
+  return generateCouplings<p>(N, M, std::normal_distribution<Real>(0, σp), r);
+}
+
+Tensor<3> contractDown(const Tensor<3>& J, const Vector& z) {
+  return J;
+}
+
+const std::array<Eigen::IndexPair<int>, 1> ip00 = {Eigen::IndexPair<int>(0, 0)};
+const std::array<Eigen::IndexPair<int>, 1> ip20 = {Eigen::IndexPair<int>(2, 0)};
+
+template <int r>
+Tensor<3> contractDown(const Tensor<r>& J, const Vector& z) {
+  Tensor<1> zT = Eigen::TensorMap<constTensor<1>>(z.data(), {z.size()});
+  Tensor<r - 1> 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;