summaryrefslogtreecommitdiff
path: root/dynamics.hpp
blob: 7d04b1a8f3f0a3daf69f8549c734d14f36006892 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#pragma once

#include <exception>
#include <eigen3/Eigen/LU>
#include <random>

#include "p-spin.hpp"

template <class Real, int p>
Vector<Real> findMinimum(const Tensor<Real, p>& J, const Vector<Real>& z0, Real ε) {
  Vector<Real> z = z0;
  Real λ = 10;

  auto [H, dH, ddH] = hamGradHess(J, z);

  Vector<Real> g = dH - z.dot(dH) * z / (Real)z.size();
  Matrix<Real> M = ddH - (dH * z.transpose() + z.dot(dH) * Matrix<Real>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose();

  while (g.norm() / z.size() > ε && λ < 1e8) {
    Vector<Real> dz = (M + λ * (Matrix<Real>)abs(M.diagonal().array()).matrix().asDiagonal()).partialPivLu().solve(g);
    dz -= z.dot(dz) * z / (Real)z.size();
    Vector<Real> zNew = normalize(z - dz);

    auto [HNew, dHNew, ddHNew] = hamGradHess(J, zNew);

    if (HNew * 1.0001 <= H) {
      z = zNew;
      H = HNew;
      dH = dHNew;
      ddH = ddHNew;

      g = dH - z.dot(dH) * z / (Real)z.size();
      M = ddH - (dH * z.transpose() + z.dot(dH) * Matrix<Real>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose();

      λ /= 2;
    } else {
      λ *= 1.5;
    }
  }

  return z;
}

template <class Real, class Scalar, int p>
Vector<Scalar> findSaddle(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, Real ε) {
  Vector<Scalar> z = z0;

  Vector<Scalar> dH;
  Matrix<Scalar> ddH;
  std::tie(std::ignore, dH, ddH) = hamGradHess(J, z);

  Scalar zz = z.transpose() * z;
  Vector<Scalar> g = dH - (z.transpose() * dH) * z / (Real)z.size() + z * (zz - (Real)z.size());
  Matrix<Scalar> M = ddH - (dH * z.transpose() + (z.transpose() * dH) * Matrix<Scalar>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size()  + Matrix<Scalar>::Identity(z.size(), z.size()) * (zz - (Real)z.size()) + 2.0 * z * z.transpose();

  while (g.norm() / z.size() > ε) {
    Vector<Scalar> dz = M.partialPivLu().solve(g);
    dz -= z.conjugate().dot(dz) / z.squaredNorm() * z.conjugate();
    z = normalize(z - dz);

    std::tie(std::ignore, dH, ddH) = hamGradHess(J, z);

    zz = z.transpose() * z;
    g = dH - (z.transpose() * dH) * z / (Real)z.size() + z * (zz - (Real)z.size());
    M = ddH - (dH * z.transpose() + (z.transpose() * dH) * Matrix<Scalar>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size()  + Matrix<Scalar>::Identity(z.size(), z.size()) * (zz - (Real)z.size()) + 2.0 * z * z.transpose();
  }

  return z;
}

template <class Scalar, class Distribution, class Generator>
Vector<Scalar> randomVector(unsigned N, Distribution d, Generator& r) {
  Vector<Scalar> z(N);

  for (unsigned i = 0; i < N; i++) {
    z(i) = d(r);
  }

  return z;
}

template <class Real, int p, class Distribution, class Generator>
Vector<Real> randomMinimum(const Tensor<Real, p>& J, Distribution d, Generator& r, Real ε) {
  Vector<Real> zSaddle;
  bool foundSaddle = false;

  while (!foundSaddle) {
    Vector<Real> z0 = normalize(randomVector<Real>(J.dimension(0), d, r.engine()));

    try {
      zSaddle = findMinimum(J, z0, ε);
      foundSaddle = true;
    } catch (std::exception& e) {}
  }

  return zSaddle;
}

template <class Real, class Scalar, int p, class Distribution, class Generator>
Vector<Scalar> randomSaddle(const Tensor<Scalar, p>& J, Distribution d, Generator& r, Real ε) {
  Vector<Scalar> zSaddle;
  bool foundSaddle = false;

  while (!foundSaddle) {
    Vector<Scalar> z0 = normalize(randomVector<Scalar>(J.dimension(0), d, r.engine()));

    try {
      zSaddle = findSaddle(J, z0, ε);
      foundSaddle = true;
    } catch (std::exception& e) {}
  }

  return zSaddle;
}