summaryrefslogtreecommitdiff
path: root/stereographic.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'stereographic.hpp')
-rw-r--r--stereographic.hpp88
1 files changed, 88 insertions, 0 deletions
diff --git a/stereographic.hpp b/stereographic.hpp
new file mode 100644
index 0000000..7a6b411
--- /dev/null
+++ b/stereographic.hpp
@@ -0,0 +1,88 @@
+
+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};
+}