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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
|
#include <getopt.h>
#include <iomanip>
#include "pcg-cpp/include/pcg_random.hpp"
#include "randutils/randutils.hpp"
#include "eigen/Eigen/Dense"
using Rng = randutils::random_generator<pcg32>;
using Real = double;
using Vector = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
using Matrix = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
Vector normalize(const Vector& x) {
return x * sqrt(x.size() / x.squaredNorm());
}
class QuadraticModel {
private:
Matrix J;
public:
unsigned N;
QuadraticModel(unsigned N, Rng& r) : J(N, N), N(N) {
for (unsigned i = 0; i < N; i++) {
for (unsigned j = i; j < N; j++) {
J(i, j) = r.variate<Real, std::normal_distribution>(0, 1 / sqrt(N));
if (i != j) {
J(j, i) = J(i, j);
}
}
}
}
Real H(const Vector& x) const {
return 0.5 * (J * x).dot(x);
}
Vector ∇H(const Vector& x) const {
Vector ∂H = J * x;
return ∂H - (∂H.dot(x) / x.squaredNorm()) * x;
}
};
Vector gradientDescent(const QuadraticModel& M, const Vector& x₀, Real E, Real ε = 1e-14) {
Vector xₜ = x₀;
Real Hₜ = M.H(x₀);
Real α = 1;
Real m;
Vector ∇H;
while (
∇H = (Hₜ / M.N - E) * M.∇H(xₜ) / M.N, m = ∇H.squaredNorm(),
m > ε
) {
Vector xₜ₊₁;
Real Hₜ₊₁;
while (
xₜ₊₁ = normalize(xₜ - α * ∇H), Hₜ₊₁ = M.H(xₜ₊₁),
pow(Hₜ₊₁ / M.N - E, 2) > pow(Hₜ / M.N - E, 2) - α * m
) {
α /= 2;
}
xₜ = xₜ₊₁;
Hₜ = Hₜ₊₁;
α *= 1.25;
}
return xₜ;
}
Vector randomStep(const QuadraticModel& M, const Vector& x₀, Real E, Rng& r, Real Δt = 1e-4) {
Vector η(M.N);
for (Real& ηᵢ : η) {
ηᵢ = r.variate<Real, std::normal_distribution>(0, sqrt(2 * Δt));
}
η -= η.dot(x₀) * x₀ / x₀.squaredNorm();
Vector ∇H₀ = M.∇H(x₀);
η -= η.dot(∇H₀) * ∇H₀ / ∇H₀.squaredNorm();
return gradientDescent(M, normalize(x₀ + η), E);
}
int main(int argc, char* argv[]) {
unsigned N = 10;
unsigned steps = 10;
Real E = 0;
int opt;
while ((opt = getopt(argc, argv, "N:E:n:")) != -1) {
switch (opt) {
case 'N':
N = (unsigned)atof(optarg);
break;
case 'E':
E = atof(optarg);
break;
case 'n':
steps = (unsigned)atof(optarg);
break;
default:
exit(1);
}
}
Rng r;
QuadraticModel ls(N, r);
Vector x₀ = Vector::Zero(N);
x₀(0) = sqrt(N);
x₀ = gradientDescent(ls, x₀, E);
std::cout << std::setprecision(15);
Vector x = x₀;
for (unsigned step = 0; step < steps; step++) {
x = randomStep(ls, x, E, r);
std::cout << ls.H(x) / N << " " << x.dot(x₀) / N << std::endl;
}
return 0;
}
|
