From f03ac74bcbdf81cd80d63cf58f53d2c27982f726 Mon Sep 17 00:00:00 2001 From: Simon Jacobsson <47160031+sjacobsson@users.noreply.github.com> Date: Fri, 20 Dec 2024 20:27:58 +0100 Subject: [PATCH] add segre and warped segre manifold (#755) * add Segre and warped Segre manifold * added formulas to Segre manifold documentation * Setup docs, fix a few type errors, runs formatter. * Add link to Segre. --------- Co-authored-by: Simon Jacobsson Co-authored-by: Ronny Bergmann --- .github/workflows/documenter.yml | 2 +- .zenodo.json | 6 + NEWS.md | 4 +- _quarto.yml | 26 - docs/make.jl | 1 + docs/src/manifolds/segre.md | 31 ++ docs/src/misc/about.md | 1 + docs/src/references.bib | 12 +- .../config/vocabularies/Manifolds/accept.txt | 1 + src/Manifolds.jl | 7 +- src/manifolds/Segre.jl | 496 ++++++++++++++++++ src/manifolds/SegreWarpedMetric.jl | 317 +++++++++++ test/manifolds/segre.jl | 349 ++++++++++++ test/runtests.jl | 1 + tutorials/_quarto.yml | 2 +- 15 files changed, 1225 insertions(+), 31 deletions(-) delete mode 100644 _quarto.yml create mode 100644 docs/src/manifolds/segre.md create mode 100644 src/manifolds/Segre.jl create mode 100644 src/manifolds/SegreWarpedMetric.jl create mode 100644 test/manifolds/segre.jl diff --git a/.github/workflows/documenter.yml b/.github/workflows/documenter.yml index f16edf2929..73427534fe 100644 --- a/.github/workflows/documenter.yml +++ b/.github/workflows/documenter.yml @@ -17,7 +17,7 @@ jobs: version: "1.6.38" - uses: julia-actions/setup-julia@latest with: - version: "1.10" + version: "1.11" - name: Julia Cache uses: julia-actions/cache@v2 - name: Cache Quarto diff --git a/.zenodo.json b/.zenodo.json index 4076fbb329..04322d2cce 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -31,6 +31,12 @@ "name": "Weiß, Manuel", "type": "ProjectMember" }, + { + "affiliation": "KU Leuven", + "name": "Jacobsson, Simon", + "orcid": "0000-0002-1181-972X", + "type": "ProjectMember" + }, { "affiliation": "Georg-August-University Göttingen", "name": "Klingbiel, Lukas", diff --git a/NEWS.md b/NEWS.md index 97c697b15c..eb55cd16fd 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,10 +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.10.9] - unreleased +## [0.10.9] - 2024-12-20 ### Added +* the `Segre` manifold +* the `WarpedMetric` for the `Segre`manifold * The manifold `HeisenbergMatrices` as the underlying manifold of `HeisenbergGroup`. ### Changed diff --git a/_quarto.yml b/_quarto.yml deleted file mode 100644 index 44421d48c1..0000000000 --- a/_quarto.yml +++ /dev/null @@ -1,26 +0,0 @@ - -project: - title: "Manifolds.jl" - #execute-dir: project - -crossref: - fig-prefix: Figure - tbl-prefix: Table -bibliography: - - CITATION.bib -# - bib/manopt-bibliography.yaml -# csl: journal-of-the-royal-statistical-society.csl -fig-format: png - -execute: - freeze: auto - eval: true - echo: true - output: true - -format: - commonmark: - variant: -raw_html+tex_math_dollars - wrap: preserve - -jupyter: julia-1.9 diff --git a/docs/make.jl b/docs/make.jl index 3150bfad52..58d7d3319d 100755 --- a/docs/make.jl +++ b/docs/make.jl @@ -201,6 +201,7 @@ makedocs(; "Projective space" => "manifolds/projectivespace.md", "Orthogonal and Unitary Matrices" => "manifolds/generalunitary.md", "Rotations" => "manifolds/rotations.md", + "Segre" => "manifolds/segre.md", "Shape spaces" => "manifolds/shapespace.md", "Skew-Hermitian matrices" => "manifolds/skewhermitian.md", "Spectrahedron" => "manifolds/spectrahedron.md", diff --git a/docs/src/manifolds/segre.md b/docs/src/manifolds/segre.md new file mode 100644 index 0000000000..a60c11c2f2 --- /dev/null +++ b/docs/src/manifolds/segre.md @@ -0,0 +1,31 @@ +# The Segre manifold + +```@docs +Segre +``` + +```@autodocs +Modules = [Manifolds] +Pages = ["manifolds/Segre.jl"] +Order = [:function] +``` + +## [A warped metric](@id segre-warped-metric-sec) + +```@docs +WarpedMetric +``` + +```@autodocs +Modules = [Manifolds] +Pages = ["manifolds/SegreWarpedMetric.jl"] +Order = [:function] +``` + + +## Literature + +```@bibliography +Pages = ["segre.md"] +Canonical=false +``` diff --git a/docs/src/misc/about.md b/docs/src/misc/about.md index f70a1b9adc..613f0aa687 100644 --- a/docs/src/misc/about.md +++ b/docs/src/misc/about.md @@ -17,6 +17,7 @@ - [Nick Dewaele](https://github.com/Nikdwal) contributed the [Tucker manifold](../manifolds/tucker.md) - [Renée Dornig](https://github.com/r-dornig) contributed the [centered matrices](../manifolds/centeredmatrices.md) and the [essential manifold](../manifolds/essentialmanifold.md) - [David Hong](https://github.com/dahong67) contributed uniform distributions on the Stiefel and Grassmann manifolds. +- [Simon Jacobsson](https://github.com/sjacobsson) contributed the [Segre manifold](../manifolds/segre.md) including its [warped metric](@ref segre-warped-metric-sec) thereon. - [Johannes Voll Kolstø](https://github.com/johannvk) contributed the [symplectic manifold](../manifolds/symplectic.md), the [symplectic Stiefel manifold](../manifolds/symplecticstiefel.md) - [Manuel Weiß](https://github.com/manuelweisser) contributed [symmetric matrices](../manifolds/symmetric.md) diff --git a/docs/src/references.bib b/docs/src/references.bib index d01b1c751c..b7b8c2f890 100644 --- a/docs/src/references.bib +++ b/docs/src/references.bib @@ -478,6 +478,16 @@ @article{HueperMarkinaSilvaLeite:2021 # # J # ---------------------------------------------------------------------------------------- +@misc{JacobssonSwijsenVandervekenVannieuwenhoven:2024, + title={Warped geometries of Segre-Veronese manifolds}, + author={Simon Jacobsson and Lars Swijsen and Joeri Van der Veken and Nick Vannieuwenhoven}, + year={2024}, + eprint={2410.00664}, + archivePrefix={arXiv}, + primaryClass={math.NA}, + url={https://arxiv.org/abs/2410.00664}, +} + @article{JourneeBachAbsilSepulchre:2010, DOI = {10.1137/080731359}, EPRINT = {0807.4423}, @@ -951,4 +961,4 @@ @article{AastroemPetraSchmitzerSchnoerr:2017 AUTHOR = {Freddie Åström and Stefania Petra and Bernhard Schmitzer and Christoph Schnörr}, TITLE = {Image Labeling by Assignment}, JOURNAL = {Journal of Mathematical Imaging and Vision} -} \ No newline at end of file +} diff --git a/docs/styles/config/vocabularies/Manifolds/accept.txt b/docs/styles/config/vocabularies/Manifolds/accept.txt index bb1628d587..9285310b82 100644 --- a/docs/styles/config/vocabularies/Manifolds/accept.txt +++ b/docs/styles/config/vocabularies/Manifolds/accept.txt @@ -1,6 +1,7 @@ Grassmann [Jj]ulia Riemannian +Segre Stiefel [sS]ymplectic struct \ No newline at end of file diff --git a/src/Manifolds.jl b/src/Manifolds.jl index b74815e7a8..701c3631aa 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -453,6 +453,8 @@ include("manifolds/MultinomialSymmetric.jl") include("manifolds/MultinomialSymmetricPositiveDefinite.jl") include("manifolds/PositiveNumbers.jl") include("manifolds/ProjectiveSpace.jl") +include("manifolds/Segre.jl") +include("manifolds/SegreWarpedMetric.jl") include("manifolds/SkewHermitian.jl") include("manifolds/Spectrahedron.jl") include("manifolds/Stiefel.jl") @@ -660,6 +662,7 @@ export Euclidean, ProbabilitySimplex, ProjectiveSpace, Rotations, + Segre, SkewHermitianMatrices, SkewSymmetricMatrices, Spectrahedron, @@ -738,7 +741,8 @@ export AbstractMetric, ProductMetric, RealSymplecticMetric, RiemannianMetric, - StiefelSubmersionMetric + StiefelSubmersionMetric, + WarpedMetric export AbstractAtlas, RetractionAtlas # Vector transport types export AbstractVectorTransportMethod, ParallelTransport, ProjectionTransport @@ -807,6 +811,7 @@ export ×, christoffel_symbols_first, christoffel_symbols_second, christoffel_symbols_second_jacobian, + connected_by_geodesic, convert, complex_dot, decorated_manifold, diff --git a/src/manifolds/Segre.jl b/src/manifolds/Segre.jl new file mode 100644 index 0000000000..0d6cfa1f7a --- /dev/null +++ b/src/manifolds/Segre.jl @@ -0,0 +1,496 @@ +@doc raw""" + Segre{𝔽,V} <: AbstractManifold{𝔽} + +The Segre manifold + +````math + \mathcal{S} = \operatorname{Seg}(𝔽^{n_1} \times \dots \times 𝔽^{n_d}) +```` + +is the set of rank-one tensors in ``𝔽^{n_1} \otimes \dots \otimes 𝔽^{n_d}``. + +When ``𝔽 = ℝ``, the Segre manifold is a normal Riemannian covering of + +````math + \mathcal{P} = ℝ^{+} \times \mathbb{S}^{n_1 - 1} \times \dots \times \mathbb{S}^{n_d - 1} +```` + +equipped with a [warped product metric](https://en.wikipedia.org/wiki/Warped_product). The tuple ``(n_1, \dots, n_d)`` is called the _valence_ of the manifold. + +The geometry of the Segre manifold is summarized in [JacobssonSwijsenVandervekenVannieuwenhoven:2024](@cite). It is named after [Corrado Segre](https://en.wikipedia.org/wiki/Corrado_Segre)(1863–1924). + +# Constructor + Segre(n::Int...; field::AbstractNumbers=ℝ) + +Generate a valence `(n, ...)` Segre manifold. +`Segre(n)` is the same as ``\mathbb{R}^{n} \setminus \{ 0 \}``. +""" +struct Segre{𝔽,V} <: AbstractManifold{𝔽} end + +function Segre(n::Int...; field::AbstractNumbers=ℝ) + return Segre{field,(n...,)}() +end + +function check_size(M::Segre{𝔽,V}, p) where {𝔽,V} + p_size = only.(size.(p)) + M_size = [1, V...] + + if p_size != M_size + return DomainError( + p_size, + "The point $(p) can not belong to the manifold $(M), since its size $(p_size) is not equal to the manifolds representation size ($(M_size)).", + ) + end + + return nothing +end + +function check_size(M::Segre{𝔽,V}, p, X;) where {𝔽,V} + p_size = only.(size.(p)) + X_size = only.(size.(X)) + M_size = [1, V...] + + if p_size != M_size + return DomainError( + p_size, + "The point $(p) can not belong to the manifold $(M), since its size $(p_size) is not equal to the manifolds representation size ($(M_size)).", + ) + end + + if X_size != M_size + return DomainError( + X_size, + "The vector $(X) can not belong to the manifold $(M), since its size $(X_size) is not equal to the manifolds representation size ($(M_size)).", + ) + end + + return nothing +end + +@doc raw""" + closest_representative!(M::Segre{ℝ, V}, p, q) + +``\mathcal{S}`` is a ``2^d``-sheeted Riemannian covering of +````math + \mathcal{P} = ℝ^{+} \times \mathbb{S}^{n_1 - 1} \times \dots \times \mathbb{S}^{n_d - 1} +```` +with a warped product metric. +Every equivalence class ``q \in \mathcal{S}`` has ``2^d`` representatives in ``\mathcal{P}``. +`closest_representative!(M, q, p)` changes representative of `q` to the one that is closest to `p` in ``\mathcal{P}``. +""" +function closest_representative!(M::Segre{ℝ,V}, q, p) where {V} + + # Find closest representation by flipping an even number of signs. + ds = [distance(Sphere(n - 1), x, y) for (n, x, y) in zip(V, p[2:end], q[2:end])] + flips = [false, (ds .> (pi / 2))...] + nbr_flips = sum(flips) + + # This code is pretty ugly. + if isodd(nbr_flips) + if nbr_flips == length(V) + flips[argmin(ds) + 1] = false + else + is = sortperm(ds; rev=true) + + flips1 = deepcopy(flips) + flips1[is[nbr_flips] + 1] = false + q1 = deepcopy(q) + q1[flips1] = -q1[flips1] + + flips2 = deepcopy(flips) + flips2[is[nbr_flips + 1] + 1] = true + q2 = deepcopy(q) + q2[flips2] = -q2[flips2] + + spherical_angle_sum(M, p, q1) < spherical_angle_sum(M, p, q2) ? flips = flips1 : + flips = flips2 + end + end + + return q[flips] = -q[flips] +end + +@doc raw""" + connected_by_geodesic(M::Segre{ℝ, V}, p, q) + +``\mathcal{S}`` is not a complete manifold, i.e. not every pair `p` and `q` of points are connected by a geodesic in ``\mathcal{S}``. +`connected_by_geodesic(M, p, q)` returns `true` if two points, `p` and `q`, are connected by a geodesic, and otherwise returns `false`. +""" +function connected_by_geodesic(M::Segre{ℝ,V}, p, q) where {V} + closest_representative!(M, q, p) + return spherical_angle_sum(M, p, q) < pi +end + +@doc raw""" + embed(M::Segre{𝔽, V}, p) + embed!(M::Segre{𝔽, V}, q, p) + +Embed ``p ≐ (λ, x_1, …, x_d)`` in ``𝔽^{n_1 ×⋯× n_d}`` using the Kronecker product + +````math + (λ, x_1, …, x_d) ↦ λ x_1 ⊗⋯⊗ x_d. +```` +""" +embed(::Segre, p) + +function embed!(M::Segre, q, p) + return q = kron(p...) +end + +@doc raw""" + embed!(M::Segre{𝔽, V}, p, X) + +Embed tangent vector ``X = (ν, u_1, …, u_d)`` at ``p ≐ (λ, x_1, …, x_d)`` in ``𝔽^{n_1 ×⋯× n_d}`` using the Kronecker product + +````math + (ν, u_1, …, u_d) ↦ ν x_1 ⊗⋯⊗ x_d + λ u_1 ⊗ x_2 ⊗⋯⊗ x_d + … + λ x_1 ⊗⋯⊗ x_{d - 1} ⊗ u_d. +```` +""" +function embed!(::Segre{𝔽,V}, u, p, X) where {𝔽,V} + # Product rule + return u = sum([ + kron([i == j ? xdot : x for (j, (x, xdot)) in enumerate(zip(p, X))]...) for + (i, _) in enumerate(p) + ]) +end + +""" + is_point(M::Segre{ℝ, V}, p; kwargs...) + +Check whether `p` is a valid point on `M`, i.e. `p[1]` is a singleton containing a positive number and `p[i + 1]` is a point on `Sphere(V[i])`. +The tolerance can be set using the `kwargs...`. +""" +is_point(M::Segre{ℝ,V}, p; kwargs...) where {V} + +function check_point(M::Segre{ℝ,V}, p; atol=1.4901161193847656e-8, kwargs...) where {V} + if p[1][1] <= 0.0 + return DomainError(p[1][1], "$(p) has non-positive modulus.") + end + + for (x, n) in zip(p[2:end], V) + e = check_point(Sphere(n - 1)::AbstractSphere, x; atol=atol, kwargs...) + if !isnothing(e) + return e + end + end +end + +""" + is_vector(M::Segre{ℝ, V}, p, X, kwargs...) + +Check whether `X` is a tangent vector to `p` on `M`, i.e. `X` has to be of same dimension as `p` and orthogonal to `p`. +The tolerance can be set using the `kwargs...`. +""" +is_vector(M::Segre{ℝ,V}, p, v; kwargs...) where {V} + +function check_vector(M::Segre{ℝ,V}, p, X; atol=1.4901161193847656e-8, kwargs...) where {V} + for (x, xdot, n) in zip(p[2:end], X[2:end], V) + e = check_vector(Sphere(n - 1)::AbstractSphere, x, xdot; atol=atol, kwargs...) + if !isnothing(e) + return e + end + end +end + +@doc raw""" + distance(M::Segre{ℝ, V}, p, q) + +Riemannian distance between two points `p` and `q` on the Segre manifold. + +Assume ``p ≐ (λ, x_1, …, x_d)``, ``q ≐ (μ, y_1, …, y_d) ∈ \mathcal{S}`` are connected by a geodesic. Let + +````math + m = \sqrt{\sphericalangle(x_1, y_1)^2 + … + \sphericalangle(x_d, y_d)^2} +```` + +and assume ``(μ, y_1, …, y_d)`` is the representation of ``q`` that minimizes ``m``. Then + +````math + \operatorname{dist}_{\mathcal{S}}(p, q) = \sqrt{λ^2 - 2 λμ\cos(m) + μ^2}. +```` +""" +function distance(M::Segre{ℝ,V}, p, q) where {V} + closest_representative!(M, q, p) + m = spherical_angle_sum(M, p, q) + return sqrt((p[1][1] - q[1][1])^2 + 4 * p[1][1] * q[1][1] * sin(m / 2)^2) + # Equivalent to sqrt(p[1][1]^2 + q[1][1]^2 - 2 * p[1][1] * q[1][1] * cos(m)) but more stable for small m +end + +@doc raw""" + exp(M::Segre{ℝ, V}, p, X) + +Exponential map on the Segre manifold. + +Let ``p ≐ (λ, x_1, …, x_d) ∈ \mathcal{S}`` and ``X = (ν, u_1, …, u_d) ∈ T_p \mathcal{S}``. +The exponential map is given by + +````math + \operatorname{exp}_p(X) ≐ + \left( + \sqrt{(λ + ν)^2 + (λ m)^2},\\ + x_1 \cos\mathopen{\Big(} \frac{f \lVert u_1 \rVert_{x_1}}{m} \mathclose{\Big)} + \frac{u_1}{\lVert u_1 \rVert_{x_1}} \sin\mathopen{\Big(} \frac{f \lVert u_1 \rVert_{x_1}}{m} \mathclose{\Big)},\\ + …,\\ + x_d \cos\mathopen{\Big(} \frac{f \lVert u_d \rVert_{x_d}}{m} \mathclose{\Big)} + \frac{u_d}{\lVert u_d \rVert_{x_d}} \sin\mathopen{\Big(} \frac{f \lVert u_d \rVert_{x_d}}{m} \mathclose{\Big)} + \right), +```` + +where + +````math + \begin{aligned} + f &= \frac{\pi}{2} - \tan^{-1}\mathopen{\Big(} \frac{λ + ν}{λ m} \mathclose{\Big)},\\ + m &= \sqrt{\lVert u_1 \rVert_{x_1}^2 + … + \lVert u_d \rVert_{x_d}^2}. + \end{aligned} +```` + +If ``m = 0`` and ``-λ < ν``, then ``\operatorname{exp}_p(v) = p + X``. + +The formula is derived in proposition 3.1 in [JacobssonSwijsenVandervekenVannieuwenhoven:2024](@cite). +""" +exp(M::Segre{ℝ,V}, p, X) where {V} + +function exp!(::Segre{ℝ,V}, q, p, X) where {V} + m = sqrt( + sum([ + norm(Sphere(n - 1), x, xdot)^2 for (n, x, xdot) in zip(V, p[2:end], X[2:end]) + ]), + ) + + q[1][1] = sqrt((p[1][1] + X[1][1])^2 + (p[1][1] * m)^2) + + f = pi / 2 - atan((p[1][1] + X[1][1]) / (p[1][1] * m)) + if m == 0 + for (x, y) in zip(p[2:end], q[2:end]) + y .= x + end + else + for (n, x, y, xdot) in zip(V, p[2:end], q[2:end], X[2:end]) + a = norm(Sphere(n - 1), x, xdot) + y .= x * cos(a * f / m) .+ xdot * (f / m) * sinc(a * f / (m * pi)) + end + end + + return q +end + +@doc raw""" + get_coordinates(M::Segre{𝔽, V}, p, X, ::DefaultOrthonormalBasis; kwargs...) + +Get coordinates of `X` in the tangent space +``T_{(λ, x_1, …, x_d)} \mathcal{S} = \mathbb{R} × T_{x_1} \mathbb{S}^{n_1 - 1} ×…× T_{x_d} \mathbb{S}^{n_d - 1}`` +using a [`DefaultOrthonormalBasis`](@extref `ManifoldsBase.DefaultOrthonormalBasis`) on each factor. +""" +get_coordinates(M::Segre, p, X, ::DefaultOrthonormalBasis; kwargs...) + +function get_coordinates_orthonormal!( + M::Segre{ℝ,V}, + c, + p, + X, + ::RealNumbers; + kwargs..., +) where {V} + return c = vcat( + X[1], + p[1][1] * [ + get_coordinates(Sphere(n - 1), x, xdot, DefaultOrthonormalBasis(); kwargs...) for (n, x, xdot) in zip(V, p[2:end], X[2:end]) + ]..., + ) +end + +@doc raw""" + get_embedding(M::Segre{𝔽,V}) + +Return the embedding of the [`Segre`](@ref) manifold ``\mathcal{S}``, which is ``𝔽^{n_1 ×⋯× n_d}``. +""" +function get_embedding(::Segre{𝔽,V}) where {𝔽,V} + return Euclidean(prod(V)) +end + +@doc raw""" + get_vector( M::Segre{𝔽, V}, p, c, DefaultOrthonormalBasis; kwargs...) + +Get tangent vector `X` from coordinates in the tangent space +``T_{(λ, x_1, …, x_d)} \mathcal{S} = \mathbb{R} × T_{x_1} \mathbb{S}^{n_1 - 1} ×⋯× T_{x_d} \mathbb{S}^{n_d - 1}`` +using a [`DefaultOrthonormalBasis`](@extref `ManifoldsBase.DefaultOrthonormalBasis`) on each factor. +""" +get_vector(M::Segre, p, c, ::DefaultOrthonormalBasis; kwargs...) + +function get_vector_orthonormal!(M::Segre{ℝ,V}, X, p, c, ::RealNumbers; kwargs...) where {V} + c_ = deepcopy(c) + X[1] = [c_[1]] + c_ = c_[2:end] + for (i, n) in enumerate(V) + X[i + 1] = + get_vector( + Sphere(n - 1), + p[i + 1], + c_[1:(n - 1)], + DefaultOrthonormalBasis(); + kwargs..., + ) / p[1][1] + c_ = c_[n:end] + end + return X +end + +@doc raw""" + inner(M::Segre{ℝ, V}, p, X, Y,) + +Inner product between two tangent vectors ``X = (ν, u_1, …, u_d)`` and ``Y = (ξ, v_1, …, v_d)`` at ``p ≐ (λ, x_1, \dots, x_d)``. +This inner product is obtained by embedding the Segre manifold in the space of tensors equipped with the Euclidean metric: + +````math + \langle X, Y \rangle_{p} = \nu \xi + \lambda^2 (\langle u_1, v_1 \rangle_{x_1} + \dots + \langle u_d, v_d \rangle_{x_d}), +```` + +where ``ν, ξ ∈ T_{λ} ℝ^{+} = ℝ`` and ``u_i``, ``v_i ∈ T_{x_i} \mathbb{S}^{n_i - 1} ⊂ ℝ^{n_i}``. +""" +function inner(::Segre{ℝ}, p, X, Y) + return X[1][1] * Y[1][1] + p[1][1]^2 * dot(X[2:end], Y[2:end]) +end + +@doc raw""" + log(M::Segre{ℝ, V}, p, q) + +Logarithmic map on the Segre manifold. + +Let ``p ≐ (λ, x_1, …, x_d)``, ``q ≐ (μ, y_1, …, y_d) ∈ \mathcal{S}``. +Assume ``p`` and ``q`` are connected by a geodesic. + +Let + +````math + m = \sqrt{\sphericalangle(x_1, y_1)^2 + … + \sphericalangle(x_d, y_d)^2} +```` + +and assume ``(μ, y_1, …, y_d)`` is the representative of ``q`` that minimizes ``m``. Then + +````math + \operatorname{log}_p(q) = + \left( + \mu \cos{m} - \lambda, + (y_1 - ⟨x_1, y_1⟩ x_1) \frac{\mu \sphericalangle(x_1, y_1) \sin{m}}{\lambda m \sin{\sphericalangle(x_1, y_1)}}, + \dots, + (y_d - ⟨x_d, y_d⟩ x_d) \frac{\mu \sphericalangle(x_d, y_d) \sin{m}}{\lambda m \sin{\sphericalangle(x_d, y_d)}} + \right). +```` + +The formula is derived in theorem 4.4 in [JacobssonSwijsenVandervekenVannieuwenhoven:2024](@cite). +""" +log(M::Segre{ℝ,V}, p, q) where {V} + +function log!(M::Segre{ℝ,V}, X, p, q) where {V} + closest_representative!(M, q, p) + m = spherical_angle_sum(M, p, q) + + X[1][1] = q[1][1] * cos(m) - p[1][1] + for (n, xdot, x, y) in zip(V, X[2:end], p[2:end], q[2:end]) + a = distance(Sphere(n - 1), x, y) + xdot .= (y - dot(x, y) * x) * (q[1][1] / p[1][1]) * sinc(m / pi) / sinc(a / pi) + end + + return X +end + +manifold_dimension(::Segre{𝔽,V}) where {𝔽,V} = (1 + sum(V .- 1)) + +@doc raw""" + spherical_angle_sum(M::Segre{ℝ, V}, p, q) + +Let ``p ≐ (λ, x_1, …, x_d)``, ``q ≐ (μ, y_1, …, y_d) ∈ \mathcal{S}``. +Then this is + +````math + \sqrt{\sphericalangle(x_1, y_1)^2 + … + \sphericalangle(x_d, y_d)^2}, +```` + +where ``\sphericalangle(x_i, y_i)`` is the distance between ``x_i`` and ``y_i`` on the sphere ``\mathbb{S}^{n_i - 1}``. + +""" +function spherical_angle_sum(::Segre{ℝ,V}, p, q) where {V} + return sqrt( + sum([distance(Sphere(n - 1), x, y)^2 for (n, x, y) in zip(V, p[2:end], q[2:end])]), + ) +end + +@doc raw""" + rand(M::Segre{ℝ, V}; vector_at=nothing) + +If `vector_at` is `nothing`, return a random point on + +````math + ℝ^{+} × \mathbb{S}^{n_1 - 1} ×⋯× \mathbb{S}^{n_d - 1} +```` + +from a log-normal distribution on ``ℝ^{+}`` and a uniform distribution on ``\mathbb{S}^{n_1 - 1} ×⋯× \mathbb{S}^{n_d - 1}``. + +If `vector_at` is not `nothing`, return a random tangent vector from a normal distribution on the tangent space. +""" +function rand(M::Segre{ℝ,V}; vector_at=nothing, kwargs...) where {V} + if isnothing(vector_at) + return [ + rand(PositiveArrays(1); kwargs...), + [rand(Sphere(n - 1); kwargs...) for n in V]..., + ] + else + return [ + rand(PositiveArrays(1); vector_at=vector_at[1], kwargs...), + [ + rand(Sphere(n - 1); vector_at=xdot, kwargs...) for + (xdot, n) in zip(vector_at[2:end], V) + ]..., + ] + end +end + +@doc raw""" + riemann_tensor(M::Segre{ℝ, V}, p, X, Y, Z) + +Riemann tensor of the Segre manifold at ``p``. + +``\mathcal{S}`` is locally a warped product of ``ℝ^{+}`` and ``\mathbb{S}^{n_1 - 1} ×⋯× \mathbb{S}^{n_d - 1}``. +If ``p ≐ (λ, x_1, …, x_d) ∈ \mathcal{S}`` and ``X``, ``Y``, ``Z ∈ T_p (\mathbb{S}^{n_1 - 1} ×⋯× \mathbb{S}^{n_d - 1}) ⊂ T_p \mathcal{S}``, then + +````math + R_{\mathcal{S}}(X, Y) Z = R_{\mathbb{S}^{n_1 - 1} ×⋯× \mathbb{S}^{n_d - 1}}(X, Y) Z + λ^{-2}(⟨X, Z⟩_p Y - ⟨Y, Z⟩_p X). +```` + +``R_{\mathcal{S}}`` is zero in the remaining (orthogonal) directions. +""" +function riemann_tensor(M::Segre{ℝ,V}, p, X, Y, Z) where {V} + return [ + [0.0], + [ + riemann_tensor(Sphere(n - 1), x, xdot1, xdot2, xdot3) for + (n, x, xdot1, xdot2, xdot3) in zip(V, p[2:end], X[2:end], Y[2:end], Z[2:end]) + ]..., + ] + + (1 / p[1][1]^2) * ( + inner(M, p, [[0.0], X[2:end]...], [[0.0], Z[2:end]...]) * [[0.0], Y[2:end]...] - + inner(M, p, [[0.0], Y[2:end]...], [[0.0], Z[2:end]...]) * [[0.0], X[2:end]...] + ) +end + +@doc raw""" + sectional_curvature(M::Segre{ℝ, V}, p, u, v) + +Sectional curvature of the Segre manifold at ``p``. + +``\mathcal{S}`` is locally a warped product of ``ℝ^{+}`` and ``\mathbb{S}^{n_1 - 1} ×⋯× \mathbb{S}^{n_d - 1}`` +If ``p ≐ (λ, x_1, …, x_d) ∈ \mathcal{S}``, ``u_i ∈ T_{x_i} \mathbb{S}^{n_i - 1}``, and ``v_j ∈ T_{x_j} \mathbb{S}^{n_j - 1}``, then + +````math + K_{\mathcal{S}}(u_i, v_j) = \frac{\delta_{i j} - 1}{\lambda^2}. +```` + +``K_{\mathcal{S}}`` is zero in the remaining (orthogonal) directions. +""" +function sectional_curvature(M::Segre{ℝ,V}, p, X, Y) where {V} + return inner(M, p, riemann_tensor(M, p, X, Y, Y), X) / + (inner(M, p, X, X) * inner(M, p, Y, Y) - inner(M, p, X, Y)^2) +end + +function Base.show(io::IO, ::Segre{𝔽,V}) where {𝔽,V} + return print(io, "Segre($(join(V, ", ")); field=$(𝔽))") +end diff --git a/src/manifolds/SegreWarpedMetric.jl b/src/manifolds/SegreWarpedMetric.jl new file mode 100644 index 0000000000..ce53b3ea1b --- /dev/null +++ b/src/manifolds/SegreWarpedMetric.jl @@ -0,0 +1,317 @@ +@doc raw""" + WarpedMetric{A} <: AbstractMetric + +The ``A``-warped metric on the Segre manifold ``\mathcal{S}`` is a generalization of the Euclidean metric on ``\mathcal{S}``. +We denote this manifold by ``\mathcal{S}_A``. + +Similarly to ``\mathcal{S}``, when ``𝔽 = ℝ``, ``\mathcal{S}_A`` is a normal Riemannian covering of the product manifold + +````math + ℝ^{+} × \mathbb{S}^{n_1 - 1} ×⋯× \mathbb{S}^{n_d - 1} +```` + +with a [warped product metric](https://en.wikipedia.org/wiki/Warped_product), but the warping function now depends on the _warping factor_ ``A``. +``A = 1`` corresponds to the usual Segre manifold. +The Segre manifold is a cone in the sense that if ``p \in \mathcal{S}``, then ``r p \in \mathcal{S}`` for all ``r \neq 0``. The tangent subspace at ``p`` defined ``\mathrm{d} (r p) / \mathrm{d} r`` is called the _radial_ direction. ``A < 1`` puts less weight on the directions orthogonal to the radial direction compared to ``\mathcal{S}``, while ``A > 1`` puts more weight on those directions. + +The geometry is summarized in [JacobssonSwijsenVandervekenVannieuwenhoven:2024](@cite). + +# Constructor + WarpedMetric(A::Real) + +Generate a warped product metric with warping factor `A`. +""" +struct WarpedMetric{A} <: AbstractMetric end + +function WarpedMetric(A::Real) + return WarpedMetric{A}() +end + +function connected_by_geodesic( + M::MetricManifold{ℝ,Segre{ℝ,V},WarpedMetric{A}}, + p, + q, +) where {V,A} + return connected_by_geodesic(M.manifold, p, q) +end + +function closest_representative!( + M::MetricManifold{ℝ,Segre{ℝ,V},WarpedMetric{A}}, + q, + p, +) where {V,A} + return closest_representative!(M.manifold, q, p) +end + +@doc raw""" + distance(M::MetricManifold{ℝ, Segre{ℝ,V}, WarpedMetric{A}}, p, q) + +Riemannian distance between two points `p` and `q` on the warped Segre manifold. + +Assume ``p ≐ (λ, x_1,…, x_d)``, ``q ≐ (μ, y_1,…, y_d) ∈ \mathcal{S}_A`` are connected by a geodesic. Let + +````math + m = \sqrt{\sphericalangle(x_1, y_1)^2 +… + \sphericalangle(x_d, y_d)^2} +```` + +and assume ``(μ, y_1,…, y_d)`` is the representation of ``q`` that minimizes ``m``. Then + +````math + \operatorname{dist}_{\mathcal{S}_A}(p, q) = \sqrt{λ^2 - 2 λ μ \cos(A m) + μ^2}. +```` +""" +function distance(M::MetricManifold{ℝ,Segre{ℝ,V},WarpedMetric{A}}, p, q) where {V,A} + closest_representative!(M, q, p) + m = spherical_angle_sum(M, p, q) + return sqrt((p[1][1] - q[1][1])^2 + 4 * p[1][1] * q[1][1] * sin(A * m / 2)^2) + # Equivalent to sqrt(p[1][1]^2 + q[1][1]^2 - 2 * p[1][1] * q[1][1] * cos(A * m)) but more stable for small m +end + +@doc raw""" + exp(M::MetricManifold{ℝ, Segre{ℝ, V}, WarpedMetric{A}}, p, X) + +Exponential map on the warped Segre manifold. + +Let ``p ≐ (λ, x_1,…, x_d) ∈ \mathcal{S}_A`` and ``X = (ν, u_1,…, u_d) ∈ T_p \mathcal{S}_A``. +Then the exponential map is given by + +````math + \operatorname{exp}_p(X) ≐ + \left( + \sqrt{(λ + ν)^2 + (λ A m)^2},\\ + x_1 \cos\mathopen{\Big(} \frac{f \lVert u_1 \rVert_{x_1}}{A m} \mathclose{\Big)} + \frac{u_1}{\lVert u_1 \rVert_{x_1}} \sin\mathopen{\Big(} \frac{f \lVert u_1 \rVert_{x_1}}{A m} \mathclose{\Big)},\\ + …,\\ + x_d \cos\mathopen{\Big(} \frac{f \lVert u_d \rVert_{x_d}}{A m} \mathclose{\Big)} + \frac{u_d}{\lVert u_d \rVert_{x_d}} \sin\mathopen{\Big(} \frac{f \lVert u_d \rVert_{x_d}}{A m} \mathclose{\Big)} + \right), +```` + +where + +````math + \begin{aligned} + f &= \frac{\pi}{2} - \tan^{-1}\mathopen{\Big(} \frac{λ + ν}{λ A m} \mathclose{\Big)},\\ + m &= \sqrt{\lVert u_1 \rVert_{x_1}^2 + … + \lVert u_d \rVert_{x_d}^2}. + \end{aligned} +```` + +If ``m = 0`` and ``-λ < ν``, then ``\operatorname{exp}_p(v) = p + X``. + +The formula is derived in proposition 3.1 in [JacobssonSwijsenVandervekenVannieuwenhoven:2024](@cite). +""" +exp(M::MetricManifold{ℝ,Segre{ℝ,V},WarpedMetric{A}}, p, X) where {V,A} + +function exp!(::MetricManifold{ℝ,Segre{ℝ,V},WarpedMetric{A}}, q, p, X) where {V,A} + m = sqrt( + sum([ + norm(Sphere(n - 1), x, xdot)^2 for (n, x, xdot) in zip(V, p[2:end], X[2:end]) + ]), + ) + + q[1][1] = sqrt((p[1][1] + X[1][1])^2 + (p[1][1] * A * m)^2) + + f = pi / 2 - atan((p[1][1] + X[1][1]) / (p[1][1] * A * m)) + if m == 0 + for (x, y) in zip(p[2:end], q[2:end]) + y .= x + end + else + for (n, x, y, xdot) in zip(V, p[2:end], q[2:end], X[2:end]) + a = norm(Sphere(n - 1), x, xdot) + y .= + x * cos(a * f / (A * m)) .+ + xdot * (f / (A * m)) * sinc(a * f / (A * m * pi)) + end + end + + return q +end + +@doc raw""" + get_coordinates(M::Segre{𝔽, V}, p, v, ::DefaultOrthonormalBasis; kwargs...) + +Get coordinates of `X` in the tangent space ``T_{(λ, x_1,…, x_d)} \mathcal{S}_A = \mathbb{R} × T_{x_1} \mathbb{S}^{n_1 - 1} ×⋯× T_{x_d} \mathbb{S}^{n_d - 1}`` using a [`DefaultOrthonormalBasis`](@extref `ManifoldsBase.DefaultOrthonormalBasis`) on each factor. +""" +get_coordinates( + M::MetricManifold{𝔽,Segre{𝔽,V},WarpedMetric{A}}, + p, + X, + ::DefaultOrthonormalBasis; + kwargs..., +) where {𝔽,V,A} + +function get_coordinates_orthonormal!( + M::MetricManifold{𝔽,Segre{𝔽,V},WarpedMetric{A}}, + c, + p, + X, + ::RealNumbers; + kwargs..., +) where {𝔽,V,A} + return c = vcat( + X[1], + A * + p[1][1] * + [ + get_coordinates(Sphere(n - 1), x, xdot, DefaultOrthonormalBasis(); kwargs...) for (n, x, xdot) in zip(V, p[2:end], X[2:end]) + ]..., + ) +end + +@doc raw""" + get_vector(M::Segre{𝔽, V}, p, c, ::DefaultOrthonormalBasis; kwargs...) + +Get tangent vector `X` from coordinates in the tangent space ``T_{(λ, x_1,…, x_d)} \mathcal{S}_A = \mathbb{R} × T_{x_1} \mathbb{S}^{n_1 - 1} ×⋯× T_{x_d} \mathbb{S}^{n_d - 1}`` using a [`DefaultOrthonormalBasis`](@extref `ManifoldsBase.DefaultOrthonormalBasis`) on each factor. +""" +get_vector( + M::MetricManifold{𝔽,Segre{𝔽,V},WarpedMetric{A}}, + p, + c, + ::DefaultOrthonormalBasis; + kwargs..., +) where {V,A,𝔽} + +function get_vector_orthonormal!( + M::MetricManifold{ℝ,Segre{ℝ,V},WarpedMetric{A}}, + X, + p, + c, + ::RealNumbers; + kwargs..., +) where {V,A} + c_ = deepcopy(c) + X[1] = [c_[1]] + c_ = c_[2:end] + for (i, n) in enumerate(V) + X[i + 1] = + get_vector( + Sphere(n - 1), + p[i + 1], + c_[1:(n - 1)], + DefaultOrthonormalBasis(); + kwargs..., + ) / (A * p[1][1]) + c_ = c_[n:end] + end + + return X +end + +@doc raw""" + inner(M::MetricManifold{ℝ, Segre{ℝ, V}, WarpedMetric{A}}, p, X, Y) + +Inner product between two tangent vectors ``X = (ν, u_1,…, u_d)`` and ``Y = (ξ, v_1,…, v_d)`` at ``p \doteq (λ, x_1,…, x_d)``: +````math + ⟨X, Y⟩_{p} = ν ξ + (A λ)^2 (⟨ u_1, v_1 ⟩_{x_1} +… + ⟨u_d, v_d⟩_{x_d}), +```` +where ``ν``, ``ξ ∈ T_{λ} ℝ^{+} = ℝ`` and ``u_i``, ``v_i ∈ T_{x_i} \mathbb{S}^{n_i - 1} \subset ℝ^{n_i}``. +""" +function inner(M::MetricManifold{ℝ,Segre{ℝ,V},WarpedMetric{A}}, p, X, Y) where {V,A} + return X[1][1] * Y[1][1] + (A * p[1][1])^2 * dot(X[2:end], Y[2:end]) +end + +@doc raw""" + log(M::MetricManifold{ℝ, Segre{ℝ,V}, WarpedMetric{A}}, p, q) + +Logarithmic map on the warped Segre manifold. + +Let ``p ≐ (λ, x_1,…, x_d)``, ``q ≐ (μ, y_1,…, y_d) ∈ \mathcal{S}_A``. +Assume ``p`` and ``q`` are connected by a geodesic. +Let + +````math + m = \sqrt{\sphericalangle(x_1, y_1)^2 +… + \sphericalangle(x_d, y_d)^2} +```` + +and assume ``(μ, y_1,…, y_d)`` is the representative of ``q`` that minimizes ``m``. Then + +````math + \operatorname{log}_p(q) = + \left( + \mu \cos{m} - \lambda, + (y_1 - ⟨x_1, y_1⟩ x_1) \frac{\mu \sphericalangle(x_1, y_1) \sin(A m)}{\lambda A m \sin{\sphericalangle(x_1, y_1)}}, + \dots, + (y_d - ⟨x_d, y_d⟩ x_d) \frac{\mu \sphericalangle(x_d, y_d) \sin(A m)}{\lambda A m \sin{\sphericalangle(x_d, y_d)}} + \right). +```` + +The formula is derived in theorem 4.4 in [JacobssonSwijsenVandervekenVannieuwenhoven:2024](@cite). +""" +log(M::MetricManifold{ℝ,Segre{ℝ,V},WarpedMetric{A}}, p, q) where {V,A} + +function log!(M::MetricManifold{ℝ,Segre{ℝ,V},WarpedMetric{A}}, X, p, q) where {V,A} + closest_representative!(M, q, p) + m = spherical_angle_sum(M, p, q) + + X[1][1] = q[1][1] * cos(A * m) - p[1][1] + for (n, xdot, x, y) in zip(V, X[2:end], p[2:end], q[2:end]) + a = distance(Sphere(n - 1), x, y) + xdot .= (y - dot(x, y) * x) * (q[1][1] / p[1][1]) * sinc(A * m / pi) / sinc(a / pi) + end + + return X +end + +@doc raw""" + riemann_tensor(M::MetricManifold{ℝ, Segre{ℝ,V}, WarpedMetric{A}}, p, X, Y) + +Riemann tensor of the warped Segre manifold at ``p``. + +``\mathcal{S}_A`` is locally a warped product of ``ℝ^{+}`` and ``\mathbb{S}^{n_1 - 1} ×⋯× \mathbb{S}^{n_d - 1}``. If ``p ≐ (λ, x_1,…, x_d) ∈ \mathcal{S}_A`` and ``X``, ``Y``, ``Z ∈ T_{(x_1,…, x_d)} (\mathbb{S}^{n_1 - 1} ×⋯× \mathbb{S}^{n_d - 1}) \subset T_p \mathcal{S}_A`` then +````math + R_{\mathcal{S}_A}(X, Y) Z = R_{\mathbb{S}^{n_1 - 1} ×⋯× \mathbb{S}^{n_d - 1}}(X, Y) Z + λ^{-2}(⟨ X, Z ⟩_{p} Y - ⟨ Y, Z ⟩_{p} X). +```` +``R_{\mathcal{S}_A}`` is zero in the remaining (orthogonal) directions. +""" +function riemann_tensor( + M::MetricManifold{ℝ,Segre{ℝ,V},WarpedMetric{A}}, + p, + X, + Y, + Z, +) where {V,A} + return [ + [0.0], + [ + riemann_tensor(Sphere(n - 1), x, xdot1, xdot2, xdot3) for + (n, x, xdot1, xdot2, xdot3) in zip(V, p[2:end], X[2:end], Y[2:end], Z[2:end]) + ]..., + ] + + (1 / p[1][1]^2) * ( + inner(M, p, [[0.0], X[2:end]...], [[0.0], Z[2:end]...]) * [[0.0], Y[2:end]...] - + inner(M, p, [[0.0], Y[2:end]...], [[0.0], Z[2:end]...]) * [[0.0], X[2:end]...] + ) +end + +@doc raw""" + sectional_curvature(M::MetricManifold{ℝ, Segre{ℝ,V}, WarpedMetric{A}}, p, X, Y) + +Sectional curvature of the warped Segre manifold at ``p``. + +``\mathcal{S}_A`` is locally a warped product of ``ℝ^{+}`` and ``\mathbb{S}^{n_1 - 1} ×⋯× \mathbb{S}^{n_d - 1}`` +If ``p = (λ, x_1,…, x_d) ∈ \mathcal{S}``, ``u_i ∈ T_{x_i} \mathbb{S}^{n_i - 1}``, and ``v_j ∈ T_{x_j} \mathbb{S}^{n_j - 1}``, then +````math + K_{\mathcal{S}_A}(u_i, v_j) = \frac{A^{-2} \delta_{i j} - 1}{λ^{2}}. +```` +``K_{\mathcal{S}_A}`` is zero in the remaining (orthogonal) directions. +""" +function sectional_curvature( + M::MetricManifold{ℝ,Segre{ℝ,V},WarpedMetric{A}}, + p, + X, + Y, +) where {V,A} + return inner(M, p, riemann_tensor(M, p, X, Y, Y), X) / + (inner(M, p, X, X) * inner(M, p, Y, Y) - inner(M, p, X, Y)^2) +end + +function show(io::IO, ::WarpedMetric{A}) where {A} + return print(io, "WarpedMetric($A)") +end + +function spherical_angle_sum( + M::MetricManifold{ℝ,Segre{ℝ,V},WarpedMetric{A}}, + p, + q, +) where {V,A} + return spherical_angle_sum(M.manifold, p, q) +end diff --git a/test/manifolds/segre.jl b/test/manifolds/segre.jl new file mode 100644 index 0000000000..df30e8bdab --- /dev/null +++ b/test/manifolds/segre.jl @@ -0,0 +1,349 @@ +using Manifolds, Test, Random, LinearAlgebra, FiniteDifferences + +@testset "Segre Manifold" begin + # Manifolds to test + Ms = [ + Segre(10), + Segre(7, 2), + Segre(7, 9, 9), + Segre(9, 3, 6, 6), + MetricManifold(Segre(10), WarpedMetric(1.20)), + MetricManifold(Segre(2, 9), WarpedMetric(1.13)), + MetricManifold(Segre(9, 6, 10), WarpedMetric(0.87)), + MetricManifold(Segre(9, 3, 8, 10), WarpedMetric(1.40)), + ] + + # Vs[i] is the valence of Ms[i] + Vs = [(10,), (7, 2), (7, 9, 9), (9, 3, 6, 6), (10,), (2, 9), (9, 6, 10), (9, 3, 8, 10)] + + # n ≥ k, for same n,k X is in TpM and can be scaled by l + unit_p(n, k) = 1 / sqrt(k) .* [ones(k)..., zeros(n - k)...] + unit_X(n, k; l=1.0) = l / sqrt(n - k) .* [zeros(k)..., ones(n - k)...] + unit_c(n, k) = normalize([mod(k * i^i, 31) - 30 / 2 for i in 1:n]) # pseudo-rng + + # ps[i] is a point on Ms[i] + ps = [ + [[0.5], unit_p(10, 4)], + [[0.6], unit_p(7, 3), unit_p(2, 1)], + [[0.7], unit_p(7, 3), unit_p(9, 5), unit_p(9, 4)], + [[0.8], unit_p(9, 3), unit_p(3, 1), unit_p(6, 4), unit_p(6, 3)], + [[0.9], unit_p(10, 4)], + [[1.0], unit_p(2, 1), unit_p(9, 5)], + [[1.1], unit_p(9, 3), unit_p(6, 5), unit_p(10, 4)], + [[1.2], unit_p(9, 3), unit_p(3, 1), unit_p(8, 4), unit_p(10, 4)], + ] + + # qs[i] is a point on Ms[i] that is connected to ps[i] by a geodesic and uses the closest representative to ps[i] + # (tricked by only switching only one entry to zero) + qs = [ + [[0.1], unit_p(10, 3)], + [[0.2], unit_p(7, 2), unit_p(2, 2)], + [[0.3], unit_p(7, 2), unit_p(9, 4), unit_p(9, 3)], + [[0.4], unit_p(9, 3), unit_p(3, 1), unit_p(6, 3), unit_p(6, 2)], + [[0.5], unit_p(10, 3)], + [[0.6], unit_p(2, 2), unit_p(9, 5)], + [[0.7], unit_p(9, 2), unit_p(6, 4), unit_p(10, 3)], + [[0.8], unit_p(9, 2), unit_p(3, 2), unit_p(8, 3), unit_p(10, 3)], + ] + + # Xs[i] is a tangent vector to Ms[i] at ps[i] + Xs = [ + [[0.5], unit_X(10, 4)], + [[0.6], unit_X(7, 3), unit_X(2, 1)], + [[0.7], unit_X(7, 3), unit_X(9, 5), unit_X(9, 4)], + [[0.8], unit_X(9, 3), unit_X(3, 1), unit_X(6, 4), unit_X(6, 3)], + [[0.9], unit_X(10, 4)], + [[1.0], unit_X(2, 1), unit_X(9, 5)], + [[1.1], unit_X(9, 3), unit_X(6, 5), unit_X(10, 4)], + [[1.2], unit_X(9, 3), unit_X(3, 1), unit_X(8, 4), unit_X(10, 4)], + ] + + # Ys[i] is a tangent vector to Ms[i] at ps[i] such that exp(Ms[i], ps[i], t * vs[i]) is the closes representative to ps[i] for t in [-1, 1] + Ys = [ + [[0.5], unit_X(10, 5)], + [[0.6], unit_X(7, 4), unit_X(2, 1)], + [[0.7], unit_X(7, 4), unit_X(9, 6), unit_X(9, 5)], + [[0.8], unit_X(9, 4), unit_X(3, 1), unit_X(6, 5), unit_X(6, 4)], + [[0.9], unit_X(10, 5)], + [[1.0], unit_X(2, 1), unit_X(9, 6)], + [[1.1], unit_X(9, 4), unit_X(6, 5), unit_X(10, 5)], + [[1.2], unit_X(9, 4), unit_X(3, 1), unit_X(8, 5), unit_X(10, 5)], + ] + + # cs[i] is coordinates for a tangent vector at ps[i] + cs = [ + unit_c(10, 1), + unit_c(8, 2), + unit_c(23, 3), + unit_c(21, 4), + unit_c(10, 5), + unit_c(10, 6), + unit_c(23, 7), + unit_c(27, 8), + ] + + # When testing that exp(Ms[i], ps[i], t * Xs[i]) is an extremum of the length functional, we take a directional derivative along dcs[i] + dcs = [ + unit_c(3 * 10, 9), + unit_c(3 * 8, 10), + unit_c(3 * 23, 11), + unit_c(3 * 21, 12), + unit_c(3 * 10, 13), + unit_c(3 * 10, 14), + unit_c(3 * 23, 15), + unit_c(3 * 27, 16), + ] + + for (M, V, p, q, X, Y, c, dc) in zip(Ms, Vs, ps, qs, Xs, Ys, cs, dcs) + @testset "Manifold $M" begin + @testset "Segre" begin + get_manifold(::Segre{ℝ,V}) where {V} = Segre{ℝ,V}() + get_manifold(::MetricManifold{ℝ,Segre{ℝ,V},WarpedMetric{A}}) where {V,A} = + Segre{ℝ,V}() + @test Segre(V...) == get_manifold(M) + end + + @testset "manifold_dimension" begin + @test manifold_dimension(M) == 1 + sum(V .- 1) + end + + @testset "is_point" begin + @test is_point(M, p) + @test is_point(M, q) + @test_throws DomainError is_point( + M, + [[1.0, 0.0], p[2:end]...]; + error=:error, + ) + @test_throws DomainError is_point( + M, + [p[1], [1.0], p[3:end]...]; + error=:error, + ) + @test_throws DomainError is_point(M, [[-1.0], p[2:end]...]; error=:error) + @test_throws DomainError is_point(M, [p[1], 2 * p[2:end]...]; error=:error) + end + + @testset "is_vector" begin + @test is_vector(M, p, X; error=:error) + @test is_vector(M, p, Y; error=:error) + @test_throws DomainError is_vector( + M, + [[1.0, 0.0], p[2:end]...], + X, + false, + true, + ) + @test_throws DomainError is_vector( + M, + p, + [[1.0, 0.0], X[2:end]...], + false, + true, + ) + @test_throws DomainError is_vector(M, p, p, false, true) + end + + @testset "connected_by_geodesic" begin + @test connected_by_geodesic(M, p, q) + end + + Random.seed!(1) + @testset "rand" begin + @test is_point(M, rand(M)) + @test is_vector(M, p, rand(M, vector_at=p)) + end + + @testset "get_embedding" begin + @test get_embedding(M) == Euclidean(prod(V)) + end + + @testset "embed!" begin + # points + p_ = zeros(prod(V)) + p__ = zeros(prod(V)) + embed!(M, p_, p) + embed!(M, p__, [p[1], [-x for x in p[2:end]]...]) + @test is_point(get_embedding(M), p_) + @test isapprox(p_, (-1)^length(V) * p__) + + # vectors + X_ = zeros(prod(V)) + embed!(M, X_, p, X) + @test is_vector(get_embedding(M), p_, X_) + end + + @testset "get_coordinates" begin + @test isapprox(X, get_vector(M, p, get_coordinates(M, p, X))) + @test isapprox(c, get_coordinates(M, p, get_vector(M, p, c))) + + # Coordinates are ON + @test isapprox( + dot(c, get_coordinates(M, p, X)), + inner(M, p, X, get_vector(M, p, c)), + ) + end + + @testset "exp" begin + # Zero vector + p_ = exp(M, p, zeros.(size.(X))) + @test is_point(M, p_) + @test isapprox(p, p_; atol=1e-5) + + # Tangent vector in the scaling direction + p_ = exp(M, p, [X[1], zeros.(size.(X[2:end]))...]) + @test is_point(M, p_) + @test isapprox([p[1] + X[1], p[2:end]...], p_; atol=1e-5) + + # Generic tangent vector + p_ = exp(M, p, X) + @test is_point(M, p) + + geodesic_speed = + central_fdm(3, 1)(t -> distance(M, p, exp(M, p, t * X)), -1.0) + @test isapprox(geodesic_speed, -norm(M, p, X); atol=1e-5) + geodesic_speed = + central_fdm(3, 1)(t -> distance(M, p, exp(M, p, t * X)), -0.811) + @test isapprox(geodesic_speed, -norm(M, p, X); atol=1e-5) + geodesic_speed = + central_fdm(3, 1)(t -> distance(M, p, exp(M, p, t * X)), -0.479) + @test isapprox(geodesic_speed, -norm(M, p, X); atol=1e-5) + geodesic_speed = + central_fdm(3, 1)(t -> distance(M, p, exp(M, p, t * X)), 0.181) + @test isapprox(geodesic_speed, norm(M, p, X); atol=1e-5) + geodesic_speed = + central_fdm(3, 1)(t -> distance(M, p, exp(M, p, t * X)), 0.703) + @test isapprox(geodesic_speed, norm(M, p, X); atol=1e-5) + geodesic_speed = + central_fdm(3, 1)(t -> distance(M, p, exp(M, p, t * X)), 1.0) + @test isapprox(geodesic_speed, norm(M, p, X); atol=1e-5) + + # Geodesics are (locally) length-minizing. So let B_a be a one-parameter + # family of curves such that B_0 is a geodesic. Then the derivative of + # length(B_a) at a = 0 should be 0, and the second derivative should be + # nonnegative. + + n = manifold_dimension(M) + c0 = 0.0 * c + c1 = 0.25 * c + c2 = 0.5 * c + c3 = 0.75 * c + c4 = 1.0 * c + + function curve_length(d::Vector{Float64}) + @assert(length(d) == 3 * n) + + # Control points + d1 = d[1:n] + d2 = d[(n + 1):(2 * n)] + d3 = d[(2 * n + 1):(3 * n)] + + # Bezier curve from 0 to v with control points y1, ..., y4 + function b(t) + return ( + (1 - t)^4 * c0 + + 4 * t * (1 - t)^3 * (c1 + d1) + + 6 * t^2 * (1 - t)^2 * (c2 + d2) + + 4 * t^3 * (1 - t) * (c3 + d3) + + t^4 * c4 + ) + end + + # Length of curve on manifold + ps_ = [exp(M, p, get_vector(M, p, b(t))) for t in 0.0:1e-3:1.0] + ds = [ + distance(M, p1, p2) for + (p1, p2) in zip(ps_[1:(end - 1)], ps_[2:end]) + ] + return sum(ds) + end + + f = a -> curve_length(a * dc) + @test isapprox(central_fdm(3, 1)(f, 0.0), 0.0; atol=1e-5) + @test central_fdm(3, 2)(f, 0.0) >= 0.0 + end + + @testset "log" begin + # Same point + X_ = log(M, p, p) + @test is_vector(M, p, X_) + @test isapprox(zeros.(size.(X)), X_; atol=1e-5) + + # Scaled point + X_ = log(M, p, [q[1], p[2:end]...]) + @test is_vector(M, p, X_) + @test isapprox(X_, [q[1] - p[1], zeros.(size.(q[2:end]))...]; atol=1e-5) + + # Generic tangent vector + X_ = log(M, p, q) + @test is_vector(M, p, X_) + end + + @testset "norm" begin + @test isapprox(norm(M, p, log(M, p, q)), distance(M, p, q)) + end + + @testset "sectional_curvature" begin + # Test that sectional curvature is difference between circumference + # and 2 pi r for small circles. + + # Orthonormalize + X_ = X / norm(M, p, X) + Y_ = Y - inner(M, p, X_, Y) * X_ + Y_ = Y_ / norm(M, p, Y_) + + r = 1e-2 + ps_ = [ + exp(M, p, r * (cos(theta) * X_ + sin(theta) * Y_)) for + theta in 0.0:1e-3:(2 * pi) + ] + ds = [distance(M, p1, p2) for (p1, p2) in zip(ps_, [ps_[2:end]..., ps_[1]])] + C = sum(ds) + K = 3 * (2 * pi * r - C) / (pi * r^3) # https://en.wikipedia.org/wiki/Bertrand%E2%80%93Diguet%E2%80%93Puiseux_theorem + + @test isapprox(K, sectional_curvature(M, p, X, Y); rtol=1e-2, atol=1e-2) + end + + @testset "riemann_tensor" begin + @test isapprox( + sectional_curvature(M, p, X, Y), + inner(M, p, riemann_tensor(M, p, X, Y, Y), X) / + (inner(M, p, X, X) * inner(M, p, Y, Y) - inner(M, p, X, Y)^2), + ) + end + end + end + + # Test a point that does not use the closest representative + @testset "log" begin + M = Ms[4] + p = ps[4] + q = qs[4] + q_ = [q[1], q[2], q[3], q[4], -q[5]] + @test is_vector(M, p, log(M, p, q_)) + + M = Ms[8] + p = ps[8] + q = qs[8] + q_ = [q[1], q[2], q[3], q[4], -q[5]] + @test is_vector(M, p, log(M, p, q_)) + end + + # Test the formulas for sectional curvature stated in the docs + @testset "sectional_curvature" begin + M = Segre(3, 3) + p = [[1.3], [1.0, 0.0, 0.0], [1.0, 0.0, 0.0]] + X = [[0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]] + Y = [[0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]] + Z = [[0.0], [0.0, 0.0, 0.0], [0.0, 1.0, 0.0]] + V = [[1.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] + @test isapprox(0.0, sectional_curvature(M, p, X, Y)) + @test isapprox(-1 / p[1][1]^2, sectional_curvature(M, p, X, Z)) + @test isapprox(0.0, sectional_curvature(M, p, X, V)) + + M = MetricManifold(Segre(3, 3, 3), WarpedMetric(1.23)) + @test isapprox((1.23^(-2) - 1) / p[1][1]^2, sectional_curvature(M, p, X, Y)) + @test isapprox(-1 / p[1][1]^2, sectional_curvature(M, p, X, Z)) + @test isapprox(0.0, sectional_curvature(M, p, X, V)) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 799c5cbc74..3e6f8a0655 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -162,6 +162,7 @@ end include_test("manifolds/probability_simplex.jl") include_test("manifolds/projective_space.jl") include_test("manifolds/rotations.jl") + include_test("manifolds/segre.jl") include_test("manifolds/shape_space.jl") include_test("manifolds/skewhermitian.jl") include_test("manifolds/spectrahedron.jl") diff --git a/tutorials/_quarto.yml b/tutorials/_quarto.yml index 0966a12e28..4ae305952e 100644 --- a/tutorials/_quarto.yml +++ b/tutorials/_quarto.yml @@ -24,4 +24,4 @@ format: variant: -raw_html+tex_math_dollars+pipe_tables+footnotes wrap: preserve -jupyter: julia-1.10 +jupyter: julia-1.11