Sectional curvature (#711)

* Sectional curvature

* formatting

* update CI settings

* fix again

* address review

* Upper-bounding NonlinearSolve

* Add NonlinearSolve as a direct dep

* tweaking

* more curvature

* add atol to failing test

* bump tolerance

* upper bound for PythonCall too

* change default vector transport method for some unitary matrix manifolds

* fix news

.github/workflows/documenter.yml

@@ -14,7 +14,7 @@ jobs:
       - uses: actions/checkout@v4
       - uses: quarto-dev/quarto-actions/setup@v2
         with:
-          version: 1.3.361
+          version: "1.4.551"
       - uses: julia-actions/setup-julia@latest
         with:
           version: "1.10"

NEWS.md

@@ -10,6 +10,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 
 * using `DocumenterInterLinks` for links to other Julia packages documentation.
+* Implementation of `sectional_curvature`, `sectional_curvature_min` and `sectional_curvature_max` for several manifolds.
+* `sectional_curvature_matrix` function and a tutorial on coordinate-free curvature.
+
+### Changed
+
+* `default_vector_transport_method` for `GeneralUnitaryMatrices` other than `Rotations` was changed to `ProjectionTransport`.
 
 ### Fixed
 

Project.toml

@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>", "Antoine Levitt <[email protected]>"]
-version = "0.9.14"
+version = "0.9.15"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -16,6 +16,7 @@ Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 MatrixEquations = "99c1a7ee-ab34-5fd5-8076-27c950a045f4"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
+PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
@@ -52,12 +53,13 @@ HybridArrays = "0.4"
 Kronecker = "0.4, 0.5"
 LinearAlgebra = "1.6"
 ManifoldDiff = "0.3.7"
-ManifoldsBase = "0.15.6"
+ManifoldsBase = "0.15.8"
 Markdown = "1.6"
 MatrixEquations = "2.2"
-NonlinearSolve = " < 3.8"
+NonlinearSolve = "1 - 3.7.3"
 OrdinaryDiffEq = "6.31"
 Plots = "1"
+PythonCall = "0.9 - 0.9.15"
 Quaternions = "0.5, 0.6, 0.7"
 Random = "1.6"
 RecipesBase = "1.1"

docs/make.jl

@@ -104,6 +104,7 @@ makedocs(;
             "work in charts" => "tutorials/working-in-charts.md",
             "perform Hand gesture analysis" => "tutorials/hand-gestures.md",
             "integrate on manifolds and handle probability densities" => "tutorials/integration.md",
+            "explore curvature without coordinates" => "tutorials/exploring-curvature.md",
         ],
         "Manifolds" => [
             "Basic manifolds" => [

docs/src/features/utilities.md

@@ -8,18 +8,12 @@ The following terms introduce a nicer notation for some operations, for example
 in
 ````
 
-## Fallback for the exponential map: Solving the corresponding ODE
-
-When additionally loading [`NLSolve.jl`](https://github.com/JuliaNLSolvers/NLsolve.jl) the following fallback for the exponential map is available.
+## Public documentation
 
-```@autodocs
-Modules = [Manifolds]
-Pages = ["nlsolve.jl"]
-Order = [:type, :function]
+```@docs
+sectional_curvature_matrix
 ```
 
-## Public documentation
-
 ### Specific exception types
 
 For some manifolds it is useful to keep an extra index, at which point on the manifold, the error occurred as well as to collect all errors that occurred on a manifold. This page contains the manifold-specific error messages this package introduces.

docs/src/misc/notation.md

@@ -49,6 +49,7 @@ Within the documented functions, the utf8 symbols are used whenever possible, as
 | ``\mathcal P_{t_1\gets t_0}^cX`` | parallel transport along the curve ``c``| ``\mathcal P^cX=\mathcal P_{1\gets 0}^cX`` | of the vector ``X`` from ``p=c(0)`` to ``c(1)``
 | ``p`` | a point on ``\mathcal M`` | ``p_1, p_2, \ldots,q`` | for 3 points one might use ``x,y,z`` |
 | ``\operatorname{retr}_pX``|  a retraction | |
+| ``\kappa_p(X, Y)`` | sectional curvature | |
 | ``ξ`` | a set of tangent vectors | ``\{X_1,\ldots,X_n\}`` | |
 | ``J_{2n} \in ℝ^{2n×2n}`` | the [`SymplecticElement`](@ref) | | |
 | ``T_p \mathcal M`` | the tangent space at ``p`` | | |

docs/src/references.bib

@@ -194,6 +194,16 @@ @article{LeBrignantPuechmorel:2019
 #
 # C
 # ----------------------------------------------------------------------------------------
+
+@book{CheegerEbin:2008,
+	address = {Providence, R.I},
+	title = {Comparison {Theorems} in {Riemannian} {Geometry}},
+	isbn = {978-0-8218-4417-5},
+	publisher = {American Mathematical Society},
+	author = {Cheeger, Jeffrey and Ebin, David G.},
+	month = may,
+	year = {2008},
+}
 @book{Chikuse:2003,
     DOI    = {10.1007/978-0-387-21540-2},
     YEAR   = {2003},
@@ -341,6 +351,18 @@ @article{GaoSonAbsilStykel:2021
     TITLE   = {Riemannian Optimization on the Symplectic Stiefel Manifold},
     JOURNAL = {SIAM Journal on Optimization}
 }
+@article{Ge:2014,
+	title = {{DDVV}-type inequality for skew-symmetric matrices and {Simons}-type inequality for {Riemannian} submersions},
+	volume = {251},
+	issn = {0001-8708},
+	url = {https://www.sciencedirect.com/science/article/pii/S000187081300385X},
+	doi = {10.1016/j.aim.2013.10.010},
+	journal = {Advances in Mathematics},
+	author = {Ge, Jianquan},
+	month = jan,
+	year = {2014},
+	pages = {62--86},
+}
 @inproceedings{Giles:2008,
 	address = {Berlin, Heidelberg},
 	series = {Lecture {Notes} in {Computational} {Science} and {Engineering}},

src/Manifolds.jl

@@ -146,6 +146,9 @@ import ManifoldsBase:
     retract_softmax!,
     riemann_tensor,
     riemann_tensor!,
+    sectional_curvature,
+    sectional_curvature_max,
+    sectional_curvature_min,
     set_component!,
     submanifold,
     submanifold_component,
@@ -906,6 +909,10 @@ export ×,
     riemannian_Hessian!,
     riemann_tensor,
     riemann_tensor!,
+    sectional_curvature,
+    sectional_curvature_matrix,
+    sectional_curvature_max,
+    sectional_curvature_min,
     set_component!,
     sharp,
     sharp!,

src/manifolds/Euclidean.jl

@@ -753,6 +753,33 @@ function riemann_tensor!(::Euclidean, Xresult, p, X, Y, Z)
     return fill!(Xresult, 0)
 end
 
+@doc raw"""
+    sectional_curvature(::Euclidean, p, X, Y)
+
+Sectional curvature of [`Euclidean`](@ref) manifold `M` is 0.
+"""
+function sectional_curvature(::Euclidean, p, X, Y)
+    return 0.0
+end
+
+@doc raw"""
+    sectional_curvature_max(::Euclidean)
+
+Sectional curvature of [`Euclidean`](@ref) manifold `M` is 0.
+"""
+function sectional_curvature_max(::Euclidean)
+    return 0.0
+end
+
+@doc raw"""
+    sectional_curvature_min(M::Euclidean)
+
+Sectional curvature of [`Euclidean`](@ref) manifold `M` is 0.
+"""
+function sectional_curvature_min(::Euclidean)
+    return 0.0
+end
+
 function Base.show(io::IO, M::Euclidean{N,𝔽}) where {N<:Tuple,𝔽}
     size = get_parameter(M.size)
     return print(io, "Euclidean($(join(size, ", ")); field=$(𝔽), parameter=:field)")

src/manifolds/GeneralUnitaryMatrices.jl

@@ -182,6 +182,8 @@ function default_approximation_method(::GeneralUnitaryMatrices{<:Any,ℝ}, ::typ
     return GeodesicInterpolationWithinRadius(π / 2 / √2)
 end
 
+default_vector_transport_method(::GeneralUnitaryMatrices) = ProjectionTransport()
+
 embed(::GeneralUnitaryMatrices, p) = p
 
 @doc raw"""

src/manifolds/Hyperbolic.jl

@@ -415,3 +415,42 @@ function riemann_tensor!(M::Hyperbolic, W, p, X, Y, Z)
     W .= inner(M, p, Z, X) .* Y .- inner(M, p, Z, Y) .* X
     return W
 end
+
+@doc raw"""
+    sectional_curvature(::Hyperbolic, p, X, Y)
+
+Sectional curvature of [`Hyperbolic`](@ref) `M` is -1 if dimension is > 1 and 0 otherwise.
+"""
+function sectional_curvature(M::Hyperbolic, p, X, Y)
+    if manifold_dimension(M) > 1
+        return -1.0
+    else
+        return 0.0
+    end
+end
+
+@doc raw"""
+    sectional_curvature_max(::Hyperbolic)
+
+Sectional curvature of [`Hyperbolic`](@ref) `M` is -1 if dimension is > 1 and 0 otherwise.
+"""
+function sectional_curvature_max(M::Hyperbolic)
+    if manifold_dimension(M) > 1
+        return -1.0
+    else
+        return 0.0
+    end
+end
+
+@doc raw"""
+    sectional_curvature_min(M::Hyperbolic)
+
+Sectional curvature of [`Hyperbolic`](@ref) `M` is -1 if dimension is > 1 and 0 otherwise.
+"""
+function sectional_curvature_min(M::Hyperbolic)
+    if manifold_dimension(M) > 1
+        return -1.0
+    else
+        return 0.0
+    end
+end

src/manifolds/Rotations.jl

@@ -62,6 +62,8 @@ function angles_4d_skew_sym_matrix(A)
     return sqrt(halfb + sqrtdisc), sqrt(halfb - sqrtdisc)
 end
 
+default_vector_transport_method(::Rotations) = ParallelTransport()
+
 # from https://github.com/JuliaManifolds/Manifolds.jl/issues/453#issuecomment-1046057557
 function _get_tridiagonal_elements(trian)
     N = size(trian, 1)
@@ -423,6 +425,35 @@ function riemannian_Hessian!(M::Rotations, Y, p, G, H, X)
     return Y
 end
 
+@doc raw"""
+    sectional_curvature_max(::Rotations)
+
+Sectional curvature of [`Rotations`](@ref) `M` is equal to 0 for `Rotations(1)` and
+`Rotations(2)`, less than or equal to 1/8 for `Rotations(3)` and less than or equal to 1/4
+for higher-dimensional rotations manifolds.
+
+For reference, see [Ge:2014](@cite), Lemma 2.5 and [CheegerEbin:2008](@cite), Corollary 3.19.
+"""
+function sectional_curvature_max(M::Rotations)
+    N = manifold_dimension(M)
+    if N <= 1
+        return 0.0
+    elseif N == 3
+        return 1 / 8
+    else
+        return 1 / 4
+    end
+end
+
+@doc raw"""
+    sectional_curvature_min(M::Rotations)
+
+Sectional curvature of [`Rotations`](@ref) `M` is greater than or equal to 0.
+"""
+function sectional_curvature_min(::Rotations)
+    return 0.0
+end
+
 Distributions.support(d::NormalRotationDistribution) = MPointSupport(d.manifold)
 
 @doc raw"""

src/manifolds/Sphere.jl

@@ -571,6 +571,48 @@ function riemann_tensor!(M::AbstractSphere{ℝ}, Xresult, p, X, Y, Z)
     return Xresult
 end
 
+@doc raw"""
+    sectional_curvature(::AbstractSphere, p, X, Y)
+
+Sectional curvature of [`AbstractSphere`](@ref) `M` is 1 if dimension is greater than 1
+and 0 otherwise.
+"""
+function sectional_curvature(M::AbstractSphere, p, X, Y)
+    if manifold_dimension(M) > 1
+        return 1.0
+    else
+        return 0.0
+    end
+end
+
+@doc raw"""
+    sectional_curvature_max(::AbstractSphere)
+
+Sectional curvature of [`AbstractSphere`](@ref) `M` is 1 if dimension is greater than 1
+and 0 otherwise.
+"""
+function sectional_curvature_max(M::AbstractSphere)
+    if manifold_dimension(M) > 1
+        return 1.0
+    else
+        return 0.0
+    end
+end
+
+@doc raw"""
+    sectional_curvature_min(M::AbstractSphere)
+
+Sectional curvature of [`AbstractSphere`](@ref) `M` is 1 if dimension is greater than 1
+and 0 otherwise.
+"""
+function sectional_curvature_min(M::AbstractSphere)
+    if manifold_dimension(M) > 1
+        return 1.0
+    else
+        return 0.0
+    end
+end
+
 @doc raw"""
     volume_density(M::AbstractSphere{ℝ}, p, X)
 

src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl

@@ -454,6 +454,30 @@ function riemann_tensor!(::SymmetricPositiveDefinite, Xresult, p, X, Y, Z)
     return Xresult
 end
 
+"""
+    sectional_curvature_min(M::SymmetricPositiveDefinite)
+
+Return minimum sectional curvature of [`SymmetricPositiveDefinite`](@ref) manifold,
+that is 0 for SPD(1) and SPD(2) and -0.25 otherwise.
+"""
+function sectional_curvature_min(M::SymmetricPositiveDefinite)
+    if manifold_dimension(M) < 2
+        return 0.0
+    else
+        return -0.25
+    end
+end
+
+"""
+    sectional_curvature_max(M::SymmetricPositiveDefinite)
+
+Return minimum sectional curvature of [`SymmetricPositiveDefinite`](@ref) manifold,
+that is 0.
+"""
+function sectional_curvature_max(::SymmetricPositiveDefinite)
+    return 0.0
+end
+
 """
     volume_density(::SymmetricPositiveDefinite, p, X)