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;
}
|
