Sectional curvature (#184)

* sectional curvature

* sectional curvature on power and product manifolds

* sectional curvature for product and power manifolds

* Update src/ManifoldsBase.jl

Co-authored-by: Ronny Bergmann <[email protected]>

* make tolerance a keyword argument

* change name

* update CI

---------

Co-authored-by: Ronny Bergmann <[email protected]>

.github/workflows/ci.yml

@@ -24,5 +24,8 @@ jobs:
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v4
         with:
+          token: ${{ secrets.CODECOV_TOKEN }}
+          file: ./lcov.info
+          name: codecov-umbrella
           fail_ci_if_error: false
         if: ${{ matrix.os =='ubuntu-latest' }}

NEWS.md

@@ -5,6 +5,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.15.8] 13/03/2024
+
+### Added
+
+* `sectional_curvature` , `sectional_curvature_max` and `sectional_curvature_min` functions for obtaining information about sectional curvature of a manifold.
+
 ## [0.15.7] 24/01/2024
 
 ### Fixed

Project.toml

@@ -1,7 +1,7 @@
 name = "ManifoldsBase"
 uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
 authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>", "Antoine Levitt <[email protected]>"]
-version = "0.15.7"
+version = "0.15.8"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/DefaultManifold.jl

@@ -179,6 +179,9 @@ function riemann_tensor!(::DefaultManifold, Xresult, p, X, Y, Z)
     return fill!(Xresult, 0)
 end
 
+sectional_curvature_max(::DefaultManifold) = 0.0
+sectional_curvature_min(::DefaultManifold) = 0.0
+
 Weingarten!(::DefaultManifold, Y, p, X, V) = fill!(Y, 0)
 
 zero_vector(::DefaultManifold, p) = zero(p)

src/ManifoldsBase.jl

@@ -140,6 +140,25 @@ Compute the angle between tangent vectors `X` and `Y` at point `p` from the
 function angle(M::AbstractManifold, p, X, Y)
     return acos(real(inner(M, p, X, Y)) / norm(M, p, X) / norm(M, p, Y))
 end
+
+"""
+    are_linearly_independent(M::AbstractManifold, p, X, Y)
+
+Check is vectors `X`, `Y` tangent at `p` to `M` are linearly independent.
+"""
+function are_linearly_independent(
+    M::AbstractManifold,
+    p,
+    X,
+    Y;
+    atol::Real = sqrt(eps(number_eltype(X))),
+)
+    norm_X = norm(M, p, X)
+    norm_Y = norm(M, p, Y)
+    innerXY = inner(M, p, X, Y)
+    return norm_X > atol && norm_Y > atol && !isapprox(abs(innerXY), norm_X * norm_Y)
+end
+
 """
     base_manifold(M::AbstractManifold, depth = Val(-1))
 
@@ -922,6 +941,52 @@ function riemann_tensor(M::AbstractManifold, p, X, Y, Z)
     return riemann_tensor!(M, Xresult, p, X, Y, Z)
 end
 
+@doc raw"""
+    sectional_curvature(M::AbstractManifold, p, X, Y)
+
+Compute the sectional curvature of a manifold ``\mathcal M`` at a point ``p \in \mathcal M``
+on two linearly independent tangent vectors at ``p``.
+The formula reads
+```math
+
+    \kappa_p(X, Y) = \frac{⟨R(X, Y, Y), X⟩_p}{\lVert X \rVert^2_p \lVert Y \rVert^2_p - ⟨X, Y⟩^2_p}
+
+```
+where ``R(X, Y, Y)`` is the [`riemann_tensor`](@ref) on ``\mathcal M``.
+
+# Input
+
+* `M`:   a manifold ``\mathcal M``
+* `p`:   a point ``p \in \mathcal M``
+* `X`:   a tangent vector ``X \in T_p \mathcal M``
+* `Y`:   a tangent vector ``Y \in T_p \mathcal M``
+
+"""
+function sectional_curvature(M::AbstractManifold, p, X, Y)
+    R = riemann_tensor(M, p, X, Y, Y)
+    return inner(M, p, R, X) / ((norm(M, p, X)^2 * norm(M, p, Y)^2) - (inner(M, p, X, Y)^2))
+end
+
+@doc raw"""
+    sectional_curvature_max(M::AbstractManifold)
+
+Upper bound on sectional curvature of manifold `M`. The formula reads
+```math
+\omega = \operatorname{sup}_{p\in\mathcal M, X\in T_p\mathcal M, Y\in T_p\mathcal M, ⟨X, Y⟩ ≠ 0} \kappa_p(X, Y)
+```
+"""
+sectional_curvature_max(M::AbstractManifold)
+
+@doc raw"""
+    sectional_curvature_min(M::AbstractManifold)
+
+Lower bound on sectional curvature of manifold `M`. The formula reads
+```math
+\omega = \operatorname{inf}_{p\in\mathcal M, X\in T_p\mathcal M, Y\in T_p\mathcal M, ⟨X, Y⟩ ≠ 0} \kappa_p(X, Y)
+```
+"""
+sectional_curvature_min(M::AbstractManifold)
+
 """
     size_to_tuple(::Type{S}) where S<:Tuple
 
