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
|
Vector stereographicToEuclidean(const Vector& ζ) {
unsigned N = ζ.size() + 1;
Vector z(N);
Scalar a = ζ.dot(ζ);
std::complex<double> b = 2.0 * sqrt(N) / (1.0 + a);
for (unsigned i = 0; i < N - 1; i++) {
z(i) = b * ζ(i);
}
z(N - 1) = b * (a - 1.0) / 2.0;
return z;
}
Vector euclideanToStereographic(const Vector& z) {
unsigned N = z.size();
Vector ζ(N - 1);
for (unsigned i = 0; i < N - 1; i++) {
ζ(i) = z(i) / (sqrt(N) - z(N - 1));
}
return ζ;
}
Matrix stereographicJacobian(const Vector& ζ) {
unsigned N = ζ.size() + 1;
Matrix J(N - 1, N);
Scalar b = ζ.dot(ζ);
for (unsigned i = 0; i < N - 1; i++) {
for (unsigned j = 0; j < N - 1; j++) {
J(i, j) = -4.0 * ζ(i) * ζ(j);
if (i == j) {
J(i, j) += 2.0 * (1.0 + b);
}
J(i, j) *= sqrt(N) / pow(1.0 + b, 2);
}
J(i, N - 1) = 4.0 * sqrt(N) * ζ(i) / pow(1.0 + b, 2);
}
return J;
}
Matrix stereographicMetric(const Vector& ζ) {
unsigned N = ζ.size();
Matrix g(N, N);
double a = ζ.cwiseAbs2().sum();
Scalar b = ζ.dot(ζ);
for (unsigned i = 0; i < N; i++) {
for (unsigned j = 0; j < N; j++) {
g(i, j) = 16.0 * (std::conj(ζ(i)) * ζ(j) * (1.0 + a) -
std::real(std::conj(ζ(i) * ζ(j)) * (1.0 + b)));
if (i == j) {
g(i, j) += 4.0 * std::abs(1.0 + b);
}
g(i, j) *= (N + 1) / std::norm(pow(b, 2));
}
}
return g;
}
std::tuple<std::complex<double>, Vector, Matrix> stereographicHamGradHess(const Tensor3& J,
const Vector& ζ) {
Vector grad;
Matrix hess;
Matrix jacobian = stereographicJacobian(ζ);
Vector z = stereographicToEuclidean(ζ);
std::complex<double> hamiltonian;
Vector gradZ;
Matrix hessZ;
std::tie(hamiltonian, gradZ, hessZ) = hamGradHess(J, z);
Matrix g = stereographicMetric(ζ);
Matrix gj = g.llt().solve(jacobian);
grad = gj * gradZ;
hess = gj * hessZ * gj.transpose();
return {hamiltonian, grad, hess};
}
|
