summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2021-11-09 00:44:01 +0100
committerJaron Kent-Dobias <jaron@kent-dobias.com>2021-11-09 00:44:01 +0100
commitd448af5010b664025c816dc2c6e383ac7bea3491 (patch)
tree68ef15565f7058d6029a0cce2e7d79724ff8cc46
parent1369382dcab07078e48e929e76802cdc4a1ceeca (diff)
downloadcode-d448af5010b664025c816dc2c6e383ac7bea3491.tar.gz
code-d448af5010b664025c816dc2c6e383ac7bea3491.tar.bz2
code-d448af5010b664025c816dc2c6e383ac7bea3491.zip
Generalized energy function to multiple tensor arguments in anticipation of the mixed p-spin.
-rw-r--r--p-spin.hpp39
1 files changed, 13 insertions, 26 deletions
diff --git a/p-spin.hpp b/p-spin.hpp
index 15d2525..6111b75 100644
--- a/p-spin.hpp
+++ b/p-spin.hpp
@@ -12,19 +12,26 @@ Vector<typename Derived::Scalar> normalize(const Eigen::MatrixBase<Derived>& z)
return z * sqrt((Real)z.size() / (typename Derived::Scalar)(z.transpose() * z));
}
-template <class Scalar, int p>
-std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>> hamGradHess(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
+template <class Scalar, int p, int ... ps>
+std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>> hamGradHess(const Tensor<Scalar, p>& J, const Tensor<Scalar, ps>& ... Js, const Vector<Scalar>& z) {
Matrix<Scalar> Jz = contractDown(J, z); // Contracts J into p - 2 copies of z.
Vector<Scalar> Jzz = Jz * z;
Scalar Jzzz = Jzz.transpose() * z;
Real pBang = factorial(p);
- Matrix<Scalar> hessian = ((p - 1) * p / pBang) * Jz;
- Vector<Scalar> gradient = (p / pBang) * Jzz;
- Scalar hamiltonian = Jzzz / pBang;
+ Matrix<Scalar> ddH = ((p - 1) * p / pBang) * Jz;
+ Vector<Scalar> dH = (p / pBang) * Jzz;
+ Scalar H = Jzzz / pBang;
+
+ if constexpr (sizeof...(Js) > 0) {
+ auto [Hs, dHs, ddHs] = hamGradHess(Js..., z);
+ H += Hs;
+ dH += dHs;
+ ddH += ddHs;
+ }
- return {hamiltonian, gradient, hessian};
+ return {H, dH, ddH};
}
template <class Scalar, int p>
@@ -71,26 +78,6 @@ Vector<Scalar> zDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) {
return -dH.conjugate() + (dH.dot(z) / z.squaredNorm()) * z.conjugate();
}
-template <class Scalar, int p>
-std::tuple<Real, Vector<Scalar>> WdW(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
- Vector<Scalar> dH;
- Matrix<Scalar> ddH;
- std::tie(std::ignore, dH, ddH) = hamGradHess(J, z);
-
- Vector<Scalar> dzdt = zDot(z, dH);
-
- Real a = z.squaredNorm();
- Scalar A = (Scalar)(z.transpose() * dzdt) / a;
- Scalar B = dH.dot(z) / a;
-
- Real W = dzdt.squaredNorm();
- Vector<Scalar> dW = ddH * (dzdt - A * z.conjugate())
- + 2 * (conj(A) * B * z).real()
- - conj(B) * dzdt - conj(A) * dH.conjugate();
-
- return {W, dW};
-}
-
template <class Scalar>
Matrix<Scalar> dzDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) {
Real z² = z.squaredNorm();