Exploring curvature without coordinates

This part of documentation covers exploration of curvature of manifolds $\mathcal{M}$. There are multiple ways to describe curvature: Christoffel symbols, Riemann tensor, Ricci tensor, sectional curvature, and many other. They are usually considered only in coordinates but there is a way to demonstrate curvature in coordinate-free way.

Sectional curvature matrix

Curvature of a manifold can be explored using the sectional_curvature_matrix function. Note that Riemann tensor and sectional curvature are equivalently full specifications of curvature in a manifold, see [CE08], Eq. (1.12).

Let’s take the SymmetricPositiveDefinite manifold as our first example. It has nonpositive sectional curvature:

using Manifolds
using LinearAlgebra
M = SymmetricPositiveDefinite(3)
p = rand(M)
cm = sectional_curvature_matrix(M, p, DefaultOrthonormalBasis())

6×6 Matrix{Float64}:
  0.0          -0.25         -0.25         …  -2.29564e-19   9.57203e-21
 -0.25          0.0          -0.125           -0.125         9.60226e-21
 -0.25         -0.125         0.0             -0.125        -0.25
 -8.11985e-23  -0.25         -1.36129e-23     -0.25          6.99265e-23
 -2.29564e-19  -0.125        -0.125            0.0          -0.25
  9.57203e-21   9.60226e-21  -0.25         …  -0.25          0.0

We can verify that the curvature is consistent with an approximation based on the Bertrand–Diguet–Puiseux theorem, which relies only on an ONB, exponential map and distance calculation:

cm_bdp = Manifolds.estimated_sectional_curvature_matrix(M, p, DefaultOrthonormalBasis(); r=1e-3, N_pts=100000)
println(norm(cm - cm_bdp))

0.004362482234959643

This approximation converges quite slowly with N_pts and is prone to numerical errors at low values of r and large values of N_pts.

You can also take the vectors from the basis and see what kind of planes they correspond to. It may be easier to see for the identity matrix as the base point.

p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
V = get_vectors(M, p, get_basis(M, p, DefaultOrthonormalBasis()))
cm = sectional_curvature_matrix(M, p, DefaultOrthonormalBasis())
for X in V
    println(exp(M, p, X))
end

[2.718281828459045 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
[1.260591836521356 0.7675231451261162 0.0; 0.7675231451261162 1.2605918365213566 0.0; 0.0 0.0 1.0]
[1.260591836521356 0.0 0.7675231451261162; 0.0 1.0 0.0; 0.7675231451261162 0.0 1.2605918365213566]
[1.0 0.0 0.0; 0.0 2.718281828459045 0.0; 0.0 0.0 1.0]
[1.0 0.0 0.0; 0.0 1.260591836521356 0.7675231451261162; 0.0 0.7675231451261162 1.2605918365213566]
[1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 2.718281828459045]

The flat planes correspond to directions where the matrix changes independently. In other cases sectional curvature indicates hyperbolic characteristic of a submanifold.

Sectional curvature can be either larger or smaller than entries in the matrix on other planes. Consider for example the manifold of rotation matrices in four dimensions, and a function that computes plane of maximum curvature using random search.

function max_curvature(M::AbstractManifold, p)
    mc = -Inf
    X = zero_vector(M, p)
    Y = zero_vector(M, p)
    for _ in 1:10000
        X_c = rand(M; vector_at=p)
        Y_c = rand(M; vector_at=p)
        sc = sectional_curvature(M, p, X_c, Y_c)
        if sc > mc
            mc = sc
            X .= X_c
            Y .= Y_c
        end
    end
    return mc, X, Y
end

M = Rotations(4)
p = Matrix(I(4) * 1.0)
println(sectional_curvature_matrix(M, p, DefaultOrthonormalBasis()))
mc, X, Y = max_curvature(M, p)
println(mc)
println(X)
println(Y)

[0.0 0.12500000000000003 0.12500000000000003 0.0 0.12500000000000003 0.12500000000000003; 0.12500000000000003 0.0 0.12500000000000003 0.12500000000000003 0.0 0.12500000000000003; 0.12500000000000003 0.12500000000000003 0.0 0.12500000000000003 0.12500000000000003 0.0; 0.0 0.12500000000000003 0.12500000000000003 0.0 0.12500000000000003 0.12500000000000003; 0.12500000000000003 0.0 0.12500000000000003 0.12500000000000003 0.0 0.12500000000000003; 0.12500000000000003 0.12500000000000003 0.0 0.12500000000000003 0.12500000000000003 0.0]
0.24526794169749644
[0.0 1.9811573509348712 0.03747356881939395 -0.09723213240472224; -1.9811573509348712 0.0 0.1688225140838157 -0.375367925089595; -0.03747356881939395 -0.1688225140838157 0.0 -2.1556650349434903; 0.09723213240472224 0.375367925089595 2.1556650349434903 0.0]
[0.0 0.4204378372126399 -1.7067515919002116 0.25433452471379986; -0.4204378372126399 0.0 -0.026561719645489295 -1.7494201866754369; 1.7067515919002116 0.026561719645489295 0.0 -0.264237403993821; -0.25433452471379986 1.7494201866754369 0.264237403993821 0.0]

In the planes corresponding to orthonormal basis, the maximum sectional curvature is 0.125 but the true upper bound is 0.25.

Literature