diff options
Diffstat (limited to 'stereographic.hpp')
-rw-r--r-- | stereographic.hpp | 88 |
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}; +} |