diff --git a/.github/workflows/documenter.yml b/.github/workflows/documenter.yml index 3536a6c71f..3c9d5d4372 100644 --- a/.github/workflows/documenter.yml +++ b/.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" diff --git a/NEWS.md b/NEWS.md index 5236d33b90..faf84201c4 100644 --- a/NEWS.md +++ b/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 diff --git a/Project.toml b/Project.toml index a7a7c8c0ff..f9d8718cd4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Manifolds" uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann ", "Antoine Levitt "] -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" diff --git a/docs/make.jl b/docs/make.jl index d29adc83d2..838774f93c 100755 --- a/docs/make.jl +++ b/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" => [ diff --git a/docs/src/features/utilities.md b/docs/src/features/utilities.md index afc4c9f34e..e766f30ea2 100644 --- a/docs/src/features/utilities.md +++ b/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. diff --git a/docs/src/misc/notation.md b/docs/src/misc/notation.md index 23a55d0f24..0b6033a884 100644 --- a/docs/src/misc/notation.md +++ b/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`` | | | diff --git a/docs/src/references.bib b/docs/src/references.bib index 3f2cc96041..e377fd0737 100644 --- a/docs/src/references.bib +++ b/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}}, diff --git a/src/Manifolds.jl b/src/Manifolds.jl index 78cd5438fa..317c663cef 100644 --- a/src/Manifolds.jl +++ b/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!, diff --git a/src/manifolds/Euclidean.jl b/src/manifolds/Euclidean.jl index 1a4ce3f2ba..0098faa3a5 100644 --- a/src/manifolds/Euclidean.jl +++ b/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)") diff --git a/src/manifolds/GeneralUnitaryMatrices.jl b/src/manifolds/GeneralUnitaryMatrices.jl index 667daeea49..4bd682dd96 100644 --- a/src/manifolds/GeneralUnitaryMatrices.jl +++ b/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""" diff --git a/src/manifolds/Hyperbolic.jl b/src/manifolds/Hyperbolic.jl index 3f734df4cc..f7be05f38b 100644 --- a/src/manifolds/Hyperbolic.jl +++ b/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 diff --git a/src/manifolds/Rotations.jl b/src/manifolds/Rotations.jl index 849c880f6f..72e27c0649 100644 --- a/src/manifolds/Rotations.jl +++ b/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""" diff --git a/src/manifolds/Sphere.jl b/src/manifolds/Sphere.jl index 920d058df6..d5266d0d52 100644 --- a/src/manifolds/Sphere.jl +++ b/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) diff --git a/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl b/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl index 8fa47f017b..f223c94f68 100644 --- a/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl +++ b/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) diff --git a/src/utils.jl b/src/utils.jl index 7ce85504a3..77bc1e01de 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -277,3 +277,90 @@ _eps_safe(::Type{T}) where {T<:Real} = eps(T) _eps_safe(::Type{T}) where {T<:Number} = eps(real(T)) max_eps(xs...) = maximum(_eps_safe ∘ eltype, xs) + +""" + sectional_curvature_matrix(M::AbstractManifold, p, B::AbstractBasis) + +Compute the matrix of sectional curvatures of manifold `M` at point `p`. +Entry `(i, j)` corresponds to sectional curvature of the surface spanned by vectors +`i` and `j` from basis `B`. +""" +function sectional_curvature_matrix(M::AbstractManifold, p, B::AbstractBasis) + V = get_vectors(M, p, get_basis(M, p, B)) + N = length(V) + result = zeros(N, N) + for (i, e_i) in enumerate(V) + for (j, e_j) in enumerate(V) + if i < j + result[i, j] = sectional_curvature(M, p, e_i, e_j) + result[j, i] = result[i, j] + end + end + end + return result +end + +@doc raw""" + estimated_sectional_curvature(M::AbstractManifold, p, X, Y; r::Real=1e-3, N::Int=10000) + +Approximate sectional curvature of manifold `M` in the plane spanned by vectors `X` and `Y` +from tangent space at `p` using a circle on `M` of radius `r` divided into `N` segments. + +The approximation is derived from the [Bertrand–Diguet–Puiseux theorem](https://en.wikipedia.org/wiki/Bertrand%E2%80%93Diguet%E2%80%93Puiseux_theorem) +which states that +````math +\kappa_p(X, Y) = \lim_{r \to 0^+} 3\frac{2\pi r-C(r)}{\pi r^3}, +```` +where ``C(r)`` is the circumference of the circle of radius ``r`` around `p` in submanifold +of `M` spanned by `X` and `Y`. The circumference calculation method has a tendency to +return curvature values larger than the exact ones. +""" +function estimated_sectional_curvature( + M::AbstractManifold, + p, + X, + Y; + r::Real=1e-3, + N::Int=10000, +) + circumference = 0.0 + p_i = similar(p) + p_ip1 = similar(p) + for i in 1:N + θ_i = 2π * (i - 1) / N + θ_ip1 = 2π * (i) / N + exp!(M, p_i, p, r .* (sin(θ_i) .* X .+ cos(θ_i) .* Y)) + exp!(M, p_ip1, p, r .* (sin(θ_ip1) .* X .+ cos(θ_ip1) .* Y)) + + circumference += distance(M, p_i, p_ip1) + end + return 3 * (2π * r - circumference) / (π * r^3) +end + +""" + estimated_sectional_curvature_matrix(M::AbstractManifold, p, B::AbstractBasis; r::Real=1e-3, N::Int=10000) + +Estimate the matrix of sectional curvatures of manifold `M` at point `p` using +`estimated_sectional_curvature`. Entry `(i, j)`` corresponds to sectional curvature of the +surface spanned by vectors `i` and `j` from basis `B`. +""" +function estimated_sectional_curvature_matrix( + M::AbstractManifold, + p, + B::AbstractBasis; + r::Real=1e-3, + N_pts::Int=10000, +) + V = get_vectors(M, p, get_basis(M, p, B)) + N = length(V) + result = zeros(N, N) + for (i, e_i) in enumerate(V) + for (j, e_j) in enumerate(V) + if i < j + result[i, j] = estimated_sectional_curvature(M, p, e_i, e_j; r=r, N=N_pts) + result[j, i] = result[i, j] + end + end + end + return result +end diff --git a/test/groups/power_group.jl b/test/groups/power_group.jl index d2aa761b94..ede5304e18 100644 --- a/test/groups/power_group.jl +++ b/test/groups/power_group.jl @@ -30,7 +30,7 @@ include("group_utils.jl") pts, X_pts, X_pts; - atol=2e-9, + atol=2e-8, test_diff=true, test_log_from_identity=true, test_exp_from_identity=true, diff --git a/test/manifolds/euclidean.jl b/test/manifolds/euclidean.jl index 64e856177a..0099f030ec 100644 --- a/test/manifolds/euclidean.jl +++ b/test/manifolds/euclidean.jl @@ -296,8 +296,12 @@ using FiniteDifferences @test size(RT) == (2, 2, 2, 2) @test norm(RT) ≈ 0.0 atol = 1e-16 - @test riemann_tensor(M, p, [1, 2], [1, 3], [1, 4]) == [0, 0] @test !Manifolds.check_chart_switch(M, A, i, p) + + @test riemann_tensor(M, p, [1, 2], [1, 3], [1, 4]) == [0, 0] + @test sectional_curvature(M, p, [1.0, 0.0], [0.0, 1.0]) == 0.0 + @test sectional_curvature_max(M) == 0.0 + @test sectional_curvature_min(M) == 0.0 end @testset "Induced Basis and local metric for EuclideanMetric" begin struct DefaultManifold <: AbstractManifold{ℝ} end diff --git a/test/manifolds/hyperbolic.jl b/test/manifolds/hyperbolic.jl index f4970ca3d1..07e18b4d67 100644 --- a/test/manifolds/hyperbolic.jl +++ b/test/manifolds/hyperbolic.jl @@ -347,23 +347,22 @@ include("../header.jl") end @testset "other metric" begin M = Hyperbolic(2) - @test riemann_tensor( - M, - [0.0, 0.0, -1.0], - [-1.0, 0.0, 0.0], - [0.0, -1.0, 0.0], - [-0.5, -0.5, 0.0], - ) == [0.5, -0.5, 0.0] - X = [0.0, 0.0, 0.0] - @test riemann_tensor!( - M, - X, - [0.0, 0.0, -1.0], - [-1.0, 0.0, 0.0], - [0.0, -1.0, 0.0], - [-0.5, -0.5, 0.0], - ) === X - @test X == [0.5, -0.5, 0.0] + p = [0.0, 0.0, -1.0] + X = [-1.0, 0.0, 0.0] + Y = [0.0, -1.0, 0.0] + Z = [-0.5, -0.5, 0.0] + @test riemann_tensor(M, p, X, Y, Z) == [0.5, -0.5, 0.0] + Xout = [0.0, 0.0, 0.0] + @test riemann_tensor!(M, Xout, p, X, Y, Z) === Xout + @test Xout == [0.5, -0.5, 0.0] + + @test sectional_curvature(M, p, X, Y) == -1.0 + @test sectional_curvature_max(M) == -1.0 + @test sectional_curvature_min(M) == -1.0 + M1 = Hyperbolic(1) + @test sectional_curvature(M1, p, X, Y) == 0.0 + @test sectional_curvature_max(M1) == 0.0 + @test sectional_curvature_min(M1) == 0.0 end @testset "ManifoldDiff" begin # ManifoldDiff diff --git a/test/manifolds/multinomial_spd.jl b/test/manifolds/multinomial_spd.jl index 70a6027a22..ea37d06023 100644 --- a/test/manifolds/multinomial_spd.jl +++ b/test/manifolds/multinomial_spd.jl @@ -34,6 +34,6 @@ include("../header.jl") @testset "Random" begin q = zeros(3, 3) M = MultinomialSymmetricPositiveDefinite(3) - @test is_point(M, rand!(MersenneTwister(), M, q)) + @test is_point(M, rand!(MersenneTwister(), M, q); atol=1e-15) end end diff --git a/test/manifolds/rotations.jl b/test/manifolds/rotations.jl index dfd4aee46e..a36b646597 100644 --- a/test/manifolds/rotations.jl +++ b/test/manifolds/rotations.jl @@ -17,6 +17,7 @@ include("../header.jl") TEST_FLOAT32 && push!(types, Matrix{Float32}) TEST_STATIC_SIZED && push!(types, MMatrix{2,2,Float64,4}) retraction_methods = [PolarRetraction(), QRRetraction()] + @test default_vector_transport_method(M) === ParallelTransport() inverse_retraction_methods = [PolarInverseRetraction(), QRInverseRetraction()] @@ -275,6 +276,14 @@ include("../header.jl") ] end + @testset "sectional curvature" begin + @test sectional_curvature_min(Rotations(3)) == 0.0 + @test sectional_curvature_max(Rotations(1)) == 0.0 + @test sectional_curvature_max(Rotations(2)) == 0.0 + @test sectional_curvature_max(Rotations(3)) == 1 / 8 + @test sectional_curvature_max(Rotations(4)) == 1 / 4 + end + @testset "field parameter" begin M = Rotations(2; parameter=:field) @test is_flat(M) diff --git a/test/manifolds/sphere.jl b/test/manifolds/sphere.jl index 9031742fb2..e81c340079 100644 --- a/test/manifolds/sphere.jl +++ b/test/manifolds/sphere.jl @@ -255,7 +255,18 @@ using ManifoldsBase: TFVector @testset "other metric" begin M = Sphere(2) - @test riemann_tensor(M, [0, 1, 0], [1, 0, 0], [1, 0, 4], [0, 0, 1]) == [4, 0, 0] + p = [0, 1, 0] + X = [1, 0, 0] + Y = [1, 0, 4] + Z = [0, 0, 1] + @test riemann_tensor(M, p, X, Y, Z) == [4, 0, 0] + @test sectional_curvature(M, p, X, Y) == 1.0 + @test sectional_curvature_max(M) == 1.0 + @test sectional_curvature_min(M) == 1.0 + M1 = Sphere(1) + @test sectional_curvature(M1, p, X, Y) == 0.0 + @test sectional_curvature_max(M1) == 0.0 + @test sectional_curvature_min(M1) == 0.0 end @testset "ManifoldDiff" begin diff --git a/test/manifolds/unitary_matrices.jl b/test/manifolds/unitary_matrices.jl index b0fc87cb06..ccd5a889f3 100644 --- a/test/manifolds/unitary_matrices.jl +++ b/test/manifolds/unitary_matrices.jl @@ -19,6 +19,7 @@ using Quaternions @test is_point(M, rand(MersenneTwister(), M); error=:error) @test abs(rand(MersenneTwister(), OrthogonalMatrices(1))[]) == 1 @test is_vector(M, p, rand(MersenneTwister(), M; vector_at=p)) + @test default_vector_transport_method(M) === ProjectionTransport() end @testset "Unitary Matrices" begin diff --git a/tutorials/_quarto.yml b/tutorials/_quarto.yml index c98d981b6e..5befbaf19d 100644 --- a/tutorials/_quarto.yml +++ b/tutorials/_quarto.yml @@ -24,4 +24,4 @@ format: variant: -raw_html+tex_math_dollars wrap: preserve -jupyter: julia-1.10 \ No newline at end of file +jupyter: julia-1.10 diff --git a/tutorials/exploring-curvature.qmd b/tutorials/exploring-curvature.qmd new file mode 100644 index 0000000000..8d13e4f03c --- /dev/null +++ b/tutorials/exploring-curvature.qmd @@ -0,0 +1,95 @@ +--- +title: Exploring curvature without coordinates +--- + +```{julia} +#| echo: false +#| code-fold: true +#| output: false +using Pkg; +cd(@__DIR__) +Pkg.activate("."); # for reproducibility use the local tutorial environment. +Pkg.develop(PackageSpec(; path=(@__DIR__) * "/../")) +using Markdown +``` + +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`](@ref) function. +Note that Riemann tensor and sectional curvature are equivalently full specifications of curvature in a manifold, see [CheegerEbin:2008](@cite), Eq. (1.12). + +Let's take the [`SymmetricPositiveDefinite`](@ref) manifold as our first example. +It has nonpositive sectional curvature: + +```{julia} +using Manifolds +using LinearAlgebra +M = SymmetricPositiveDefinite(3) +p = rand(M) +cm = sectional_curvature_matrix(M, p, DefaultOrthonormalBasis()) +``` + +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: +```{julia} +cm_bdp = Manifolds.estimated_sectional_curvature_matrix(M, p, DefaultOrthonormalBasis(); r=1e-3, N_pts=100000) +println(norm(cm - cm_bdp)) +``` + +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. +```{julia} +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 +``` + +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. + +```{julia} +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) +``` + +In the planes corresponding to orthonormal basis, the maximum sectional curvature is 0.125 but the true upper bound is 0.25. + +## Literature + +```@bibliography +Pages = ["exploring-curvature.md"] +Canonical=false +``` \ No newline at end of file