@@ -1185,6 +1250,9 @@ export ×,
     retract!,
     riemann_tensor,
     riemann_tensor!,
+    sectional_curvature,
+    sectional_curvature_max,
+    sectional_curvature_min,
     vector_space_dimension,
     vector_transport_along,
     vector_transport_along!,

src/PowerManifold.jl

@@ -1364,6 +1364,64 @@ function riemann_tensor!(M::PowerManifoldNestedReplacing, Xresult, p, X, Y, Z)
     return Xresult
 end
 
+@doc raw"""
+    sectional_curvature(M::AbstractPowerManifold, p, X, Y)
+
+Compute the sectional curvature of a power manifold manifold ``\mathcal M`` at a point
+``p \in \mathcal M`` on two linearly independent tangent vectors at ``p``. It may be 0 for
+ if projections of `X` and `Y` on subspaces corresponding to component manifolds
+are not linearly independent.
+"""
+function sectional_curvature(M::AbstractPowerManifold, p, X, Y)
+    curvature = zero(number_eltype(X))
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        p_i = _read(M, rep_size, p, i)
+        X_i = _read(M, rep_size, X, i)
+        Y_i = _read(M, rep_size, Y, i)
+        if are_linearly_independent(M.manifold, p_i, X_i, Y_i)
+            curvature += sectional_curvature(M.manifold, p_i, X_i, Y_i)
+        end
+    end
+    return curvature
+end
+
+@doc raw"""
+    sectional_curvature_max(M::AbstractPowerManifold)
+
+Upper bound on sectional curvature of [`AbstractPowerManifold`](@ref) `M`. It is the maximum
+of sectional curvature of the wrapped manifold and 0 in case there are two or more component
+manifolds, as the sectional curvature corresponding to the plane spanned by vectors
+`(X_1, 0, ... 0)` and `(0, X_2, 0, ..., 0)` is 0.
+"""
+function sectional_curvature_max(M::AbstractPowerManifold)
+    d = prod(power_dimensions(M))
+    mscm = sectional_curvature_max(M.manifold)
+    if d > 1
+        return max(mscm, zero(mscm))
+    else
+        return mscm
+    end
+end
+
+@doc raw"""
+    sectional_curvature_min(M::AbstractPowerManifold)
+
+Lower bound on sectional curvature of [`AbstractPowerManifold`](@ref) `M`. It is the minimum
+of sectional curvature of the wrapped manifold and 0 in case there are two or more component
+manifolds, as the sectional curvature corresponding to the plane spanned by vectors
+`(X_1, 0, ... 0)` and `(0, X_2, 0, ..., 0)` is 0.
+"""
+function sectional_curvature_min(M::AbstractPowerManifold)
+    d = prod(power_dimensions(M))
+    mscm = sectional_curvature_max(M.manifold)
+    if d > 1
+        return min(mscm, zero(mscm))
+    else
+        return mscm
+    end
+end
+
 """
     set_component!(M::AbstractPowerManifold, q, p, idx...)
 

src/ProductManifold.jl

@@ -890,6 +890,63 @@ function riemann_tensor!(M::ProductManifold, Xresult, p, X, Y, Z)
     return Xresult
 end
 
+@doc raw"""
+    sectional_curvature(M::ProductManifold, p, X, Y)
+
+Compute the sectional curvature of a manifold ``\mathcal M`` at a point ``p \in \mathcal M``
+on two linearly independent tangent vectors at ``p``. It may be 0 for a product of non-flat
+manifolds if projections of `X` and `Y` on subspaces corresponding to component manifolds
+are not linearly independent.
+"""
+function sectional_curvature(M::ProductManifold, p, X, Y)
+    curvature = zero(number_eltype(X))
+    map(
+        M.manifolds,
+        submanifold_components(M, p),
+        submanifold_components(M, X),
+        submanifold_components(M, Y),
+    ) do M_i, p_i, X_i, Y_i
+        if are_linearly_independent(M_i, p_i, X_i, Y_i)
+            curvature += sectional_curvature(M_i, p_i, X_i, Y_i)
+        end
+    end
+    return curvature
+end
+
+@doc raw"""
+    sectional_curvature_max(M::ProductManifold)
+
+Upper bound on sectional curvature of [`ProductManifold`](@ref) `M`. It is the maximum of
+sectional curvatures of component manifolds and 0 in case there are two or more component
+manifolds, as the sectional curvature corresponding to the plane spanned by vectors
+`(X_1, 0)` and `(0, X_2)` is 0.
+"""
+function sectional_curvature_max(M::ProductManifold)
+    max_sc = mapreduce(sectional_curvature_max, max, M.manifolds)
+    if length(M.manifolds) > 1
+        return max(max_sc, zero(max_sc))
+    else
+        return max_sc
+    end
+end
+
+@doc raw"""
+    sectional_curvature_min(M::ProductManifold)
+
+Lower bound on sectional curvature of [`ProductManifold`](@ref) `M`. It is the minimum of
+sectional curvatures of component manifolds and 0 in case there are two or more component
+manifolds, as the sectional curvature corresponding to the plane spanned by vectors
+`(X_1, 0)` and `(0, X_2)` is 0.
+"""
+function sectional_curvature_min(M::ProductManifold)
+    min_sc = mapreduce(sectional_curvature_min, min, M.manifolds)
+    if length(M.manifolds) > 1
+        return min(min_sc, zero(min_sc))
+    else
+        return min_sc
+    end
+end
+
 """
     select_from_tuple(t::NTuple{N, Any}, positions::Val{P})
 

test/default_manifold.jl

@@ -368,6 +368,10 @@ Base.size(x::MatrixVectorTransport) = (size(x.m, 2),)
             @test riemann_tensor!(M, tv_rt, pts[1], tv1, tv2, tv1) === tv_rt
             @test tv_rt == zero(tv1)
 
+            @test sectional_curvature(M, pts[1], tv1, tv2) == 0.0
+            @test sectional_curvature_max(M) == 0.0
+            @test sectional_curvature_min(M) == 0.0
+
             q = copy(M, pts[1])
             Ts = [0.0, 1.0 / 2, 1.0]
             @testset "Geodesic interface test" begin

test/power.jl

@@ -345,6 +345,24 @@ struct TestArrayRepresentation <: AbstractPowerRepresentation end
                 [-0.1 * X, -0.1 * X],
             )
         end
+
+        @testset "sectional curvature" begin
+            Mpr = PowerManifold(M1, NestedPowerRepresentation(), 2)
+            @test sectional_curvature_max(Mpr) == 1.0
+            @test sectional_curvature_min(Mpr) == 0.0
+
+            p = [[0.0, 1.0, 0.0], [0.0, 1.0, 0.0]]
+            X1 = [[1.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
+            X2 = [[0.0, 0.0, 1.0], [0.0, 0.0, 0.0]]
+            X3 = [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]
+
+            @test sectional_curvature(Mpr, p, X1, X2) == 1.0
+            @test sectional_curvature(Mpr, p, X1, X3) == 0.0
+
+            Mss = PowerManifold(M1, NestedPowerRepresentation(), 1)
+            @test sectional_curvature_max(Mss) == 1.0
+            @test sectional_curvature_min(Mss) == 1.0
+        end
     end
 
     @testset "Type size" begin

test/product_manifold.jl

@@ -548,6 +548,25 @@ include("test_manifolds.jl")
         @test isapprox(Xresult2, Xresult)
     end
 
+    @testset "sectional curvature" begin
+        p = ArrayPartition([0.0, 1.0, 0.0], [0.0, 1.0, 0.0])
+        X1 = ArrayPartition([1.0, 0.0, 0.0], [0.0, 0.0, 0.0])
+        X2 = ArrayPartition([0.0, 0.0, 1.0], [0.0, 0.0, 0.0])
+        X3 = ArrayPartition([0.0, 0.0, 0.0], [0.0, 0.0, 1.0])
+
+        @test sectional_curvature_max(M) == 1.0
+        @test sectional_curvature_min(M) == 0.0
+
+        Mss = ProductManifold(M1, M1)
+        @test sectional_curvature_max(Mss) == 1.0
+        @test sectional_curvature_min(Mss) == 0.0
+        @test sectional_curvature(Mss, p, X1, X2) == 1.0
+        @test sectional_curvature(Mss, p, X1, X3) == 0.0
+
+        @test sectional_curvature_max(ProductManifold(M1)) == 1.0
+        @test sectional_curvature_min(ProductManifold(M1)) == 1.0
+    end
+
     @testset "× constructors" begin
         @test M1 × M2 == M
         @test M × M2 == ProductManifold(M1, M2, M2)

test/test_manifolds.jl

@@ -179,6 +179,13 @@ function ManifoldsBase.riemann_tensor!(M::TestSphere, Xresult, p, X, Y, Z)
     return Xresult
 end
 
+function ManifoldsBase.sectional_curvature_max(M::TestSphere)
+    return ifelse(manifold_dimension(M) == 1, 0.0, 1.0)
+end
+function ManifoldsBase.sectional_curvature_min(M::TestSphere)
+    return ifelse(manifold_dimension(M) == 1, 0.0, 1.0)
+end
+
 # from Absil, Mahony, Trumpf, 2013 https://sites.uclouvain.be/absil/2013-01/Weingarten_07PA_techrep.pdf
 function ManifoldsBase.Weingarten!(::TestSphere, Y, p, X, V)
     return Y .= -X * p' * V