diff --git a/.gitignore b/.gitignore index de6ee81394..68db5183b7 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,4 @@ tutorials/_freeze accuracy/*_files accuracy/_freeze accuracy/.quarto +docs/src/tutorials/hand-gestures_files/figure-commonmark diff --git a/CITATION.bib b/CITATION.bib index 064943942a..ca8c186b08 100644 --- a/CITATION.bib +++ b/CITATION.bib @@ -1,9 +1,11 @@ -@online{2106.08777, +@article{AxenBaranBergmannRzecki:2023, AUTHOR = {Seth D. Axen and Mateusz Baran and Ronny Bergmann and Krzysztof Rzecki}, + EPRINT = {2021.08777}, + EPRINTTYPE = {arXiv}, + JOURNAL = {AMS Transactions on Mathematical Software}, + NOTE = {accepted for publication}, TITLE = {Manifolds.jl: An Extensible {J}ulia Framework for Data Analysis on Manifolds}, - YEAR = {2021}, - EPRINT = {2106.08777}, - EPRINTTYPE = {arXiv} + YEAR = {2023} } @softawre{manifoldsjl-zenodo-mostrecent, AUTHOR = {Seth D. Axen and Mateusz Baran and Ronny Bergmann}, diff --git a/Project.toml b/Project.toml index f90a728cee..1e78a71ff0 100644 --- a/Project.toml +++ b/Project.toml @@ -42,14 +42,15 @@ ManifoldsRecipesBaseExt = ["Colors", "RecipesBase"] ManifoldsTestExt = "Test" [compat] +BoundaryValueDiffEq = "3" Colors = "0.12" Distributions = "0.22.6, 0.23, 0.24, 0.25" Einsum = "0.4" Graphs = "1.4" HybridArrays = "0.4" Kronecker = "0.4, 0.5" -ManifoldDiff = "0.3.3" -ManifoldsBase = "0.14.1" +ManifoldDiff = "0.3.6" +ManifoldsBase = "0.14.11" MatrixEquations = "2.2" OrdinaryDiffEq = "6.31" Plots = "1" diff --git a/README.md b/README.md index 91cda761a2..3950e2dbcf 100644 --- a/README.md +++ b/README.md @@ -56,12 +56,14 @@ If you have any questions regarding the Manifolds.jl ecosystem feel free to reac If you use `Manifolds.jl` in your work, please cite the following ```biblatex -@online{2106.08777, - Author = {Seth D. Axen and Mateusz Baran and Ronny Bergmann and Krzysztof Rzecki}, - Title = {Manifolds.jl: An Extensible {J}ulia Framework for Data Analysis on Manifolds}, - Year = {2021}, - Eprint = {2106.08777}, - Eprinttype = {arXiv}, +@article{AxenBaranBergmannRzecki:2023, + AUTHOR = {Seth D. Axen and Mateusz Baran and Ronny Bergmann and Krzysztof Rzecki}, + EPRINT = {2021.08777}, + EPRINTTYPE = {arXiv}, + JOURNAL = {AMS Transactions on Mathematical Software}, + NOTE = {accepted for publication}, + TITLE = {Manifolds.jl: An Extensible {J}ulia Framework for Data Analysis on Manifolds}, + YEAR = {2023} } ``` diff --git a/docs/Project.toml b/docs/Project.toml index a253c731c4..c467dda81b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -5,6 +5,7 @@ CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" HybridArrays = "1baab800-613f-4b0a-84e4-9cd3431bfbb9" @@ -19,7 +20,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] -BoundaryValueDiffEq = "2" +BoundaryValueDiffEq = "3" CondaPkg = "0.2" DiffEqCallbacks = "2" Distributions = "0.22.6, 0.23, 0.24, 0.25" diff --git a/docs/make.jl b/docs/make.jl index 7417755f1c..9ab2b04778 100755 --- a/docs/make.jl +++ b/docs/make.jl @@ -34,6 +34,7 @@ end # (c) load necessary packages for the docs using Plots, RecipesBase, Manifolds, ManifoldsBase, Documenter, PythonPlot +using DocumenterCitations # required for loading methods that handle differential equation solving using OrdinaryDiffEq, BoundaryValueDiffEq, DiffEqCallbacks # required for loading the manifold tests functions @@ -61,7 +62,9 @@ open(joinpath(generated_path, "contributing.md"), "w") do io end # (e) ...finally! make docs +bib = CitationBibliography(joinpath(@__DIR__, "src", "references.bib"); style=:alpha) makedocs( + bib; # for development, we disable prettyurls format=Documenter.HTML(prettyurls=false, assets=["assets/favicon.ico"]), modules=[ @@ -159,6 +162,7 @@ makedocs( "Contributing" => "misc/contributing.md", "Internals" => "misc/internals.md", "Notation" => "misc/notation.md", + "References" => "misc/references.md", ], ], ) diff --git a/docs/src/features/statistics.md b/docs/src/features/statistics.md index 9a9af46a03..51df528995 100644 --- a/docs/src/features/statistics.md +++ b/docs/src/features/statistics.md @@ -7,3 +7,8 @@ Order = [:type, :function] ``` ## Literature + +```@bibliography +Pages = ["features/statistics.md"] +Canonical=false +``` \ No newline at end of file diff --git a/docs/src/manifolds/fixedrankmatrices.md b/docs/src/manifolds/fixedrankmatrices.md index 782ea62f94..c1fb92ab7a 100644 --- a/docs/src/manifolds/fixedrankmatrices.md +++ b/docs/src/manifolds/fixedrankmatrices.md @@ -7,3 +7,8 @@ Order = [:type, :function] ``` ## Literature + +```@bibliography +Pages = ["manifolds/fixedrankmatrices.md"] +Canonical=false +``` \ No newline at end of file diff --git a/docs/src/manifolds/grassmann.md b/docs/src/manifolds/grassmann.md index 0c3c346ebd..d0e964592f 100644 --- a/docs/src/manifolds/grassmann.md +++ b/docs/src/manifolds/grassmann.md @@ -21,3 +21,11 @@ Modules = [Manifolds] Pages = ["manifolds/GrassmannProjector.jl"] Order = [:type,:function] ``` + + +## Literature + +```@bibliography +Pages = ["manifolds/grassmann.md"] +Canonical=false +``` \ No newline at end of file diff --git a/docs/src/manifolds/group.md b/docs/src/manifolds/group.md index 4ae63967f9..67ae1a5fab 100644 --- a/docs/src/manifolds/group.md +++ b/docs/src/manifolds/group.md @@ -255,4 +255,4 @@ Order = [:type, :function] Modules = [Manifolds] Pages = ["groups/connections.jl"] Order = [:type, :function] -``` +``` \ No newline at end of file diff --git a/docs/src/manifolds/hyperbolic.md b/docs/src/manifolds/hyperbolic.md index 64ad617207..e3cb4bc089 100644 --- a/docs/src/manifolds/hyperbolic.md +++ b/docs/src/manifolds/hyperbolic.md @@ -176,3 +176,10 @@ And we can again look at the corresponding geodesics, for example plot!(scene, M, [pts[1], origin], geodesic_interpolation=100) plot!(scene, M, [pts[2], origin], geodesic_interpolation=100) ``` + +## Literature + +```@bibliography +Pages = ["manifolds/hyperbolic.md"] +Canonical=false +``` \ No newline at end of file diff --git a/docs/src/manifolds/probabilitysimplex.md b/docs/src/manifolds/probabilitysimplex.md index 19237f6491..182b838393 100644 --- a/docs/src/manifolds/probabilitysimplex.md +++ b/docs/src/manifolds/probabilitysimplex.md @@ -1,13 +1,15 @@ # The probability simplex ```@autodocs -Modules = [Manifolds] +Modules = [Manifolds, Base] Pages = ["manifolds/ProbabilitySimplex.jl"] Order = [:type, :function] Private=false Public=true ``` +## Euclidean metric + ```@autodocs Modules = [Manifolds] Pages = ["manifolds/ProbabilitySimplexEuclideanMetric.jl"] @@ -24,7 +26,7 @@ An isometric embedding of interior of [`ProbabilitySimplex`](@ref) in positive o This embedding isometrically maps the Fisher-Rao metric on the open probability simplex to the sphere of radius 1 with Euclidean metric. More details can be found in Section 2.2 -of [^AyJostLeSchwachhöfer2017]. +of [AyJostLeSchwachhoefer:2017](@cite). The name derives from the notion of probability amplitudes in quantum mechanics. They are complex-valued and their squared norm corresponds to probability. This construction @@ -39,3 +41,8 @@ Public=false ``` ## Literature + +```@bibliography +Pages = ["manifolds/probabilitysimplex.md"] +Canonical=false +``` \ No newline at end of file diff --git a/docs/src/manifolds/rotations.md b/docs/src/manifolds/rotations.md index 90a297aa6d..5eafa1f249 100644 --- a/docs/src/manifolds/rotations.md +++ b/docs/src/manifolds/rotations.md @@ -37,3 +37,8 @@ Order = [:type, :function] ``` ## Literature + +```@bibliography +Pages = ["manifolds/rotations.md"] +Canonical=false +``` diff --git a/docs/src/manifolds/sphere.md b/docs/src/manifolds/sphere.md index 8417c6f9e7..8db38c7e1d 100644 --- a/docs/src/manifolds/sphere.md +++ b/docs/src/manifolds/sphere.md @@ -62,3 +62,10 @@ p3 = 1/sqrt(3) .* [1.0, -1.0, 1.0] vecs = log.(Ref(M), pts2, Ref(p3)) plot!(scene, M, pts2, vecs; wireframe = false, linewidth=1.5) ``` + +## Literature + +```@bibliography +Pages = ["manifolds/sphere.md"] +Canonical=false +``` \ No newline at end of file diff --git a/docs/src/manifolds/stiefel.md b/docs/src/manifolds/stiefel.md index 98581b7678..ca988418e2 100644 --- a/docs/src/manifolds/stiefel.md +++ b/docs/src/manifolds/stiefel.md @@ -65,3 +65,8 @@ Private = true ``` ## Literature + +```@bibliography +Pages = ["manifolds/stiefel.md"] +Canonical=false +``` \ No newline at end of file diff --git a/docs/src/manifolds/symmetricpositivedefinite.md b/docs/src/manifolds/symmetricpositivedefinite.md index 1638bf148d..a2ffc10c9d 100644 --- a/docs/src/manifolds/symmetricpositivedefinite.md +++ b/docs/src/manifolds/symmetricpositivedefinite.md @@ -13,7 +13,7 @@ The manifold can be equipped with different metrics ## Common and metric independent functions ```@autodocs -Modules = [Manifolds] +Modules = [Manifolds, ManifoldsBase] Pages = ["manifolds/SymmetricPositiveDefinite.jl"] Order = [:function] Public=true diff --git a/docs/src/manifolds/symplectic.md b/docs/src/manifolds/symplectic.md index dabd1c6fd5..7f435c10b3 100644 --- a/docs/src/manifolds/symplectic.md +++ b/docs/src/manifolds/symplectic.md @@ -1,34 +1,34 @@ -# Symplectic +# Symplectic -The [`Symplectic`](@ref) manifold, denoted $\operatorname{Sp}(2n, \mathbb{F})$, is a closed, embedded, submanifold of +The [`Symplectic`](@ref) manifold, denoted $\operatorname{Sp}(2n, \mathbb{F})$, is a closed, embedded, submanifold of $\mathbb{F}^{2n \times 2n}$ that represents transformations into symplectic subspaces which keep the canonical symplectic form over $\mathbb{F}^{2n \times 2n }$ invariant under the standard embedding inner product. -The canonical symplectic form is a non-degenerate bilinear and skew symmetric map -$\omega\colon \mathbb{F}^{2n} \times \mathbb{F}^{2n} +The canonical symplectic form is a non-degenerate bilinear and skew symmetric map +$\omega\colon \mathbb{F}^{2n} \times \mathbb{F}^{2n} \rightarrow \mathbb{F}$, given by $\omega(x, y) = x^T Q_{2n} y$ for elements $x, y \in \mathbb{F}^{2n}$, with ````math - Q_{2n} = + Q_{2n} = \begin{bmatrix} 0_n & I_n \\ -I_n & 0_n \end{bmatrix}. -```` -That means that an element $p \in \operatorname{Sp}(2n)$ must fulfill the requirement that +```` +That means that an element $p \in \operatorname{Sp}(2n)$ must fulfill the requirement that ````math \omega (p x, p y) = x^T(p^TQp)y = x^TQy = \omega(x, y), ```` leading to the requirement on $p$ that $p^TQp = Q$. -The symplectic manifold also forms a group under matrix multiplication, called the $\textit{symplectic group}$. -Since all the symplectic matrices necessarily have determinant one, the [symplectic group](https://en.wikipedia.org/wiki/Symplectic_group) -$\operatorname{Sp}(2n, \mathbb{F})$ is a subgroup of the special linear group, $\operatorname{SL}(2n, \mathbb{F})$. When the underlying +The symplectic manifold also forms a group under matrix multiplication, called the $\textit{symplectic group}$. +Since all the symplectic matrices necessarily have determinant one, the [symplectic group](https://en.wikipedia.org/wiki/Symplectic_group) +$\operatorname{Sp}(2n, \mathbb{F})$ is a subgroup of the special linear group, $\operatorname{SL}(2n, \mathbb{F})$. When the underlying field is either $\mathbb{R}$ or $\mathbb{C}$ the symplectic group with a manifold structure constitutes a Lie group, with the Lie -Algebra +Algebra ````math \mathfrak{sp}(2n,F) = \{H \in \mathbb{F}^{2n \times 2n} \;|\; Q H + H^{T} Q = 0\}. ```` -This set is also known as the [Hamiltonian matrices](https://en.wikipedia.org/wiki/Hamiltonian_matrix), which have the +This set is also known as the [Hamiltonian matrices](https://en.wikipedia.org/wiki/Hamiltonian_matrix), which have the property that $(QH)^T = QH$ and are commonly used in physics. ```@autodocs @@ -38,3 +38,8 @@ Order = [:type, :function] ``` ## Literature + +```@bibliography +Pages = ["manifolds/symplectic.md"] +Canonical=false +``` \ No newline at end of file diff --git a/docs/src/manifolds/symplecticstiefel.md b/docs/src/manifolds/symplecticstiefel.md index 6f0bd71c7b..2fda3518e8 100644 --- a/docs/src/manifolds/symplecticstiefel.md +++ b/docs/src/manifolds/symplecticstiefel.md @@ -1,20 +1,20 @@ # Symplectic Stiefel -The [`SymplecticStiefel`](@ref) manifold, denoted $\operatorname{SpSt}(2n, 2k)$, -represents canonical symplectic bases of $2k$ dimensonal symplectic subspaces of $\mathbb{R}^{2n \times 2n}$. -This means that the columns of each element $p \in \operatorname{SpSt}(2n, 2k) \subset \mathbb{R}^{2n \times 2k}$ -constitute a canonical symplectic basis of $\operatorname{span}(p)$. -The canonical symplectic form is a non-degenerate, bilinear, and skew symmetric map -$\omega_{2k}\colon \mathbb{F}^{2k} \times \mathbb{F}^{2k} +The [`SymplecticStiefel`](@ref) manifold, denoted $\operatorname{SpSt}(2n, 2k)$, +represents canonical symplectic bases of $2k$ dimensonal symplectic subspaces of $\mathbb{R}^{2n \times 2n}$. +This means that the columns of each element $p \in \operatorname{SpSt}(2n, 2k) \subset \mathbb{R}^{2n \times 2k}$ +constitute a canonical symplectic basis of $\operatorname{span}(p)$. +The canonical symplectic form is a non-degenerate, bilinear, and skew symmetric map +$\omega_{2k}\colon \mathbb{F}^{2k} \times \mathbb{F}^{2k} \rightarrow \mathbb{F}$, given by $\omega_{2k}(x, y) = x^T Q_{2k} y$ for elements $x, y \in \mathbb{F}^{2k}$, with ````math - Q_{2k} = + Q_{2k} = \begin{bmatrix} 0_k & I_k \\ -I_k & 0_k \end{bmatrix}. -```` +```` Specifically given an element $p \in \operatorname{SpSt}(2n, 2k)$ we require that ````math \omega_{2n} (p x, p y) = x^T(p^TQ_{2n}p)y = x^TQ_{2k}y = \omega_{2k}(x, y) \;\forall\; x, y \in \mathbb{F}^{2k}, @@ -29,3 +29,8 @@ Order = [:type, :function] ``` ## Literature + +```@bibliography +Pages = ["manifolds/symplecticstiefel.md"] +Canonical=false +``` \ No newline at end of file diff --git a/docs/src/misc/internals.md b/docs/src/misc/internals.md index 3b0673f9de..3d5d10b2b4 100644 --- a/docs/src/misc/internals.md +++ b/docs/src/misc/internals.md @@ -14,9 +14,19 @@ Manifolds.nzsign Manifolds.realify Manifolds.realify! Manifolds.select_from_tuple +Manifolds.symmetrize +Manifolds.symmetrize! Manifolds.unrealify! Manifolds.usinc Manifolds.usinc_from_cos Manifolds.vec2skew! Manifolds.ziptuples ``` + +## Types in Extensions + +```@autodocs +Modules = [Manifolds] +Pages = ["../ext/ManifoldsOrdinaryDiffEqDiffEqCallbacksExt.jl"] +Order = [:type, :function] +``` \ No newline at end of file diff --git a/docs/src/misc/notation.md b/docs/src/misc/notation.md index 67fcf7259c..768c426c78 100644 --- a/docs/src/misc/notation.md +++ b/docs/src/misc/notation.md @@ -25,9 +25,12 @@ Within the documented functions, the utf8 symbols are used whenever possible, as | ``F`` | a fiber | | see [`VectorBundleFibers`](@ref) | | ``\mathbb F`` | a field, usually ``\mathbb F \in \{\mathbb R,\mathbb C, \mathbb H\}``, i.e. the real, complex, and quaternion numbers, respectively. | |field a manifold or a basis is based on | | ``\gamma`` | a geodesic | ``\gamma_{p;q}``, ``\gamma_{p,X}`` | connecting two points ``p,q`` or starting in ``p`` with velocity ``X``. | -| ``∇ f(p)`` | gradient of function ``f \colon \mathcal{M} \to \mathbb{R}`` at ``p \in \mathcal{M}`` | | | +| ``\operatorname{grad} f(p)`` | (Riemannian) gradient of function ``f \colon \mathcal{M} \to \mathbb{R}`` at ``p \in \mathcal{M}`` | | | +| ``\nabla f(p)`` | (Euclidean) gradient of function ``f \colon \mathcal{M} \to \mathbb{R}`` at ``p \in \mathcal{M}`` but thought of as evaluated in the embedding | `G` | | | ``\circ`` | a group operation | | | ``\cdot^\mathrm{H}`` | Hermitian or conjugate transposed for both complex or quaternion matrices| | +| ``\operatorname{Hess} f(p)`` | (Riemannian) Hessian of function ``f \colon T_p\mathcal{M} \to T_p\mathcal M`` (i.e. the 1-1-tensor form) at ``p \in \mathcal{M}`` | | | +| ``\nabla^2 f(p)`` | (Euclidean) Hessian of function ``f`` in the embedding | `H` | | | ``e`` | identity element of a group | | | ``I_k`` | identity matrix of size ``k\times k`` | | | ``k`` | indices | ``i,j`` | | @@ -37,6 +40,8 @@ Within the documented functions, the utf8 symbols are used whenever possible, as | ``\mathcal{G}`` | a (Lie) group | | | ``\log_p q`` | logarithmic map at ``p \in \mathcal M`` of a point ``q \in \mathcal M`` | ``\log_p(q)`` | | | ``\mathcal M`` | a manifold | ``\mathcal M_1, \mathcal M_2,\ldots,\mathcal N`` | | +| ``N_p \mathcal M`` | the normal space of the tangent space ``T_p \mathcal M`` in some embedding ``\mathcal E`` that should be clear from context | | | +| ``V`` | a normal vector from ``N_p \mathcal M`` | ``W`` | | | ``\operatorname{Exp}`` | the matrix exponential | | | ``\operatorname{Log}`` | the matrix logarithm | | | ``\mathcal P_{q\gets p}X`` | parallel transport | | of the vector ``X`` from ``T_p\mathcal M`` to ``T_q\mathcal M`` @@ -52,7 +57,8 @@ Within the documented functions, the utf8 symbols are used whenever possible, as | ``e_i \in \mathbb R^n`` | the ``i``th unit vector | ``e_i^n`` | the space dimension (``n``) is omited, when clear from context | ``B`` | a vector bundle | | | ``\mathcal T_{q\gets p}X`` | vector transport | | of the vector ``X`` from ``T_p\mathcal M`` to ``T_q\mathcal M`` -| ``\mathcal T_{p,Y}X`` | vector transport in direction ``Y`` | | of the vector ``X`` from ``T_p\mathcal M`` to ``T_q\mathcal M``, where ``q`` is deretmined by ``Y``, for example using the exponential map or some retraction. +| ``\mathcal T_{p,Y}X`` | vector transport in direction ``Y`` | | of the vector ``X`` from ``T_p\mathcal M`` to ``T_q\mathcal M``, where ``q`` is deretmined by ``Y``, for example using the exponential map or some retraction. | | ``\operatorname{Vol}(\mathcal M)`` | volume of manifold ``\mathcal M`` | | | ``\theta_p(X)`` | volume density for vector ``X`` tangent at point ``p`` | | +| ``\mathcal W`` | the Weingarten map ``\mathcal W: T_p\mathcal M × N_p\mathcal M → T_p\mathcal M`` | ``\mathcal W_p`` | the second notation to emphasize the dependency of the point ``p\in\mathcal M`` | | ``0_k`` | the ``k\times k`` zero matrix. | | diff --git a/docs/src/misc/references.md b/docs/src/misc/references.md new file mode 100644 index 0000000000..e96dc6abf4 --- /dev/null +++ b/docs/src/misc/references.md @@ -0,0 +1,10 @@ +# Literature + +We are slowly moving to using [`DocumenterCitations.jl`](https://github.com/JuliaDocs/DocumenterCitations.jl/). +The goal is to have all references used / mentioned in the documentation of Manifolds.jl also listed here. +If you notice a reference still defined in a footnote, please change it into a BibTeX reference and [open a PR](https://github.com/JuliaManifolds/Manifolds.jl/compare) + +Usually you will find a small reference section at the end of every documentation page that contains references for just that page. + +```@bibliography +``` \ No newline at end of file diff --git a/docs/src/references.bib b/docs/src/references.bib new file mode 100644 index 0000000000..ec6b78b79c --- /dev/null +++ b/docs/src/references.bib @@ -0,0 +1,821 @@ +# Manifold.jl +# +# Literature used within the Documentation of Manifolds.jl +# ======================================================== +# +# citekeys should be of the form AllAuthors:Year, unless there is really many authors. +# +# A +# ---------------------------------------------------------------------------------------- +@book{AbsilMahonySepulchre:2008, + AUTHOR = {Absil, P.-A. and Mahony, R. and Sepulchre, R.}, + DOI = {10.1515/9781400830244}, + NOTE = {available online at [press.princeton.edu/chapters/absil/](http://press.princeton.edu/chapters/absil/)}, + PUBLISHER = {Princeton University Press}, + TITLE = {Optimization Algorithms on Matrix Manifolds}, + YEAR = {2008}, +} +@incollection{AbsilMahonyTrumpf:2013, + AUTHOR = {P. -A. Absil and Robert Mahony and Jochen Trumpf}, + YEAR = {2013}, + DOI = {10.1007/978-3-642-40020-9_39}, + EDITOR = {Nielsen, Frank + and Barbaresco, Frédéric}, + ISBN = {978-3-642-40020-9}, + PUBLISHER = {Springer Berlin Heidelberg}, + PAGES = {361--368}, + TITLE = {An Extrinsic Look at the Riemannian Hessian}, + BOOKTITLE = {Geometric Science of Information}, +} +@article{AfsariTronVidal:2013, + DOI = {10.1137/12086282x}, + EPRINT = {1201.0925}, + EPRINTTYPE = {arXiv}, + YEAR = {2013}, + VOLUME = {51}, + NUMBER = {3}, + PAGES = {2230--2260}, + AUTHOR = {Bijan Afsari and Roberto Tron and René Vidal}, + TITLE = {On the Convergence of Gradient Descent for Finding the Riemannian Center of Mass}, + JOURNAL = {SIAM Journal on Control and Optimization} +} +@article{AndricaRohan:2013, + AUTHOR = {Andrica, D. and Rohan, R.-A.}, + ISSUE = {2}, + JOURNAL = {Balkan Journal of Geometry and Its Applications}, + PAGES = {1–10}, + TITLE = {Computing the Rodrigues coefficients of the exponential map of the Lie groups of matrices}, + URL = {https://www.emis.de/journals/BJGA/v18n2/B18-2-an.pdf}, + VOLUME = {18}, + YEAR = {2013}, +} +@article{AndruchowLarotondaRechtVarela:2014, + DOI = {10.1016/j.geomphys.2014.08.009}, + EPRINT = {1109.0520}, + EPRINTTYPE = {arXiv}, + YEAR = {2014}, + VOLUME = {86}, + PAGES = {241--257}, + AUTHOR = {E. Andruchow and G. Larotonda and L. Recht and A. Varela}, + TITLE = {The left invariant metric in the general linear group}, + JOURNAL = {Journal of Geometry and Physics} +} +@book{AyJostLeSchwachhoefer:2017, + DOI = {10.1007/978-3-319-56478-4}, + YEAR = {2017}, + AUTHOR = {Nihat Ay and Jürgen Jost and Hông Vân Lê and Lorenz Schwachhöfer}, + TITLE = {Information Geometry}, + PUBLISHER = {Springer Cham} +} +@article{AxenBaranBergmannRzecki:2023, + AUTHOR = {Seth D. Axen and Mateusz Baran and Ronny Bergmann and Krzysztof Rzecki}, + EPRINT = {2021.08777}, + EPRINTTYPE = {arXiv}, + JOURNAL = {AMS Transactions on Mathematical Software}, + NOTE = {accepted for publication}, + TITLE = {Manifolds.jl: An Extensible {J}ulia Framework for Data Analysis on Manifolds}, + YEAR = {2023} +} +# +# B +# ---------------------------------------------------------------------------------------- +@article{Bacak:2014, + AUTHOR = {Bačák, M.}, + DOI = {10.1137/140953393}, + EPRINT = {1210.2145}, + EPRINTTYPE = {arXiv}, + JOURNAL = {SIAM Journal on Optimization}, + NOTE = {arXiv: [1210.2145](https://arxiv.org/abs/1210.2145)}, + NUMBER = {3}, + PAGES = {1542--1566}, + TITLE = {Computing medians and means in Hadamard spaces}, + VOLUME = {24}, + YEAR = {2014} +} +@article{BendokatZimmermann:2021, + AUTHOR = {Bendokat, Thomas and Zimmermann, Ralf}, + Journal = {arXiv Preprint, 2108.12447}, + TITLE = {The real symplectic Stiefel and Grassmann manifolds: metrics, geodesics and applications}, + URL = {https://arxiv.org/abs/2108.12447}, + YEAR = {2021} +} +@article{BendokatZimmermannAbsil:2020, + AUTHOR = {Bendokat, Thomas and Zimmermann, Ralf and Absil, Pierre-Antoine}, + EPRINT = {2011.13699}, + EPRINTTYPE = {arXiv}, + JOURNAL = {arXiv Preprint}, + TITLE = {A Grassmann Manifold Handbook: Basic Geometry and Computational Aspects}, + URL = {https://arxiv.org/abs/2011.13699}, + YEAR = {2020} +} +@article{BergmannGousenbourger:2018, + AUTHOR = {Bergmann, Ronny and Gousenbourger, Pierre-Yves}, + DOI = {10.3389/fams.2018.00059}, + EPRINT = {1807.10090}, + EPRINTTYPE = {arXiv}, + JOURNAL = {Frontiers in Applied Mathematics and Statistics}, + TITLE = {A variational model for data fitting on manifolds by minimizing the acceleration of a Bézier curve}, + VOLUME = {4}, + YEAR = {2018} +} +@book{BinzPods:2008, + AUTHOR = {Biny, E and Pods, S}, + PUBLISHER = {American Mathematical Socienty}, + TITLE = {The Geometry of Heisenberg Groups: With Applications in Signal Theory, Optics, Quantization, and Field Quantization}, + YEAR = {2008} +} +@article{BirteaCaşuComănescu:2020, + DOI = {10.1007/s00605-020-01369-9}, + YEAR = {2020}, + VOLUME = {191}, + NUMBER = {3}, + PAGES = {465--485}, + AUTHOR = {Petre Birtea and Ioan Caçu and Dan Comănescu}, + TITLE = {Optimization on the real symplectic group}, + JOURNAL = {Monatshefte für Mathematik} +} +@article{BoyaSudarshanTilma:2003, + DOI = {10.1016/s0034-4877(03)80038-1}, + YEAR = 2003, + VOLUME = {52}, + NUMBER = {3}, + PAGES = {401--422}, + AUTHOR = {Luis J. Boya and E.C.G. Sudarshan and Todd Tilma}, + TITLE = {Volumes of compact manifolds}, + JOURNAL = {Reports on Mathematical Physics} +} + +@book{Boumal:2023, + TITLE = {An Introduction to Optimization on Smooth Manifolds}, + AUTHOR = {Boumal, Nicolas}, + YEAR = {2023}, + MONTH = mar, + EDITION = {First}, + PUBLISHER = {Cambridge University Press}, + DOI = {10.1017/9781009166164}, + URLDATE = {2023-07-13}, + ABSTRACT = {Optimization on Riemannian manifolds-the result of smooth geometry and optimization merging into one elegant modern framework-spans many areas of science and engineering, including machine learning, computer vision, signal processing, dynamical systems and scientific computing. This text introduces the differential geometry and Riemannian geometry concepts that will help students and researchers in applied mathematics, computer science and engineering gain a firm mathematical grounding to use these tools confidently in their research. Its charts-last approach will prove more intuitive from an optimizer's viewpoint, and all definitions and theorems are motivated to build time-tested optimization algorithms. Starting from first principles, the text goes on to cover current research on topics including worst-case complexity and geodesic convexity. Readers will appreciate the tricks of the trade for conducting research and for numerical implementations sprinkled throughout the book.}, + ISBN = {978-1-00-916616-4}, + URL = {https://www.nicolasboumal.net/#book} +} +@article{LeBrignantPuechmorel:2019, + DOI = {10.3390/e21010043}, + YEAR = 2019, + VOLUME = {21}, + NUMBER = {1}, + PAGES = {43}, + AUTHOR = {Alice Le Brigant and St{\'{e}}phane Puechmorel}, + TITLE = {Approximation of Densities on Riemannian Manifolds}, + JOURNAL = {Entropy} +} +# +# C +# ---------------------------------------------------------------------------------------- +@book{Chikuse:2003, + DOI = {10.1007/978-0-387-21540-2}, + YEAR = {2003}, + AUTHOR = {Yasuko Chikuse}, + PUBLISHER = {Springer New York}, + TITLE = {Statistics on Special Manifolds} +} +@inproceedings{ChakrabortyVemuri:2015, + AUTHOR = {Rudrasis Chakraborty and Baba C. Vemuri}, + BOOKTITLE = {2015 {IEEE} International Conference on Computer Vision (ICCV)}, + DOI = {10.1109/iccv.2015.481}, + TITLE = {Recursive Fréchet Mean Computation on the Grassmannian and Its Applications to Computer Vision}, + YEAR = {2015} +} +@article{ChakrabortyVemuri:2019, + DOI = {10.1214/18-aos1692}, + EPRINT = {1708.00045}, + EPRINTTYPE = {1708.00045}, + YEAR = {2019}, + VOLUME = {47}, + NUMBER = {1}, + AUTHOR = {Rudrasis Chakraborty and Baba C. Vemuri}, + TITLE = {Statistics on the Stiefel manifold: Theory and applications}, + JOURNAL = {The Annals of Statistics} +} +@incollection{ChengHoSalehianVemuri:2016, + AUTHOR = {Guang Cheng and Jeffrey Ho and Hesamoddin Salehian and Baba C. Vemuri}, + BOOKTITLE = {Riemannian Computing in Computer Vision}, + DOI = {10.1007/978-3-319-22957-7_2}, + PAGES = {21--43}, + PUBLISHER = {Springer, Cham}, + TITLE = {Recursive Computation of the Fréchet Mean on Non-positively Curved Riemannian Manifolds with Applications}, + YEAR = {2016} +} +@article{ChevallierKalungaAngulo:2017, + DOI = {10.1137/15m1053566}, + YEAR = 2017, + VOLUME = {10}, + NUMBER = {1}, + PAGES = {191--215}, + AUTHOR = {Emmanuel Chevallier and Emmanuel Kalunga and Jes{\'{u}}s Angulo}, + TITLE = {Kernel Density Estimation on Spaces of Gaussian Distributions and Symmetric Positive Definite Matrices}, + JOURNAL = {{SIAM} Journal on Imaging Sciences} +} +@article{ChevallierLiLuDunson:2022, + AUTHOR = {Chevallier, E. and Li, D. and Lu, Y. and Dunson, D. B.}, + JOURNAL = {ArXiv Preprint}, + TITLE = {Exponential-wrapped distributions on symmetric spaces}, + URL = {https://doi.org/10.48550/arXiv.2009.01983}, + YEAR = {2022} +} +# +# D +# ---------------------------------------------------------------------------------------- +@book{Devroye:1986, + DOI = {10.1007/978-1-4613-8643-8}, + YEAR = {1986}, + AUTHOR = {Luc Devroye}, + TITLE = {Non-Uniform Random Variate Generation}, + PUBLISHER = {Springer New York, NY} +} +@article{DewaeleBreidingVannieuwenhoven:2021, + AUTHOR = {Dewaele, Nick and Breiding, Paul and Vannieuwenhoven, Nick}, + EPRINT = {2106.13034}, + EPRINTTYPE = {arXiv}, + JOURNAL = {arXiv Preprint}, + TITLE = {The condition number of many tensor decompositions is invariant under Tucker compression}, + URL = {https://arxiv.org/abs/2106.13034}, + YEAR = {2021} +} +@article{DouikHassibi:2019, + DOI = {10.1109/tsp.2019.2946024}, + EPRINT = {1802.02628}, + YEAR = {2019}, + VOLUME = {67}, + NUMBER = {22}, + PAGES = {5761--5774}, + AUTHOR = {Ahmed Douik and Babak Hassibi}, + TITLE = {Manifold Optimization Over the Set of Doubly Stochastic Matrices: A Second-Order Geometry}, + JOURNAL = {IEEE Transactions on Signal Processing} +} +# +# E +# ---------------------------------------------------------------------------------------- +@article{EdelmanAriasSmith:1998, + DOI = {10.1137/s0895479895290954}, + EPRINT = {806030}, + EPRINTTYPE = {arXiv}, + YEAR = {1998}, + VOLUME = {20}, + NUMBER = {2}, + PAGES = {303--353}, + AUTHOR = {Alan Edelman and Tomás A. Arias and Steven T. Smith}, + TITLE = {The Geometry of Algorithms with Orthogonality Constraints}, + JOURNAL = {SIAM Journal on Matrix Analysis and Applications} +} +# +# F +# ---------------------------------------------------------------------------------------- +@article{FalorsideHaanDavidsonForre:2019, + AUTHOR = {Falorsi, Luca and de Haan, Pim and Davidson, Tim R. and Forré, Patrick}, + JOURNAL = {arXiv Preprint}, + URL = {https://arxiv.org/abs/1903.02958}, + TITLE = {Reparameterizing Distributions on Lie Groups}, + YEAR = {2019}, +} +@article{Fiori:2011, + DOI = {10.1137/100817115}, + YEAR = {2011}, + VOLUME = {32}, + NUMBER = {3}, + PAGES = {938--968}, + AUTHOR = {Simone Fiori}, + TITLE = {Solving Minimal-Distance Problems over the Manifold of Real-Symplectic Matrices}, + JOURNAL = {SIAM Journal on Matrix Analysis and Applications} +} +@inproceedings{FletcherVenkatasubramanianJoshi:2008, + DOI = {10.1109/cvpr.2008.4587747}, + YEAR = {2008}, + AUTHOR = {P. Thomas Fletcher and Suresh Venkatasubramanian and Sarang Joshi}, + TITLE = {Robust statistics on Riemannian manifolds via the geometric median}, + BOOKTITLE = {2008 {IEEE} Conference on Computer Vision and Pattern Recognition} +} +# +# G +# ---------------------------------------------------------------------------------------- +@article{GallierXu:2002, + AUTHOR = {Gallier, J and Xu, D}, + ISSUE = {4}, + JOURNAL = {International Journal of Robotics and Automation}, + PAGES = {1–11}, + TITLE = {Computing exponentials of skew-symmetric matrices and logarithms of orthogonal matrices}, + VOLUME = {17}, + URL = {http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.35.3205}, + YEAR = {2002} +} +@article{GaoSonAbsilStykel:2021, + DOI = {10.1137/20m1348522}, + YEAR = {2021}, + VOLUME = {31}, + NUMBER = {2}, + PAGES = {1546--1575}, + AUTHOR = {Bin Gao and Nguyen Thanh Son and P.-A. Absil and Tatjana Stykel}, + TITLE = {Riemannian Optimization on the Symplectic Stiefel Manifold}, + JOURNAL = {SIAM Journal on Optimization} +} +@incollection{GuiguiMaignantTroouvePennec:2021, + DOI = {10.1007/978-3-030-80209-7_12}, + YEAR = {2021}, + PAGES = {103--110}, + AUTHOR = {Nicolas Guigui and Elodie Maignant and Alain Trouv{\'{e}} and Xavier Pennec}, + TITLE = {Parallel Transport on Kendall Shape Spaces}, + BOOKTITLE = {Geometric Science of Information}, + PUBLISHER = {SPringer Cham} +} +# +# H +# ---------------------------------------------------------------------------------------- +@article{HanMishraJawanpuriaGao:2021, + AUTHOR = {Han, A and Mushra, B and Jawapanpuria, P and Gao, J}, + EPRINT = {2110.10464}, + EPRINTTYPE = {arXiv}, + JOURNAL = {arXive preprint}, + URL = {https://arxiv.org/abs/2110.10464}, + TITLE = {Learning with symmetric positive definite matrices via generalized Bures-Wasserstein geometry}, + YEAR = {2021}, +} +@inproceedings{HoChengSalehianVemuri:2013, + AUTHOR = {Ho, J. and Cheng, G. and Salehian, H. and Vemuri, B. C.}, + BOOKTITLE = {16th International Conference on Artificial Intelligence and Statistics}, + TITLE = {Recursive Karcher expectation estimators and geometric law of large numbers}, + YEAR = {2013}, + URL={http://proceedings.mlr.press/v31/ho13a.pdf} +} +@article{HosseiniUschmajew:2017, + DOI = {10.1137/16m1069298}, + YEAR = {2017}, + VOLUME = {27}, + NUMBER = {1}, + PAGES = {173--189}, + AUTHOR = {Seyedehsomayeh Hosseini and André Uschmajew}, + TITLE = {A Riemannian Gradient Sampling Algorithm for Nonsmooth Optimization on Manifolds}, + JOURNAL = {SIAM J. Optim.} +} +@article{HuangGallivanAbsil:2015, + DOI = {10.1137/140955483}, + YEAR = {2015}, + VOLUME = {25}, + NUMBER = {3}, + PAGES = {1660--1685}, + AUTHOR = {Wen Huang and K. A. Gallivan and P.-A. Absil}, + TITLE = {A Broyden Class of Quasi-Newton Methods for Riemannian Optimization}, + JOURNAL = {SIAM Journal on Optimization} +} +@article{HueperMarkinaSilvaLeite:2021, + DOI = {10.3934/jgm.2020031}, + YEAR = {2021}, + VOLUME = {13}, + NUMBER = {1}, + PAGES = {55}, + AUTHOR = {Knut Hüper and Irina Markina and F{\'{a}}tima Silva Leite}, + TITLE = {A Lagrangian approach to extremal curves on Stiefel manifolds}, + JOURNAL = {Journal of Geometric Mechanics} +} +# +# I +# ---------------------------------------------------------------------------------------- + +# +# J +# ---------------------------------------------------------------------------------------- +@article{JourneeBachAbsilSepulchre:2010, + DOI = {10.1137/080731359}, + EPRINT = {0807.4423}, + EPRINTTYPE = {arXiv}, + YEAR = {2010}, + VOLUME = {20}, + NUMBER = {5}, + PAGES = {2327--2351}, + AUTHOR = {M. Journée and F. Bach and P.-A. Absil and R. Sepulchre}, + TITLE = {Low-Rank Optimization on the Cone of Positive Semidefinite Matrices}, + JOURNAL = {SIAM Journal on Optimization} +} +# +# K +# ---------------------------------------------------------------------------------------- +@article{KanekoFioriTanaka:2013, + DOI = {10.1109/tsp.2012.2226167}, + YEAR = 2013, + VOLUME = {61}, + NUMBER = {4}, + PAGES = {883--894}, + AUTHOR = {Tetsuya Kaneko and Simone Fiori and Toshihisa Tanaka}, + TITLE = {Empirical Arithmetic Averaging Over the Compact Stiefel Manifold}, + JOURNAL = {IEEE Transactions on Signal Processing} +} +@article{Karcher:1977, + AUTHOR = {Karcher, H.}, + DOI = {10.1002/cpa.3160300502}, + JOURNAL = {Communications on Pure and Applied Mathematics}, + NUMBER = {5}, + PAGES = {509--541}, + TITLE = {Riemannian center of mass and mollifier smoothing}, + VOLUME = {30}, + YEAR = {1977} +} +@article{Kendall:1984, + DOI = {10.1112/blms/16.2.81}, + YEAR = {1984}, + VOLUME = {16}, + NUMBER = {2}, + PAGES = {81--121}, + AUTHOR = {David G. Kendall}, + TITLE = {Shape Manifolds, Procrustean Metrics, and Complex Projective Spaces}, + JOURNAL = {Bulletin of the London Mathematical Society} +} +@article{Kendall:1989, + DOI = {10.1214/ss/1177012582}, + YEAR = {1989}, + VOLUME = {4}, + NUMBER = {2}, + PAGES = {87–99}, + AUTHOR = {David G. Kendall}, + TITLE = {A Survey of the Statistical Theory of Shape}, + JOURNAL = {Statistical Sciences} +} +@article{KochLubich:2010, + DOI = {10.1137/09076578x}, + YEAR = {2010}, + VOLUME = {31}, + NUMBER = {5}, + PAGES = {2360--2375}, + AUTHOR = {Othmar Koch and Christian Lubich}, + TITLE = {Dynamical Tensor Approximation}, + JOURNAL = {{SIAM} Journal on Matrix Analysis and Applications} +} +@article{KressnerSteinlechnerVandereycken:2013, + DOI = {10.1007/s10543-013-0455-z}, + YEAR = 2013, + VOLUME = {54}, + NUMBER = {2}, + PAGES = {447--468}, + AUTHOR = {Daniel Kressner and Michael Steinlechner and Bart Vandereycken}, + TITLE = {Low-rank tensor completion by Riemannian optimization}, + JOURNAL = {{BIT} Numerical Mathematics} +} +# +# L +# ---------------------------------------------------------------------------------------- +julia> doi2bib("10.1080/10618600.2018.1549052"; abbreviate=false) +@article{LangreneMarin:2019, + DOI = {10.1080/10618600.2018.1549052}, + YEAR = 2019, + VOLUME = {28}, + NUMBER = {3}, + PAGES = {596--608}, + AUTHOR = {Nicolas Langren{\'{e}} and Xavier Warin}, + TITLE = {Fast and Stable Multivariate Kernel Density Estimation by Fast Sum Updating}, + JOURNAL = {Journal of Computational and Graphical Statistics} +} +@article{DeLathauwerDeMoorVanderwalle:2000, + DOI = {10.1137/s0895479896305696}, + YEAR = 2000, + VOLUME = {21}, + NUMBER = {4}, + PAGES = {1253--1278}, + AUTHOR = {Lieven De Lathauwer and Bart De Moor and Joos Vandewalle}, + TITLE = {A Multilinear Singular Value Decomposition}, + JOURNAL = {SIAM Journal on Matrix Analysis and Applications} +} +@article{Lin:2019, + AUTHOR = {Lin, Zhenhua}, + TITLE = {Riemannian Geometry of Symmetric Positive Definite Matrices via Cholesky Decomposition}, + JOURNAL = {SIAM Journal on Matrix Analysis and Applications}, + VOLUME = {40}, + NUMBER = {4}, + PAGES = {1353-1370}, + YEAR = {2019}, + DOI = {10.1137/18M1221084}, + EPRINT = {1908.09326}, + EPRINTTYPE = {arXiv}, +} +@article{LorenziPennec:2013, + DOI = {10.1007/s10851-013-0470-3}, + EPRINT = {00870489}, + EPRINTTYPE = {HAL}, + YEAR = {2013}, + VOLUME = {50}, + NUMBER = {1-2}, + PAGES = {5--17}, + AUTHOR = {Marco Lorenzi and Xavier Pennec}, + TITLE = {Efficient Parallel Transport of Deformations in Time Series of Images: From Schild's to Pole Ladder}, + JOURNAL = {Journal of Mathematical Imaging and Vision} +} +# +# M +# ---------------------------------------------------------------------------------------- +@article{MalagoMontruccioPistone:2018, + DOI = {10.1007/s41884-018-0014-4}, + YEAR = {2018}, + VOLUME = {1}, + NUMBER = {2}, + PAGES = {137--179}, + AUTHOR = {Luigi Malagó and Luigi Montrucchio and Giovanni Pistone}, + TITLE = {Wasserstein Riemannian geometry of Gaussian densities}, + JOURNAL = {Information Geometry} +} +@article{Marsaglia:1972, + DOI = {10.1214/aoms/1177692644}, + YEAR = {1972}, + VOLUME = {43}, + NUMBER = {2}, + PAGES = {645--646}, + AUTHOR = {George Marsaglia}, + TITLE = {Choosing a Point from the Surface of a Sphere}, + JOURNAL = {Annals of Mathematical Statistics} +} +@article{MartinNeff:2016, + DOI = {10.3934/jgm.2016010}, + EPRINT = {1409.7849}, + EPRINTTYPE = {arXiv}, + YEAR = {2016}, + VOLUME = {8}, + NUMBER = {3}, + PAGES = {323--357}, + AUTHOR = {Patrizio Neff and Robert J. Martin}, + TITLE = {Minimal geodesics on {GL}(n) for left-invariant, right-O(n)-invariant Riemannian metrics}, + JOURNAL = {J. Geom. Mech.} +} +@article{MassartAbsil:2020, + DOI = {10.1137/18m1231389}, + YEAR = {2020}, + VOLUME = {41}, + NOTE = {Preprint: [sites.uclouvain.be/absil/2018.06](https://sites.uclouvain.be/absil/2018.06)}, + NUMBER = {1}, + PAGES = {171--198}, + AUTHOR = {Estelle Massart and P.-A. Absil}, + TITLE = {Quotient Geometry with Simple Geodesics for the Manifold of Fixed-Rank Positive-Semidefinite Matrices}, + JOURNAL = {SIAM Journal on Matrix Analysis and Applications} +} +@inproceedings{MuralidharanFlecther:2012, + DOI = {10.1109/cvpr.2012.6247780}, + YEAR = 2012, + AUTHOR = {P. Muralidharan and P. T. Fletcher}, + TITLE = {Sasaki metrics for analysis of longitudinal data on manifolds}, + BOOKTITLE = {2012 {IEEE} Conference on Computer Vision and Pattern Recognition} +} + +# +# N +# ---------------------------------------------------------------------------------------- +@article{Nguyen:2023, + DOI = {10.1007/s10957-023-02242-z}, + EPRINT = {2009.10159}, + EPRINTTYPE = {arXiv}, + YEAR = {2023}, + MONTH = jun, + PUBLISHER = {Springer Science and Business Media LLC}, + VOLUME = {198}, + NUMBER = {1}, + PAGES = {135--164}, + AUTHOR = {Du Nguyen}, + TITLE = {Operator-Valued Formulas for Riemannian Gradient and Hessian and Families of Tractable Metrics in Riemannian Optimization}, + JOURNAL = {Journal of Optimization Theory and Applications} +} +# +# P +# ---------------------------------------------------------------------------------------- +@article{Pennec:2006, + DOI = {10.1007/s10851-006-6228-4}, + YEAR = {2006}, + VOLUME = {25}, + NUMBER = {1}, + PAGES = {127--154}, + AUTHOR = {Xavier Pennec}, + TITLE = {Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric Measurements}, + JOURNAL = {Journal of Mathematical Imaging and Vision} +} +@article{Pennec:2018, + AUTHOR ={Xavier Pennec}, + EPRINT = {1805.11436}, + EPRINTTYPE = {arXiv}, + Journal = {arXiv Preprint}, + TITLE = {Parallel Transport with Pole Ladder: a Third Order Scheme in Affine Connection Spaces which is Exact in Affine Symmetric Spaces.}, + YEAR = {2018}, + URL = {https://arxiv.org/abs/1805.11436}, +} +@incollection{PennecArsigny:2012, + DOI = {10.1007/978-3-642-30232-9_7}, + EPRINT = {00699361}, + EPRINTTYPE = {HAL}, + YEAR = {2012}, + PAGES = {123--166}, + AUTHOR = {Xavier Pennec and Vincent Arsigny}, + TITLE = {Exponential Barycenters of the Canonical Cartan Connection and Invariant Means on Lie Groups}, + BOOKTITLE = {Matrix Information Geometry}, + PUBLISHER = {Springer, Berlin, Heidelberg} +} +@incollection{PennecLorenzi:2020, + DOI = {10.1016/b978-0-12-814725-2.00012-1}, + EDITORS = {X. Pennec and S. Sommer and T. Fletcher}, + YEAR = {2020}, + PAGES = {169--229}, + AUTHOR = {Xavier Pennec and Marco Lorenzi}, + TITLE = {Beyond Riemannian geometry: The affine connection setting for transformation groups}, + BOOKTITLE = {Riemannian Geometric Statistics in Medical Image Analysis}, + PUBLISHER = {Elsevier} +} +# +# R +# ---------------------------------------------------------------------------------------- +@inproceedings{Rentmeesters:2011, + DOI = {10.1109/cdc.2011.6161280}, + YEAR = {2011}, + AUTHOR = {Quentin Rentmeesters}, + PAGES = {7141–7146}, + TITLE = {A gradient method for geodesic data fitting on some symmetric Riemannian manifolds}, + BOOKTITLE = {IEEE Conference on Decision and Control and European Control Conference} +} +@phdthesis{RicoMartinez:1988, + AUTHOR = {Rico Martinez, Jose Maria}, + SCHOOL = {University of FLorida}, + TITLE = {Representations of the Euclidean group and its applications to the kinematics of spatial chains}, + URL = {https://ufdc.ufl.edu/aa00040974/00001}, + YEAR = {1988} +} +# +# S +# ---------------------------------------------------------------------------------------- +@inproceedings{SalehianEtAl:2015, + AUTHOR = {Salehian, Hesamoddin and Chakraborty, Rudrasis and and Ofori, Edward and Vaillancourt, David and Vemuri, Baba C.}, + BOOKTITLE = {5th MICCAI workshop on Mathematical Foundations of Computational Anatomy}, + MONTH = {10}, + TITLE = {An efficient recursive estimator of the Fréchet mean on hypersphere with applications to Medical Image Analysis}, + URL = {https://www-sop.inria.fr/asclepios/events/MFCA15/Papers/MFCA15_4_2.pdf}, + YEAR = {2015}, +} +@article{Sasaki:1958, + DOI = {10.2748/tmj/1178244668}, + YEAR = {1958}, + VOLUME = {10}, + NUMBER = {3}, + AUTHOR = {Shigeo Sasaki}, + TITLE = {On the differential geometry of tangent bundles of Riemannian manifolds}, + JOURNAL = {Tohoku Math. J.} +} +@book{Suhubi:2013, + AUTHOR = {Suhubi, E}, + PUBLISHER = {Academic Press}, + TITLE = {Exterior Analysis: Using Applications of Differential Forms}, + YEAR = {2013} +} +@book{SrivastavaKlassen:2016, + AUTHOR = {Anuj Srivastava and Eric P. Klassen}, + DOI = {10.1007/978-1-4939-4020-2}, + ISBN = {978-1-4939-4018-9}, + PUBLISHER = {Springer NEw York}, + TITLE = {Functional and Shape Data Analysis}, + YEAR = {2016} +} +# +# T +# ---------------------------------------------------------------------------------------- +@misc{Tornier:2020, + AUTHOR = {Tornier, Stephan}, + TITLE = {Haar Measures}, + JOURNAL = {arXiv Preprint}, + URL = {https://arxiv.org/abs/2006.10956}, + YEAR = {2020}, +} + +@article{TronDaniilidis:2017, + DOI = {10.1137/16m1091332}, + YEAR = {2017}, + VOLUME = {10}, + NUMBER = {3}, + PAGES = {1416--1445}, + AUTHOR = {Roberto Tron and Kostas Daniilidis}, + TITLE = {The Space of Essential Matrices as a Riemannian Quotient Manifold}, + JOURNAL = {SIAM J. Imaging Sci.} +} +# +# V +# ---------------------------------------------------------------------------------------- +@article{Vandereycken:2013, + AUTHOR = {Vandereycken, Bart}, + PUBLISHER = {Society for Industrial \& Applied Mathematics (SIAM)}, + YEAR = {2013}, + DOI = {10.1137/110845768}, + JOURNAL = {SIAM Journal on Optimization}, + NUMBER = {2}, + PAGES = {1214--1236}, + TITLE = {Low-rank matrix completion by Riemannian optimization}, + VOLUME = {23} +} +@article{VannieuwenhovenVanderbrilMeerbergen:2012, + DOI = {10.1137/110836067}, + YEAR = {2012}, + VOLUME = {34}, + NUMBER = {2}, + PAGES = {A1027--A1052}, + AUTHOR = {Nick Vannieuwenhoven and Raf Vandebril and Karl Meerbergen}, + TITLE = {A New Truncation Strategy for the Higher-Order Singular Value Decomposition}, + JOURNAL = {SIAM Journal on Scientific Computing} +} + +# +# W +# ---------------------------------------------------------------------------------------- +@article{WangSunFiori:2018, + DOI = {10.1002/mma.4890}, + YEAR = {2018}, + VOLUME = {41}, + NUMBER = {11}, + PAGES = {4273--4286}, + AUTHOR = {Jing Wang and Huafei Sun and Simone Fiori}, + TITLE = {A Riemannian-steepest-descent approach for optimization on the real symplectic group}, + JOURNAL = {Mathematical Methods in the Applied Science} +} +@article{West:1979, + DOI = {10.1145/359146.359153}, + YEAR = {1979}, + VOLUME = {22}, + NUMBER = {9}, + PAGES = {532--535}, + AUTHOR = {D. H. D. West}, + TITLE = {Updating mean and variance estimates}, + JOURNAL = {Communications of the ACM} +} +# +# X +# ---------------------------------------------------------------------------------------- + +# +# Y +# ---------------------------------------------------------------------------------------- +@article{YeWongLim:2021, + DOI = {10.1007/s10107-021-01640-3}, + YEAR = {2021}, + VOLUME = {194}, + NUMBER = {1-2}, + PAGES = {621--660}, + AUTHOR = {Ke Ye and Ken Sze-Wai Wong and Lek-Heng Lim}, + TITLE = {Optimization on flag manifolds}, + JOURNAL = {Mathematical Programming} +} +# +# Z +# ---------------------------------------------------------------------------------------- +@article{Zhu:2016, + DOI = {10.1007/s10589-016-9883-4}, + YEAR = {2016}, + VOLUME = {67}, + NUMBER = {1}, + PAGES = {73--110}, + AUTHOR = {Xiaojing Zhu}, + TITLE = {A Riemannian conjugate gradient method for optimization on the Stiefel manifold}, + JOURNAL = {Computational Optimization and Applications} +} +@article{ZhuDuan:2018, + DOI = {10.1007/s11590-018-1341-z}, + YEAR = {2018}, + VOLUME = {13}, + NUMBER = {5}, + PAGES = {1069--1083}, + AUTHOR = {Xiaojing Zhu and Chunyan Duan}, + TITLE = {On matrix exponentials and their approximations related to optimization on the Stiefel manifold}, + JOURNAL = {Optimization Letters} +} +@article{Zimmermann:2017, + DOI = {10.1137/16m1074485}, + EPRINT = {1604.05054}, + EPRINTTYPE = {arXiv}, + YEAR = {2017}, + VOLUME = {38}, + NUMBER = {2}, + PAGES = {322--342}, + AUTHOR = {Ralf Zimmermann}, + TITLE = {A Matrix-Algebraic Algorithm for the Riemannian Logarithm on the Stiefel Manifold under the Canonical Metric}, + JOURNAL = {SIAM J. Matrix Anal. Appl.} +} +@article{ZimmermannHueper:2022, + AUTHOR = {Zimmermann, Ralf and Hüper, Knut}, + TITLE = {Computing the Riemannian Logarithm on the Stiefel Manifold: Metrics, Methods, and Performance}, + JOURNAL = {SIAM Journal on Matrix Analysis and Applications}, + VOLUME = {43}, + NUMBER = {2}, + PAGES = {953-980}, + YEAR = {2022}, + DOI = {10.1137/21M1425426}, + EPRINT = {2103.12046}, + EPRINTTYPE = {arXiv} +} +# +# Å +# ---------------------------------------------------------------------------------------- +@article{AastroemPetraSchmitzerSchnoerr:2017, + DOI = {10.1007/s10851-016-0702-4}, + EPRINT = {1603.05285}, + EPRINTTYPE = {arXiv}, + YEAR = {2017}, + VOLUME = {58}, + NUMBER = {2}, + PAGES = {211--238}, + 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/src/Manifolds.jl b/src/Manifolds.jl index cb0001b925..259c4a7b17 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -165,6 +165,8 @@ import ManifoldsBase: _vector_transport_to!, vee, vee!, + Weingarten, + Weingarten!, zero_vector, zero_vector!, CotangentSpace, @@ -176,7 +178,9 @@ import ManifoldDiff: jacobi_field, jacobi_field!, riemannian_gradient, - riemannian_gradient! + riemannian_gradient!, + riemannian_Hessian, + riemannian_Hessian! using Base.Iterators: repeated using Distributions @@ -481,7 +485,7 @@ Base.in(p, M::AbstractManifold; kwargs...) = is_point(M, p, false; kwargs...) X ∈ TangentSpaceAtPoint(M,p) Check whether `X` is a tangent vector from (in) the tangent space $T_p\mathcal M$, i.e. -the [`TangentSpaceAtPoint`](@ref Manifolds.TangentSpaceAtPoint) at `p` on the [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold) `M`. +the [`TangentSpaceAtPoint`](@ref) at `p` on the [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold) `M`. This method uses [`is_vector`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions.html#ManifoldsBase.is_vector) deactivating the error throw option. """ function Base.in(X, TpM::TangentSpaceAtPoint; kwargs...) @@ -493,12 +497,7 @@ end Volume of manifold `M` defined through integration of Riemannian volume element in a chart. Note that for many manifolds there is no universal agreement over the exact ranges over -which the integration should happen. For details see [^BoyaSudarshanTilma2003]. - -[^BoyaSudarshanTilma2003]: - > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports - > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003, - > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1) +which the integration should happen. For details see [BoyaSudarshanTilma:2003](@cite). """ manifold_volume(::AbstractManifold) @@ -508,14 +507,9 @@ manifold_volume(::AbstractManifold) Volume density function of manifold `M`, i.e. determinant of the differential of exponential map `exp(M, p, X)`. Determinant can be understood as computed in a basis, from the matrix of the linear operator said differential corresponds to. Details are available in Section 4.1 -of [^ChevallierLiLuDunson2022]. +of [ChevallierLiLuDunson:2022](@cite). Note that volume density is well-defined only for `X` for which `exp(M, p, X)` is injective. - -[^ChevallierLiLuDunson2022]: - > E. Chevallier, D. Li, Y. Lu, and D. B. Dunson, “Exponential-wrapped distributions on - > symmetric spaces.” arXiv, Oct. 09, 2022. - > doi: [10.48550/arXiv.2009.01983](https://doi.org/10.48550/arXiv.2009.01983). """ volume_density(::AbstractManifold, p, X) @@ -857,6 +851,8 @@ export ×, retract!, riemannian_gradient, riemannian_gradient!, + riemannian_Hessian, + riemannian_Hessian!, riemann_tensor, riemann_tensor!, set_component!, @@ -886,6 +882,8 @@ export ×, vertical_component, vertical_component!, volume_density, + Weingarten, + Weingarten!, zero_vector, zero_vector! # Lie group types & functions diff --git a/src/groups/connections.jl b/src/groups/connections.jl index 5f8f5acd93..8b68ff6c2e 100644 --- a/src/groups/connections.jl +++ b/src/groups/connections.jl @@ -3,13 +3,7 @@ AbstractCartanSchoutenConnection Abstract type for Cartan-Schouten connections, that is connections whose geodesics -going through group identity are one-parameter subgroups. See[^Pennec2020] for details. - -[^Pennec2020]: - > X. Pennec and M. Lorenzi, “5 - Beyond Riemannian geometry: The affine connection - > setting for transformation groups,” in Riemannian Geometric Statistics in Medical Image - > Analysis, X. Pennec, S. Sommer, and T. Fletcher, Eds. Academic Press, 2020, pp. 169–229. - > doi: [10.1016/B978-0-12-814725-2.00012-1](https://doi.org/10.1016/B978-0-12-814725-2.00012-1). +going through group identity are one-parameter subgroups. See [PennecLorenzi:2020](@cite) for details. """ abstract type AbstractCartanSchoutenConnection <: AbstractAffineConnection end @@ -50,7 +44,7 @@ const CartanSchoutenZeroGroup{𝔽,M} = ConnectionManifold{𝔽,M,CartanSchouten exp(M::ConnectionManifold{𝔽,<:AbstractDecoratorManifold{𝔽},<:AbstractCartanSchoutenConnection}, p, X) where {𝔽} Compute the exponential map on the [`ConnectionManifold`](@ref) `M` with a Cartan-Schouten -connection. See Sections 5.3.2 and 5.3.3 of [^Pennec2020] for details. +connection. See Sections 5.3.2 and 5.3.3 of [PennecLorenzi:2020](@cite) for details. """ function exp( M::ConnectionManifold{ @@ -109,7 +103,7 @@ end log(M::ConnectionManifold{𝔽,<:AbstractDecoratorManifold{𝔽},<:AbstractCartanSchoutenConnection}, p, q) where {𝔽} Compute the logarithmic map on the [`ConnectionManifold`](@ref) `M` with a Cartan-Schouten -connection. See Sections 5.3.2 and 5.3.3 of [^Pennec2020] for details. +connection. See Sections 5.3.2 and 5.3.3 of [PennecLorenzi:2020](@cite) for details. """ function log( M::ConnectionManifold{ @@ -144,7 +138,7 @@ end parallel_transport_to(M::CartanSchoutenMinusGroup, p, X, q) Transport tangent vector `X` at point `p` on the group manifold `M` with the -[`CartanSchoutenMinus`](@ref) connection to point `q`. See [^Pennec2020] for details. +[`CartanSchoutenMinus`](@ref) connection to point `q`. See [PennecLorenzi:2020](@cite) for details. """ function parallel_transport_to(M::CartanSchoutenMinusGroup, p, X, q) return inverse_translate_diff(M.manifold, q, p, X, LeftForwardAction()) @@ -158,7 +152,7 @@ end vector_transport_to(M::CartanSchoutenPlusGroup, p, X, q) Transport tangent vector `X` at point `p` on the group manifold `M` with the -[`CartanSchoutenPlus`](@ref) connection to point `q`. See [^Pennec2020] for details. +[`CartanSchoutenPlus`](@ref) connection to point `q`. See [PennecLorenzi:2020](@cite) for details. """ parallel_transport_to(M::CartanSchoutenPlusGroup, p, X, q) @@ -170,7 +164,7 @@ end parallel_transport_direction(M::CartanSchoutenZeroGroup, ::Identity, X, d) Transport tangent vector `X` at identity on the group manifold with the -[`CartanSchoutenZero`](@ref) connection in the direction `d`. See [^Pennec2020] for details. +[`CartanSchoutenZero`](@ref) connection in the direction `d`. See [PennecLorenzi:2020](@cite) for details. """ function parallel_transport_direction(M::CartanSchoutenZeroGroup, p::Identity, X, d) dexp_half = exp_lie(M.manifold, d / 2) diff --git a/src/groups/general_linear.jl b/src/groups/general_linear.jl index 9e69dbf2ff..844d3b429c 100644 --- a/src/groups/general_linear.jl +++ b/src/groups/general_linear.jl @@ -65,20 +65,7 @@ The exponential map is ```` where ``\operatorname{Exp}(⋅)`` denotes the matrix exponential, and ``⋅^\mathrm{H}`` is -the conjugate transpose. [^AndruchowLarotondaRechtVarela2014][^MartinNeff2016] - -[^AndruchowLarotondaRechtVarela2014]: - > Andruchow E., Larotonda G., Recht L., and Varela A.: - > “The left invariant metric in the general linear group”, - > Journal of Geometry and Physics 86, pp. 241-257, 2014. - > doi: [10.1016/j.geomphys.2014.08.009](https://doi.org/10.1016/j.geomphys.2014.08.009), - > arXiv: [1109.0520v1](https://arxiv.org/abs/1109.0520v1). -[^MartinNeff2016]: - > Martin, R. J. and Neff, P.: - > “Minimal geodesics on GL(n) for left-invariant, right-O(n)-invariant Riemannian metrics”, - > Journal of Geometric Mechanics 8(3), pp. 323-357, 2016. - > doi: [10.3934/jgm.2016010](https://doi.org/10.3934/jgm.2016010), - > arXiv: [1409.7849v2](https://arxiv.org/abs/1409.7849v2). +the conjugate transpose [AndruchowLarotondaRechtVarela:2014](@cite) [MartinNeff:2016](@cite). """ function exp(M::GeneralLinear, p, X) q = similar(p) diff --git a/src/groups/general_unitary_groups.jl b/src/groups/general_unitary_groups.jl index c56311b86c..2ceba9969f 100644 --- a/src/groups/general_unitary_groups.jl +++ b/src/groups/general_unitary_groups.jl @@ -62,18 +62,7 @@ exp_lie(::GeneralUnitaryMultiplicationGroup{2,ℝ}, X) exp_lie(G::SpecialOrthogonal{4}, X) Compute the group exponential map on the [`Orthogonal`](@ref)`(4)` or the [`SpecialOrthogonal`](@ref) group. -The algorithm used is a more numerically stable form of those proposed in [^Gallier2002], [^Andrica2013]. - -[^Gallier2002]: - > Gallier J.; Xu D.; Computing exponentials of skew-symmetric matrices - > and logarithms of orthogonal matrices. - > International Journal of Robotics and Automation (2002), 17(4), pp. 1-11. - > [pdf](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.35.3205). -[^Andrica2013]: - > Andrica D.; Rohan R.-A.; Computing the Rodrigues coefficients of the - > exponential map of the Lie groups of matrices. - > Balkan Journal of Geometry and Its Applications (2013), 18(2), pp. 1-2. - > [pdf](https://www.emis.de/journals/BJGA/v18n2/B18-2-an.pdf). +The algorithm used is a more numerically stable form of those proposed in [GallierXu:2002](@cite), [AndricaRohan:2013](@cite). """ exp_lie(::GeneralUnitaryMultiplicationGroup{4,ℝ}, X) diff --git a/src/groups/group.jl b/src/groups/group.jl index 5ef80603f1..b91587236b 100644 --- a/src/groups/group.jl +++ b/src/groups/group.jl @@ -115,7 +115,7 @@ abstract type ActionDirection end @doc raw""" LeftForwardAction() -Left action of a group on a manifold. For an action ``α: G × X → X`` it is characterized by +Left action of a group on a manifold. For an action ``α: G × X → X`` it is characterized by ```math α(g, α(h, x)) = α(gh, x) ``` @@ -126,7 +126,7 @@ struct LeftForwardAction <: ActionDirection end @doc raw""" LeftBackwardAction() -Left action of a group on a manifold. For an action ``α: X × G → X`` it is characterized by +Left action of a group on a manifold. For an action ``α: X × G → X`` it is characterized by ```math α(α(x, h), g) = α(x, gh) ``` @@ -141,7 +141,7 @@ const LeftAction = LeftForwardAction """ RightForwardAction() -Right action of a group on a manifold. For an action ``α: G × X → X`` it is characterized by +Right action of a group on a manifold. For an action ``α: G × X → X`` it is characterized by ```math α(g, α(h, x)) = α(hg, x) ``` @@ -154,7 +154,7 @@ struct RightForwardAction <: ActionDirection end """ RightBackwardAction() -Right action of a group on a manifold. For an action ``α: X × G → X`` it is characterized by +Right action of a group on a manifold. For an action ``α: X × G → X`` it is characterized by ```math α(α(x, h), g) = α(x, hg) ``` @@ -947,7 +947,7 @@ the group $\mathcal{G}$, and $𝔤$ is its Lie algebra, the group exponential is \exp : 𝔤 → \mathcal{G}, ```` such that for $t,s ∈ ℝ$, $γ(t) = \exp (t X)$ defines a one-parameter subgroup with the -following properties. Note that one-parameter subgroups are commutative (see [^Suhubi2013], +following properties. Note that one-parameter subgroups are commutative (see [Suhubi:2013](@cite), section 3.5), even if the Lie group itself is not commutative. ````math @@ -972,10 +972,6 @@ the Lie exponential is the numeric/matrix exponential. Since this function also depends on the group operation, make sure to implement the corresponding trait version `exp_lie(::TraitList{<:IsGroupManifold}, G, X)`. - -[^Suhubi2013]: - > E. Suhubi, Exterior Analysis: Using Applications of Differential Forms, 1st edition. - > Amsterdam: Academic Press, 2013. """ exp_lie(G::AbstractManifold, X) @trait_function exp_lie(M::AbstractDecoratorManifold, X) diff --git a/src/groups/heisenberg.jl b/src/groups/heisenberg.jl index 2c7cfa2608..e155b9337f 100644 --- a/src/groups/heisenberg.jl +++ b/src/groups/heisenberg.jl @@ -2,7 +2,7 @@ @doc raw""" HeisenbergGroup{n} <: AbstractDecoratorManifold{ℝ} -Heisenberg group `HeisenbergGroup(n)` is the group of ``(n+2) × (n+2)`` matrices [^BinzPods2008] +Heisenberg group `HeisenbergGroup(n)` is the group of ``(n+2) × (n+2)`` matrices [BinzPods:2008](@cite) ```math \begin{bmatrix} 1 & \mathbf{a} & c \\ @@ -15,10 +15,6 @@ where ``I_n`` is the ``n×n`` unit matrix, ``\mathbf{a}`` is a row vector of len The group operation is matrix multiplication. The left-invariant metric on the manifold is used. - -[^BinzPods2008]: - > E. Binz and S. Pods, The Geometry of Heisenberg Groups: With Applications in Signal - > Theory, Optics, Quantization, and Field Quantization. American Mathematical Soc., 2008. """ struct HeisenbergGroup{n} <: AbstractDecoratorManifold{ℝ} end diff --git a/src/groups/rotation_action.jl b/src/groups/rotation_action.jl index a601d3e568..4ae5fc77af 100644 --- a/src/groups/rotation_action.jl +++ b/src/groups/rotation_action.jl @@ -314,7 +314,7 @@ end Compute optimal alignment for the left [`ColumnwiseMultiplicationAction`](@ref), i.e. the group element ``O^{*}`` that, when it acts on `p`, returns the point closest to `q`. Details -of computation are described in Section 2.2.1 of [^Srivastava2016]. +of computation are described in Section 2.2.1 of [SrivastavaKlassen:2016](@cite). The formula reads ```math @@ -325,13 +325,6 @@ U K V^{\mathrm{T}} & \text{otherwise} ``` where ``U \Sigma V^{\mathrm{T}}`` is the SVD decomposition of ``p q^{\mathrm{T}}`` and ``K`` is the unit diagonal matrix with the last element on the diagonal replaced with -1. - -# References - -[^Srivastava2016]: - > A. Srivastava and E. P. Klassen, Functional and Shape Data Analysis. Springer New York, 2016. - > ISBN: 978-1-4939-4018-9. - > doi: [10.1007/978-1-4939-4020-2](https://doi.org/10.1007/978-1-4939-4020-2). """ function optimal_alignment(A::LeftColumnwiseMultiplicationAction, p, q) is_point(A.manifold, p, true) diff --git a/src/groups/special_euclidean.jl b/src/groups/special_euclidean.jl index 4733702056..5b563f0ff1 100644 --- a/src/groups/special_euclidean.jl +++ b/src/groups/special_euclidean.jl @@ -151,13 +151,10 @@ R & t \\ ```` This function embeds $\mathrm{SE}(n)$ in the general linear group $\mathrm{GL}(n+1)$. -It is an isometric embedding and group homomorphism [^RicoMartinez1988]. +It is an isometric embedding and group homomorphism [RicoMartinez:1988](@cite). See also [`screw_matrix`](@ref) for matrix representations of the Lie algebra. -[^RicoMartinez1988]: - > Rico Martinez, J. M., “Representations of the Euclidean group and its applications - > to the kinematics of spatial chains,” PhD Thesis, University of Florida, 1988. """ function affine_matrix(G::SpecialEuclidean{n}, p) where {n} pis = submanifold_components(G, p) diff --git a/src/groups/unitary.jl b/src/groups/unitary.jl index 3534e0a0b0..127c6850da 100644 --- a/src/groups/unitary.jl +++ b/src/groups/unitary.jl @@ -4,7 +4,7 @@ The group of unitary matrices ``\mathrm{U}(n, 𝔽)``, either complex (when 𝔽=ℂ) or quaternionic (when 𝔽=ℍ) -The group consists of all points ``p ∈ 𝔽^{n × n}`` where ``p^\mathrm{H}p = pp^\mathrm{H} = I``. +The group consists of all points ``p ∈ 𝔽^{n × n}`` where ``p^{\mathrm{H}}p = pp^{\mathrm{H}} = I``. The tangent spaces are if the form diff --git a/src/manifolds/CenteredMatrices.jl b/src/manifolds/CenteredMatrices.jl index 49e67a46dd..49c9342bd7 100644 --- a/src/manifolds/CenteredMatrices.jl +++ b/src/manifolds/CenteredMatrices.jl @@ -127,3 +127,16 @@ project!(::CenteredMatrices, Y, p, X) = (Y .= X .- mean(X, dims=1)) function Base.show(io::IO, ::CenteredMatrices{m,n,𝔽}) where {m,n,𝔽} return print(io, "CenteredMatrices($(m), $(n), $(𝔽))") end + +@doc raw""" + Y = Weingarten(M::CenteredMatrices, p, X, V) + Weingarten!(M::CenteredMatrices, Y, p, X, V) + +Compute the Weingarten map ``\mathcal W_p`` at `p` on the [`CenteredMatrices`](@ref) `M` with respect to the +tangent vector ``X \in T_p\mathcal M`` and the normal vector ``V \in N_p\mathcal M``. + +Since this a flat space by itself, the result is always the zero tangent vector. +""" +Weingarten(::CenteredMatrices, p, X, V) + +Weingarten!(::CenteredMatrices, Y, p, X, V) = fill!(Y, 0) diff --git a/src/manifolds/CholeskySpace.jl b/src/manifolds/CholeskySpace.jl index 5261c0f0d5..40f8cbd56f 100644 --- a/src/manifolds/CholeskySpace.jl +++ b/src/manifolds/CholeskySpace.jl @@ -3,17 +3,13 @@ The manifold of lower triangular matrices with positive diagonal and a metric based on the cholesky decomposition. The formulae for this manifold -are for example summarized in Table 1 of [^Lin2019]. +are for example summarized in Table 1 of [Lin:2019](@cite). # Constructor CholeskySpace(n) Generate the manifold of $n× n$ lower triangular matrices with positive diagonal. - -[^Lin2019]: - > Lin, Zenhua: _Riemannian Geometry of Symmetric Positive Definite Matrices via - > Cholesky Decomposition_, arXiv: [1908.09326](https://arxiv.org/abs/1908.09326). """ struct CholeskySpace{N} <: AbstractManifold{ℝ} end diff --git a/src/manifolds/Elliptope.jl b/src/manifolds/Elliptope.jl index feae2c820b..0f52dcce53 100644 --- a/src/manifolds/Elliptope.jl +++ b/src/manifolds/Elliptope.jl @@ -31,20 +31,13 @@ endowed with the [`Euclidean`](@ref) metric from the embedding, i.e. from the $ This manifold was for example -investigated in[^JourneeBachAbsilSepulchre2010]. +investigated in[JourneeBachAbsilSepulchre:2010](@cite). # Constructor Elliptope(n,k) generates the manifold $\mathcal E(n,k) \subset ℝ^{n × n}$. - -[^JourneeBachAbsilSepulchre2010]: - > Journée, M., Bach, F., Absil, P.-A., and Sepulchre, R.: - > “Low-Rank Optimization on the Cone of Positive Semidefinite Matrices”, - > SIAM Journal on Optimization (20)5, pp. 2327–2351, 2010. - > doi: [10.1137/080731359](https://doi.org/10.1137/080731359), - > arXiv: [0807.4423](http://arxiv.org/abs/0807.4423). """ struct Elliptope{N,K} <: AbstractDecoratorManifold{ℝ} end diff --git a/src/manifolds/EssentialManifold.jl b/src/manifolds/EssentialManifold.jl index f43961c010..a68d8dea8e 100644 --- a/src/manifolds/EssentialManifold.jl +++ b/src/manifolds/EssentialManifold.jl @@ -40,7 +40,7 @@ E = (R'_1)^T [T'_2 - T'_1]_{×} R'_2, where the poses of two cameras ``(R_i', T_i'), i=1,2``, are contained in the space of rigid body transformations $SE(3)$ and the operator $[⋅]_{×}\colon ℝ^3 \to \operatorname{SkewSym}(3)$ -denotes the matrix representation of the cross product operator. For more details see [^TronDaniilidis2017]. +denotes the matrix representation of the cross product operator. For more details see [TronDaniilidis:2017](@cite). # Constructor EssentialManifold(is_signed=true) @@ -48,12 +48,6 @@ denotes the matrix representation of the cross product operator. For more detail Generate the manifold of essential matrices, either the signed (`is_signed=true`) or unsigned (`is_signed=false`) variant. -[^TronDaniilidis2017]: - > Tron R.; Daniilidis K.; The Space of Essential Matrices as a Riemannian Quotient - > AbstractManifold. - > SIAM Journal on Imaging Sciences (2017), - > DOI: [10.1137/16M1091332](https://doi.org/10.1137/16M1091332), - > PDF: [https://www.cis.upenn.edu/~kostas/mypub.dir/tron17siam.pdf](https://www.cis.upenn.edu/~kostas/mypub.dir/tron17siam.pdf). """ struct EssentialManifold <: AbstractPowerManifold{ℝ,Rotations{3},NestedPowerRepresentation} is_signed::Bool @@ -110,7 +104,7 @@ $θ$ and $t_{\text{opt}}$ is the minimizer of the cost function ````math f(t) = f_1 + f_2, \quad f_i = \frac{1}{2} θ^2_i(t), \quad θ_i(t)=d(R_{p_i},R_z(t)R_{b_i}) \text{ for } i=1,2, ```` -where $d(⋅,⋅)$ is the distance function in $SO(3)$[^TronDaniilidis2017]. +where $d(⋅,⋅)$ is the distance function in $SO(3)$ [TronDaniilidis:2017](@cite). """ distance(M::EssentialManifold, p, q) = norm(M, p, log(M, p, q)) @@ -122,7 +116,7 @@ Compute the exponential map on the [`EssentialManifold`](@ref) from `p` into dir ````math \text{exp}_p(X) =\text{exp}_g( \tilde X), \quad g \in \text(SO)(3)^2, ```` -where $\tilde X$ is the horizontal lift of $X$[^TronDaniilidis2017]. +where $\tilde X$ is the horizontal lift of $X$[TronDaniilidis:2017](@cite). """ exp(::EssentialManifold, ::Any...) @@ -168,7 +162,7 @@ as \log_p (S_z(t_{\text{opt}})q) = [\text{Log}(R_{p_i}^T R_z(t_{\text{opt}})R_{b_i})]_{i=1,2}, ```` where $\text{Log}$ is the [`logarithm`](@ref log(::Rotations, ::Any...)) on $SO(3)$. For more -details see [^TronDaniilidis2017]. +details see [TronDaniilidis:2017](@cite). """ log(M::EssentialManifold, ::Any, ::Any) @@ -221,7 +215,7 @@ f(t) = f_1 + f_2, \quad f_i = \frac{1}{2} θ^2_i(t), \quad θ_i(t)=d(R_{p_i},R_z for the given values. This is done by finding the discontinuity points $t_{d_i}, i=1,2$ of its derivative and using Newton's method to minimize the function over the intervals $[t_{d_1},t_{d_2}]$ and $[t_{d_2},t_{d_1}+2π]$ separately. Then, the minimizer for which $f$ is minimal is chosen and given back together with the minimal value. -For more details see Algorithm 1 in [^TronDaniilidis2017]. +For more details see Algorithm 1 in [TronDaniilidis:2017](@cite). """ function dist_min_angle_pair(p, q) #compute rotations @@ -321,7 +315,7 @@ This function computes the point $t_{\text{di}}$ for which the first derivative f(t) = f_1 + f_2, \quad f_i = \frac{1}{2} θ^2_i(t), \quad θ_i(t)=d(R_{p_i},R_z(t)R_{b_i}) \text{ for } i=1,2, ```` does not exist. This is the case for $\sin(θ_i(t_{\text{di}})) = 0$. For more details see Proposition 9 -and its proof, as well as Lemma 1 in [^TronDaniilidis2017]. +and its proof, as well as Lemma 1 in [TronDaniilidis:2017](@cite). """ function dist_min_angle_pair_discontinuity_distance(q) c1 = q[1, 1] + q[2, 2] @@ -338,7 +332,7 @@ end @doc raw""" dist_min_angle_pair_compute_df_break(t_break, q) -This function computes the derivatives of each term $f_i, i=1,2,$ at discontinuity point `t_break`. For more details see [^TronDaniilidis2017]. +This function computes the derivatives of each term $f_i, i=1,2,$ at discontinuity point `t_break`. For more details see [TronDaniilidis:2017](@cite). """ function dist_min_angle_pair_compute_df_break(t_break, q) c = cos(t_break) @@ -357,7 +351,7 @@ This function computes the minimizer of the function ````math f(t) = f_1 + f_2, \quad f_i = \frac{1}{2} θ^2_i(t), \quad θ_i(t)=d(R_{p_i},R_z(t)R_{b_i}) \text{ for } i=1,2, ```` -in the interval $[$`t_low`, `t_high`$]$ using Newton's method. For more details see [^TronDaniilidis2017]. +in the interval $[$`t_low`, `t_high`$]$ using Newton's method. For more details see [TronDaniilidis:2017](@cite). """ function dist_min_angle_pair_df_newton(m1, Φ1, c1, m2, Φ2, c2, t_min, t_low, t_high) tol_dist = sqrt(eps(eltype(t_min))) @@ -400,7 +394,7 @@ end @doc raw""" manifold_dimension(M::EssentialManifold{is_signed, ℝ}) -Return the manifold dimension of the [`EssentialManifold`](@ref), which is `5`[^TronDaniilidis2017]. +Return the manifold dimension of the [`EssentialManifold`](@ref), which is `5`[TronDaniilidis:2017](@cite). """ function manifold_dimension(::EssentialManifold) return 5 @@ -483,7 +477,7 @@ Project `X` onto the vertical space $T_{\text{vp}}\text{SO}(3)^2$ with \text{vert\_proj}_p(X) = e_z^T(R_1 X_1 + R_2 X_2), ```` where $e_z$ is the third unit vector, $X_i ∈ T_{p}\text{SO}(3)$ for $i=1,2,$ and it holds $R_i = R_0 R'_i, i=1,2,$ where $R'_i$ is part of the -pose of camera $i$ $g_i = (R_i,T'_i) ∈ \text{SE}(3)$ and $R_0 ∈ \text{SO}(3)$ such that $R_0(T'_2-T'_1) = e_z$ [^TronDaniilidis2017]. +pose of camera $i$ $g_i = (R_i,T'_i) ∈ \text{SE}(3)$ and $R_0 ∈ \text{SO}(3)$ such that $R_0(T'_2-T'_1) = e_z$ [TronDaniilidis:2017](@cite). """ function vert_proj(M::EssentialManifold, p, X) return sum(vert_proj.(Ref(M.manifold), p, X)) diff --git a/src/manifolds/Euclidean.jl b/src/manifolds/Euclidean.jl index cd3b663aa1..34d5ac01ca 100644 --- a/src/manifolds/Euclidean.jl +++ b/src/manifolds/Euclidean.jl @@ -725,6 +725,19 @@ function volume_density(::Euclidean, p, X) return one(eltype(X)) end +@doc raw""" + Y = Weingarten(M::Euclidean, p, X, V) + Weingarten!(M::Euclidean, Y, p, X, V) + +Compute the Weingarten map ``\mathcal W_p`` at `p` on the [`Euclidean`](@ref) `M` with respect to the +tangent vector ``X \in T_p\mathcal M`` and the normal vector ``V \in N_p\mathcal M``. + +Since this a flat space by itself, the result is always the zero tangent vector. +""" +Weingarten(::Euclidean, p, X, V) + +Weingarten!(::Euclidean, Y, p, X, V) = fill!(Y, 0) + """ zero_vector(M::Euclidean, x) diff --git a/src/manifolds/FixedRankMatrices.jl b/src/manifolds/FixedRankMatrices.jl index ffcae9c3f1..4b19ff6ff3 100644 --- a/src/manifolds/FixedRankMatrices.jl +++ b/src/manifolds/FixedRankMatrices.jl @@ -3,7 +3,7 @@ The manifold of ``m × n`` real-valued or complex-valued matrices of fixed rank ``k``, i.e. ````math -\bigl\{ p ∈ 𝔽^{m × n}\ \big|\ \operatorname{rank}(p) = k \bigr\}, +\bigl\{ p ∈ 𝔽^{m × n}\ \big|\ \operatorname{rank}(p) = k\bigr\}, ```` where ``𝔽 ∈ \{ℝ,ℂ\}`` and the rank is the number of linearly independent columns of a matrix. @@ -29,18 +29,12 @@ T_p\mathcal M = \bigl\{ U_p M V_p^\mathrm{H} + U_X V_p^\mathrm{H} + U_p V_X^\mat where ``0_k`` is the ``k × k`` zero matrix. See [`UMVTVector`](@ref) for details. The (default) metric of this manifold is obtained by restricting the metric -on ``ℝ^{m × n}`` to the tangent bundle[^Vandereycken2013]. +on ``ℝ^{m × n}`` to the tangent bundle [Vandereycken:2013](@cite). # Constructor FixedRankMatrices(m, n, k[, field=ℝ]) Generate the manifold of `m`-by-`n` (`field`-valued) matrices of rank `k`. - -[^Vandereycken2013]: - > Bart Vandereycken: "Low-rank matrix completion by Riemannian Optimization, - > SIAM Journal on Optiomoization, 23(2), pp. 1214–1236, 2013. - > doi: [10.1137/110845768](https://doi.org/10.1137/110845768), - > arXiv: [1209.3834](https://arxiv.org/abs/1209.3834). """ struct FixedRankMatrices{M,N,K,𝔽} <: AbstractDecoratorManifold{𝔽} end function FixedRankMatrices(m::Int, n::Int, k::Int, field::AbstractNumbers=ℝ) @@ -356,12 +350,7 @@ get_embedding(::FixedRankMatrices{m,n,k,𝔽}) where {m,n,k,𝔽} = Euclidean(m, injectivity_radius(::FixedRankMatrices) Return the incjectivity radius of the manifold of [`FixedRankMatrices`](@ref), i.e. 0. -See [^HosseiniUschmajew2017]. - -[^HosseiniUschmajew2017]: - > S. Hosseini and A. Uschmajew, “A Riemannian Gradient Sampling Algorithm for Nonsmooth - > Optimization on Manifolds,” SIAM J. Optim., vol. 27, no. 1, pp. 173–189, Jan. 2017, - > doi: [10.1137/16M1069298](https://doi.org/10.1137/16M1069298). +See [HosseiniUschmajew:2017](@cite). """ function injectivity_radius(::FixedRankMatrices{m,n,k}) where {m,n,k} return 0.0 @@ -565,6 +554,28 @@ function retract_polar!( return q end +@doc raw""" + Y = riemannian_Hessian(M::FixedRankMatrices, p, G, H, X) + riemannian_Hessian!(M::FixedRankMatrices, Y, p, G, H, X) + +Compute the Riemannian Hessian ``\operatorname{Hess} f(p)[X]`` given the +Euclidean gradient ``∇ f(\tilde p)`` in `G` and the Euclidean Hessian ``∇^2 f(\tilde p)[\tilde X]`` in `H`, +where ``\tilde p, \tilde X`` are the representations of ``p,X`` in the embedding,. + +The Riemannian Hessian can be computed as stated in Remark 4.1 [Nguyen:2023](@cite) +or Section 2.3 [Vandereycken:2013](@cite), that B. Vandereycken adopted for [Manopt (Matlab)](https://www.manopt.org/reference/manopt/manifolds/fixedrank/fixedrankembeddedfactory.html). +""" +riemannian_Hessian(M::FixedRankMatrices, p, G, H, X) + +function riemannian_Hessian!(M::FixedRankMatrices, Y, p, G, H, X) + project!(M, Y, p, H) + T1 = (G * X.Vt) / Diagonal(p.S) + Y.U .+= T1 .- p.U * (p.U' * T1) + T2 = (G' * X.U) / Diagonal(p.S) + Y.Vt .+= T2 .- p.Vt' * (p.Vt * T2) + return Y +end + function Base.show(io::IO, ::FixedRankMatrices{M,N,K,𝔽}) where {M,N,K,𝔽} return print(io, "FixedRankMatrices($(M), $(N), $(K), $(𝔽))") end diff --git a/src/manifolds/Flag.jl b/src/manifolds/Flag.jl index ab6fb95f96..46ba6b6d13 100644 --- a/src/manifolds/Flag.jl +++ b/src/manifolds/Flag.jl @@ -47,7 +47,7 @@ end @doc raw""" Flag{N,d} <: AbstractDecoratorManifold{ℝ} -Flag manifold of ``d`` subspaces of ``ℝ^N``[^YeWongLim2022]. By default the manifold uses +Flag manifold of ``d`` subspaces of ``ℝ^N`` [YeWongLim:2021](@cite). By default the manifold uses the Stiefel coordinates representation, embedding it in the [`Stiefel`](@ref) manifold. The other available representation is an embedding in [`OrthogonalMatrices`](@ref). It can be utilized using [`OrthogonalPoint`](@ref) and [`OrthogonalTVector`](@ref) wrappers. @@ -64,12 +64,6 @@ Generate the manifold ``\operatorname{Flag}(n_1, n_2, ..., n_d; N)`` of subspace ``` where ``𝕍_i`` for ``i ∈ 1, 2, …, d`` are subspaces of ``ℝ^N`` of dimension ``\operatorname{dim} 𝕍_i = n_i``. - - -[^YeWongLim2022]: - > K. Ye, K. S.-W. Wong, and L.-H. Lim, “Optimization on flag manifolds,” Math. Program., - > vol. 194, no. 1, pp. 621–660, Jul. 2022, - > doi: [10.1007/s10107-021-01640-3](https://doi.org/10.1007/s10107-021-01640-3). """ struct Flag{N,dp1} <: AbstractDecoratorManifold{ℝ} subspace_dimensions::ZeroTuple{NTuple{dp1,Int}} diff --git a/src/manifolds/FlagStiefel.jl b/src/manifolds/FlagStiefel.jl index 9f0b91e0ec..40cd6aee59 100644 --- a/src/manifolds/FlagStiefel.jl +++ b/src/manifolds/FlagStiefel.jl @@ -120,7 +120,7 @@ end project(::Flag, p, X) Project vector `X` in the Euclidean embedding to the tangent space at point `p` on -[`Flag`](@ref) manifold. The formula reads[^YeWongLim2022]: +[`Flag`](@ref) manifold. The formula reads [YeWongLim:2021](@cite): ```math Y_i = X_i - (p_i p_i^{\mathrm{T}}) X_i + \sum_{j \neq i} p_j X_j^{\mathrm{T}} p_i diff --git a/src/manifolds/GeneralUnitaryMatrices.jl b/src/manifolds/GeneralUnitaryMatrices.jl index d2abb58a19..23c16724b7 100644 --- a/src/manifolds/GeneralUnitaryMatrices.jl +++ b/src/manifolds/GeneralUnitaryMatrices.jl @@ -201,20 +201,7 @@ exp_p(X) = p\mathrm{e}^X For different sizes, like ``n=2,3,4`` there is specialised implementations The algorithm used is a more numerically stable form of those proposed in -[^Gallier2002] and [^Andrica2013]. - -[^Gallier2002]: - > Gallier J.; Xu D.; Computing exponentials of skew-symmetric matrices - > and logarithms of orthogonal matrices. - > International Journal of Robotics and Automation (2002), 17(4), pp. 1-11. - > [pdf](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.35.3205). - -[^Andrica2013]: - > Andrica D.; Rohan R.-A.; Computing the Rodrigues coefficients of the - > exponential map of the Lie groups of matrices. - > Balkan Journal of Geometry and Its Applications (2013), 18(2), pp. 1-2. - > [pdf](https://www.emis.de/journals/BJGA/v18n2/B18-2-an.pdf). - +[GallierXu:2002](@cite) and [AndricaRohan:2013](@cite). """ exp(::GeneralUnitaryMatrices, p, X) @@ -683,7 +670,7 @@ manifold_dimension(::GeneralUnitaryMatrices{n,ℂ,DeterminantOneMatrices}) where manifold_volume(::GeneralUnitaryMatrices{n,ℝ,AbsoluteDeterminantOneMatrices}) where {n} Volume of the manifold of real orthogonal matrices of absolute determinant one. The -formula reads [^BoyaSudarshanTilma2003]: +formula reads [BoyaSudarshanTilma:2003](@cite): ```math \begin{cases} @@ -691,11 +678,6 @@ formula reads [^BoyaSudarshanTilma2003]: \frac{2^{k+1}(2\pi)^{k(k+1)}}{\prod_{s=1}^{k-1} (2s+1)!} & \text{ if } n = 2k+1 \end{cases} ``` - -[^BoyaSudarshanTilma2003]: - > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports - > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003, - > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1) """ function manifold_volume( ::GeneralUnitaryMatrices{n,ℝ,AbsoluteDeterminantOneMatrices}, @@ -706,7 +688,7 @@ end manifold_volume(::GeneralUnitaryMatrices{n,ℝ,DeterminantOneMatrices}) where {n} Volume of the manifold of real orthogonal matrices of determinant one. The -formula reads [^BoyaSudarshanTilma2003]: +formula reads [BoyaSudarshanTilma:2003](@cite): ```math \begin{cases} @@ -718,11 +700,6 @@ formula reads [^BoyaSudarshanTilma2003]: It differs from the paper by a factor of `sqrt(2)` due to a different choice of normalization. - -[^BoyaSudarshanTilma2003]: - > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports - > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003, - > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1) """ function manifold_volume(::GeneralUnitaryMatrices{n,ℝ,DeterminantOneMatrices}) where {n} vol = 1.0 @@ -748,16 +725,11 @@ end manifold_volume(::GeneralUnitaryMatrices{n,ℂ,AbsoluteDeterminantOneMatrices}) where {n} Volume of the manifold of complex general unitary matrices of absolute determinant one. The -formula reads [^BoyaSudarshanTilma2003]: +formula reads [BoyaSudarshanTilma:2003](@cite) ```math \sqrt{n 2^{n+1}} π^{n(n+1)/2} \prod_{k=1}^{n-1}\frac{1}{k!} ``` - -[^BoyaSudarshanTilma2003]: - > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports - > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003, - > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1) """ function manifold_volume( ::GeneralUnitaryMatrices{n,ℂ,AbsoluteDeterminantOneMatrices}, @@ -774,16 +746,11 @@ end manifold_volume(::GeneralUnitaryMatrices{n,ℂ,DeterminantOneMatrices}) where {n} Volume of the manifold of complex general unitary matrices of determinant one. The formula -reads [^BoyaSudarshanTilma2003]: +reads [BoyaSudarshanTilma:2003](@cite) ```math \sqrt{n 2^{n-1}} π^{(n-1)(n+2)/2} \prod_{k=1}^{n-1}\frac{1}{k!} ``` - -[^BoyaSudarshanTilma2003]: - > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports - > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003, - > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1) """ function manifold_volume(::GeneralUnitaryMatrices{n,ℂ,DeterminantOneMatrices}) where {n} vol = sqrt(n * 2^(n - 1)) * π^(((n - 1) * (n + 2)) // 2) @@ -907,12 +874,7 @@ end riemann_tensor(::GeneralUnitaryMatrices, p, X, Y, Z) Compute the value of Riemann tensor on the [`GeneralUnitaryMatrices`](@ref) manifold. -The formula reads[^Rentmeesters2011] ``R(X,Y)Z=\frac{1}{4}[Z, [X, Y]]``. - -[^Rentmeesters2011]: - > Q. Rentmeesters, “A gradient method for geodesic data fitting on some symmetric - > Riemannian manifolds,” in 2011 50th IEEE Conference on Decision and Control and - > European Control Conference, Dec. 2011, pp. 7141–7146. doi: [10.1109/CDC.2011.6161280](https://doi.org/10.1109/CDC.2011.6161280). +The formula reads [Rentmeesters:2011](@cite) ``R(X,Y)Z=\frac{1}{4}[Z, [X, Y]]``. """ riemann_tensor(::GeneralUnitaryMatrices, p, X, Y, Z) @@ -927,18 +889,8 @@ end Compute volume density function of a sphere, i.e. determinant of the differential of exponential map `exp(M, p, X)`. It is derived from Eq. (4.1) and Corollary 4.4 -in [^ChevallierLiLuDunson2022]. See also Theorem 4.1 in [^FalorsideHaanDavidsonForré2019], +in [ChevallierLiLuDunson:2022](@ref). See also Theorem 4.1 in [FalorsideHaanDavidsonForre:2019](@cite), (note that it uses a different convention). - -[^ChevallierLiLuDunson2022]: - > E. Chevallier, D. Li, Y. Lu, and D. B. Dunson, “Exponential-wrapped distributions on - > symmetric spaces.” arXiv, Oct. 09, 2022. - > doi: [10.48550/arXiv.2009.01983](https://doi.org/10.48550/arXiv.2009.01983). - -[^FalorsideHaanDavidsonForré2019]: - > L. Falorsi, P. de Haan, T. R. Davidson, and P. Forré, “Reparameterizing Distributions - > on Lie Groups,” arXiv:1903.02958 [cs, math, stat], Mar. 2019 - > doi: [10.48550/arXiv.1903.02958](https://doi.org/10.48550/arXiv.1903.02958) """ function volume_density(M::GeneralUnitaryMatrices{n,ℝ}, p, X) where {n} dens = one(eltype(X)) @@ -964,7 +916,8 @@ end @doc raw""" volume_density(M::GeneralUnitaryMatrices{3,ℝ}, p, X) -Compute the volume density on O(3)/SO(3). The formula reads [^FalorsideHaanDavidsonForré2019]: +Compute the volume density on O(3)/SO(3). The formula reads [FalorsideHaanDavidsonForre:2019](@cite) + ```math \frac{1-1\cos(\sqrt{2}\lVert X \rVert)}{\lVert X \rVert^2}. ``` diff --git a/src/manifolds/Grassmann.jl b/src/manifolds/Grassmann.jl index 0227ca1f30..85a176f8f6 100644 --- a/src/manifolds/Grassmann.jl +++ b/src/manifolds/Grassmann.jl @@ -61,7 +61,7 @@ see also [`ProjectorPoint`](@ref) and [`ProjectorTVector`](@ref). The manifold is named after [Hermann G. Graßmann](https://en.wikipedia.org/wiki/Hermann_Grassmann) (1809-1877). -A good overview can be found in[^BendokatZimmermannAbsil2020]. +A good overview can be found in[BendokatZimmermannAbsil:2020](@cite). # Constructor @@ -69,11 +69,6 @@ A good overview can be found in[^BendokatZimmermannAbsil2020]. Generate the Grassmann manifold $\operatorname{Gr}(n,k)$, where the real-valued case `field = ℝ` is the default. - -[^BendokatZimmermannAbsil2020]: - > T. Bendokat, R. Zimmermann, and P. -A. Absil: - > _A Grassmann Manifold Handbook: Basic Geometry and Computational Aspects_, - > arXiv preprint [2011.13699](https://arxiv.org/abs/2011.13699), 2020. """ struct Grassmann{n,k,𝔽} <: AbstractDecoratorManifold{𝔽} end diff --git a/src/manifolds/GrassmannProjector.jl b/src/manifolds/GrassmannProjector.jl index ea7571a11c..2d5984a89f 100644 --- a/src/manifolds/GrassmannProjector.jl +++ b/src/manifolds/GrassmannProjector.jl @@ -201,7 +201,7 @@ Compute the exponential map on the [`Grassmann`](@ref) as ``` where ``\operatorname{Exp}`` denotes the matrix exponential and ``[A,B] = AB-BA`` denotes the matrix commutator. -For details, see Proposition 3.2 in [^BendokatZimmermannAbsil2020]. +For details, see Proposition 3.2 in [BendokatZimmermannAbsil:2020](@cite). """ exp(M::Grassmann, p::ProjectorPoint, X::ProjectorTVector) @@ -235,7 +235,7 @@ horizontal_lift!(::Stiefel, Y, q, X::ProjectorTVector) = copyto!(Y, X.value * q) ) Compute the parallel transport of `X` from the tangent space at `p` into direction `d`, -i.e. to ``q=\exp_pd``. The formula is given in Proposition 3.5 of [^BendokatZimmermannAbsil2020] as +i.e. to ``q=\exp_pd``. The formula is given in Proposition 3.5 of [BendokatZimmermannAbsil:2020](@cite) as ```math \mathcal{P}_{q ← p}(X) = \operatorname{Exp}([d,p])X\operatorname{Exp}(-[d,p]), diff --git a/src/manifolds/GrassmannStiefel.jl b/src/manifolds/GrassmannStiefel.jl index 5019a602f8..bc67a36f94 100644 --- a/src/manifolds/GrassmannStiefel.jl +++ b/src/manifolds/GrassmannStiefel.jl @@ -294,17 +294,37 @@ function retract_qr!(::Grassmann{N,K}, q, p, X, t::Number) where {N,K} return q end +@doc raw""" + riemannian_Hessian(M::Grassmann, p, G, H, X) + +The Riemannian Hessian can be computed by adopting Eq. (6.6) [Nguyen:2023](@cite), +where we use for the [`EuclideanMetric`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.EuclideanMetric) ``α_0=α_1=1`` in their formula. +Let ``\nabla f(p)`` denote the Euclidean gradient `G`, +``\nabla^2 f(p)[X]`` the Euclidean Hessian `H`. Then the formula reads + +```math + \operatorname{Hess}f(p)[X] + = + \operatorname{proj}_{T_p\mathcal M}\Bigl( + ∇^2f(p)[X] - X p^{\mathrm{H}}∇f(p) + \Bigr). +``` + +Compared to Eq. (5.6) also the metric conversion simplifies to the identity. +""" +riemannian_Hessian(M::Grassmann, p, G, H, X) + +function riemannian_Hessian!(M::Grassmann, Y, p, G, H, X) + project!(M, Y, p, H - X * p' * G) + return Y +end + @doc raw""" riemann_tensor(::Grassmann{n,k,ℝ}, p, X, Y, Z) where {n,k} Compute the value of Riemann tensor on the real [`Grassmann`](@ref) manifold. -The formula reads[^Rentmeesters2011] +The formula reads [Rentmeesters:2011](@cite) ``R(X,Y)Z = (XY^\mathrm{T} - YX^\mathrm{T})Z + Z(Y^\mathrm{T}X - X^\mathrm{T}Y)``. - -[^Rentmeesters2011]: - > Q. Rentmeesters, “A gradient method for geodesic data fitting on some symmetric - > Riemannian manifolds,” in 2011 50th IEEE Conference on Decision and Control and - > European Control Conference, Dec. 2011, pp. 7141–7146. doi: [10.1109/CDC.2011.6161280](https://doi.org/10.1109/CDC.2011.6161280). """ riemann_tensor(::Grassmann{n,k,ℝ}, p, X, Y, Z) where {n,k} @@ -330,12 +350,8 @@ Uniform distribution on given (real-valued) [`Grassmann`](@ref) `M`. Specifically, this is the normalized Haar measure on `M`. Generated points will be of similar type as `p`. -The implementation is based on Section 2.5.1 in [^Chikuse2003]; -see also Theorem 2.2.2(iii) in [^Chikuse2003]. - -[^Chikuse2003]: - > Y. Chikuse: "Statistics on Special Manifolds", Springer New York, 2003, - > doi: [10.1007/978-0-387-21540-2](https://doi.org/10.1007/978-0-387-21540-2). +The implementation is based on Section 2.5.1 in [Chikuse:2003](@cite); +see also Theorem 2.2.2(iii) in [Chikuse:2003](@cite). """ function uniform_distribution(M::Grassmann{n,k,ℝ}, p) where {n,k} μ = Distributions.Zeros(n, k) diff --git a/src/manifolds/HyperbolicHyperboloid.jl b/src/manifolds/HyperbolicHyperboloid.jl index 681fc561db..4bd14068e7 100644 --- a/src/manifolds/HyperbolicHyperboloid.jl +++ b/src/manifolds/HyperbolicHyperboloid.jl @@ -433,17 +433,41 @@ function parallel_transport_to!(::Hyperbolic, Y, p, X, q) ) end +@doc raw""" + Y = riemannian_Hessian(M::Hyperbolic, p, G, H, X) + riemannian_Hessian!(M::Hyperbolic, Y, p, G, H, X) + +Compute the Riemannian Hessian ``\operatorname{Hess} f(p)[X]`` given the +Euclidean gradient ``∇ f(\tilde p)`` in `G` and the Euclidean Hessian ``∇^2 f(\tilde p)[\tilde X]`` in `H`, +where ``\tilde p, \tilde X`` are the representations of ``p,X`` in the embedding,. + +Let ``\mathbf{g} = \mathbf{g}^{-1} = \operatorname{diag}(1,...,1,-1)``. +Then using Remark 4.1 [Nguyen:2023](@cite) the formula reads + +```math +\operatorname{Hess}f(p)[X] += +\operatorname{proj}_{T_p\mathcal M}\bigl( + \mathbf{g}^{-1}\nabla^2f(p)[X] + X⟨p,\mathbf{g}^{-1}∇f(p)⟩_p +\bigr). +``` +""" +riemannian_Hessian(M::Hyperbolic, p, G, H, X) + +function riemannian_Hessian!(M::Hyperbolic, Y, p, G, H, X) + g = copy(G) + g[end] *= -1 # = g^{-1}G + h = copy(H) + H[end] *= -1 # = g^{-1}H + project!(M, Y, p, h .+ dot(p, g) .* X) + return Y +end @doc raw""" volume_density(M::Hyperbolic, p, X) Compute volume density function of the hyperbolic manifold. The formula reads ``(\sinh(\lVert X\rVert)/\lVert X\rVert)^(n-1)`` where `n` is the dimension of `M`. -It is derived from Eq. (4.1) in [^ChevallierLiLuDunson2022]. - -[^ChevallierLiLuDunson2022]: - > E. Chevallier, D. Li, Y. Lu, and D. B. Dunson, “Exponential-wrapped distributions on - > symmetric spaces.” arXiv, Oct. 09, 2022. - > doi: [10.48550/arXiv.2009.01983](https://doi.org/10.48550/arXiv.2009.01983). +It is derived from Eq. (4.1) in[ChevallierLiLuDunson:2022](@cite). """ function volume_density(M::Hyperbolic, p, X) Xnorm = norm(X) diff --git a/src/manifolds/KendallsPreShapeSpace.jl b/src/manifolds/KendallsPreShapeSpace.jl index adb937f14c..f069a493af 100644 --- a/src/manifolds/KendallsPreShapeSpace.jl +++ b/src/manifolds/KendallsPreShapeSpace.jl @@ -4,7 +4,7 @@ Kendall's pre-shape space of ``k`` landmarks in ``ℝ^n`` represented by n×k matrices. In each row the sum of elements of a matrix is equal to 0. The Frobenius norm of the matrix -is equal to 1 [^Kendall1984][^Kendall1989]. +is equal to 1 [Kendall:1984](@cite)[Kendall:1989](@cite). The space can be interpreted as tuples of ``k`` points in ``ℝ^n`` up to simultaneous translation and scaling of all points, so this can be thought of as a quotient manifold. @@ -88,7 +88,7 @@ manifold_dimension(::KendallsPreShapeSpace{n,k}) where {n,k} = n * (k - 1) - 1 Project point `p` from the embedding to [`KendallsPreShapeSpace`](@ref) by selecting the right element from the orthogonal section representing the quotient manifold `M`. -See Section 3.7 of [^Srivastava2016] for details. +See Section 3.7 of [SrivastavaKlassen:2016](@cite) for details. The method computes the mean of the landmarks and moves them to make their mean zero; afterwards the Frobenius norm of the landmarks (as a matrix) is normalised to fix the scaling. @@ -106,14 +106,7 @@ end Project tangent vector `X` at point `p` from the embedding to [`KendallsPreShapeSpace`](@ref) by selecting the right element from the tangent space to orthogonal section representing the -quotient manifold `M`. See Section 3.7 of [^Srivastava2016] for details. - -# References - -[^Srivastava2016]: - > A. Srivastava and E. P. Klassen, Functional and Shape Data Analysis. Springer New York, 2016. - > ISBN: 978-1-4939-4018-9. - > doi: [10.1007/978-1-4939-4020-2](https://doi.org/10.1007/978-1-4939-4020-2). +quotient manifold `M`. See Section 3.7 of [SrivastavaKlassen:2016](@cite) for details. """ project(::KendallsPreShapeSpace, p, X) diff --git a/src/manifolds/KendallsShapeSpace.jl b/src/manifolds/KendallsShapeSpace.jl index ba9e51c2b5..0bcb29db52 100644 --- a/src/manifolds/KendallsShapeSpace.jl +++ b/src/manifolds/KendallsShapeSpace.jl @@ -6,7 +6,7 @@ Kendall's shape space, defined as quotient of a [`KendallsPreShapeSpace`](@ref) (represented by n×k matrices) by the action [`ColumnwiseMultiplicationAction`](@ref). The space can be interpreted as tuples of ``k`` points in ``ℝ^n`` up to simultaneous -translation and scaling and rotation of all points [^Kendall1984][^Kendall1989]. +translation and scaling and rotation of all points [Kendall:1984](@cite)[Kendall:1989](@cite). This manifold possesses the [`IsQuotientManifold`](@ref) trait. @@ -15,15 +15,6 @@ This manifold possesses the [`IsQuotientManifold`](@ref) trait. KendallsShapeSpace(n::Int, k::Int) # References - -[^Kendall1989]: - > D. G. Kendall, “A Survey of the Statistical Theory of Shape,” Statist. Sci., vol. 4, - > no. 2, pp. 87–99, May 1989 - > doi: [10.1214/ss/1177012582](https://doi.org/10.1214/ss/1177012582). -[^Kendall1984]: - > D. G. Kendall, “Shape Manifolds, Procrustean Metrics, and Complex Projective Spaces,” - > Bull. London Math. Soc., vol. 16, no. 2, pp. 81–121, Mar. 1984 - > doi: [10.1112/blms/16.2.81](https://doi.org/10.1112/blms/16.2.81). """ struct KendallsShapeSpace{n,k} <: AbstractDecoratorManifold{ℝ} end @@ -56,12 +47,8 @@ end exp(M::KendallsShapeSpace, p, X) Compute the exponential map on [`KendallsShapeSpace`](@ref) `M`. -See [^Guigui2021] for discussion about its computation. +See [GuiguiMaignantTroouvePennec:2021](@cite) for discussion about its computation. -[^Guigui2021]: - > N. Guigui, E. Maignant, A. Trouvé, and X. Pennec, “Parallel Transport on Kendall Shape - > Spaces,” in Geometric Science of Information, Cham, 2021, pp. 103–110. - > doi: [10.1007/978-3-030-80209-7_12](https://doi.org/10.1007/978-3-030-80209-7_12). """ exp(M::KendallsShapeSpace, p, X) @@ -87,7 +74,7 @@ end horizontal_component(::KendallsShapeSpace, p, X) Compute the horizontal component of tangent vector `X` at `p` on [`KendallsShapeSpace`](@ref) -`M`. See [^Guigui2021], Section 2.3 for details. +`M`. See [GuiguiMaignantTroouvePennec:2021](@cite), Section 2.3 for details. """ horizontal_component(::KendallsShapeSpace, p, X) @@ -140,7 +127,7 @@ end Return the dimension of the [`KendallsShapeSpace`](@ref) manifold `M`. The dimension is given by ``n(k - 1) - 1 - n(n - 1)/2`` in the typical case where ``k \geq n+1``, and ``(k + 1)(k - 2) / 2`` otherwise, unless ``k`` is equal to 1, in which case the dimension -is 0. See [^Kendall1984] for a discussion of the over-dimensioned case. +is 0. See [Kendall:1984](@cite) for a discussion of the over-dimensioned case. """ function manifold_dimension(::KendallsShapeSpace{n,k}) where {n,k} if k < n + 1 # over-dimensioned case diff --git a/src/manifolds/MetricManifold.jl b/src/manifolds/MetricManifold.jl index 4b8311d029..e83539a050 100644 --- a/src/manifolds/MetricManifold.jl +++ b/src/manifolds/MetricManifold.jl @@ -754,6 +754,26 @@ function vector_transport_to!( return vector_transport_to!(M.manifold, Y, p, X, q, m) end +function Weingarten( + ::TraitList{IsDefaultMetric{G}}, + M::MetricManifold{𝔽,TM,G}, + p, + X, + V, +) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold} + return Weingarten(M.manifold, p, X, V) +end +function Weingarten!( + ::TraitList{IsDefaultMetric{G}}, + M::MetricManifold{𝔽,TM,G}, + Y, + p, + X, + V, +) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold} + return Weingarten!(M.manifold, Y, p, X, V) +end + zero_vector(M::MetricManifold, p) = zero_vector(M.manifold, p) zero_vector!(M::MetricManifold, X, p) = zero_vector!(M.manifold, X, p) @@ -804,6 +824,10 @@ for mf in [ ricci_curvature, ricci_tensor, riemann_tensor, + riemannian_gradient, + riemannian_gradient!, + riemannian_Hessian, + riemannian_Hessian!, sharp!, vector_transport_along, vector_transport_along!, @@ -811,6 +835,8 @@ for mf in [ vector_transport_direction!, vector_transport_to, vector_transport_to!, + Weingarten, + Weingarten!, ] @eval is_metric_function(::typeof($mf)) = true end diff --git a/src/manifolds/MultinomialDoublyStochastic.jl b/src/manifolds/MultinomialDoublyStochastic.jl index 7876c849e2..9e02dd6e88 100644 --- a/src/manifolds/MultinomialDoublyStochastic.jl +++ b/src/manifolds/MultinomialDoublyStochastic.jl @@ -37,20 +37,13 @@ X\mathbf{1}_n = X^{\mathrm{T}}\mathbf{1}_n = \mathbf{0}_n where $\mathbf{0}_n$ is the vector of length $n$ containing zeros. -More details can be found in Section III[^DouikHassibi2019]. +More details can be found in Section III [DouikHassibi:2019](@cite). # Constructor MultinomialDoubleStochastic(n) Generate the manifold of matrices $\mathbb R^{n×n}$ that are doubly stochastic and symmetric. - -[^DouikHassibi2019]: - > A. Douik, B. Hassibi: - > AbstractManifold Optimization Over the Set of Doubly Stochastic Matrices: A Second-Order Geometry, - > IEEE Transactions on Signal Processing 67(22), pp. 5761–5774, 2019. - > doi: [10.1109/tsp.2019.2946024](http://doi.org/10.1109/tsp.2019.2946024), - > arXiv: [1802.02628](https://arxiv.org/abs/1802.02628). """ struct MultinomialDoubleStochastic{N} <: AbstractMultinomialDoublyStochastic{N} end diff --git a/src/manifolds/MultinomialSymmetric.jl b/src/manifolds/MultinomialSymmetric.jl index 68f854534d..3ee993c120 100644 --- a/src/manifolds/MultinomialSymmetric.jl +++ b/src/manifolds/MultinomialSymmetric.jl @@ -31,20 +31,13 @@ X\mathbf{1}_n = \mathbf{0}_n where $\mathbf{0}_n$ is the vector of length $n$ containing zeros. -More details can be found in Section IV[^DouikHassibi2019]. +More details can be found in Section IV [DouikHassibi:2019](@cite). # Constructor MultinomialSymmetric(n) Generate the manifold of matrices $\mathbb R^{n×n}$ that are doubly stochastic and symmetric. - -[^DouikHassibi2019]: - > A. Douik, B. Hassibi: - > AbstractManifold Optimization Over the Set of Doubly Stochastic Matrices: A Second-Order Geometry, - > IEEE Transactions on Signal Processing 67(22), pp. 5761–5774, 2019. - > doi: [10.1109/tsp.2019.2946024](http://doi.org/10.1109/tsp.2019.2946024), - > arXiv: [1802.02628](https://arxiv.org/abs/1802.02628). """ struct MultinomialSymmetric{N} <: AbstractMultinomialDoublyStochastic{N} end diff --git a/src/manifolds/PositiveNumbers.jl b/src/manifolds/PositiveNumbers.jl index 6196e2b4ec..da9ceadc73 100644 --- a/src/manifolds/PositiveNumbers.jl +++ b/src/manifolds/PositiveNumbers.jl @@ -309,6 +309,20 @@ function Random.rand!( return pX end +@doc raw""" + riemannian_Hessian(M::SymmetricPositiveDefinite, p, G, H, X) + +The Riemannian Hessian can be computed as stated in Eq. (7.3) [Nguyen:2023](@cite). +Let ``\nabla f(p)`` denote the Euclidean gradient `G`, +``\nabla^2 f(p)[X]`` the Euclidean Hessian `H`. +Then the formula reads + +```math + \operatorname{Hess}f(p)[X] = p\bigl(∇^2 f(p)[X]\bigr)p + X\bigl(∇f(p)\bigr)p +``` +""" +riemannian_Hessian(::PositiveNumbers, p, G, H, X) = p * H * p + X * G * p + function vector_transport_direction( M::PositiveNumbers, p, diff --git a/src/manifolds/PowerManifold.jl b/src/manifolds/PowerManifold.jl index 43cf98aef4..56aaa46574 100644 --- a/src/manifolds/PowerManifold.jl +++ b/src/manifolds/PowerManifold.jl @@ -180,7 +180,7 @@ end @doc raw""" manifold_volume(M::PowerManifold) -Return the manifold volume of an [`PowerManifold`](@ref) `M`. +Return the manifold volume of an [`PowerManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.PowerManifold) `M`. """ function manifold_volume(M::PowerManifold{𝔽,<:AbstractManifold,TSize}) where {𝔽,TSize} return manifold_volume(M.manifold)^prod(size_to_tuple(TSize)) @@ -266,6 +266,33 @@ function representation_size(M::PowerManifold{𝔽,<:AbstractManifold,TSize}) wh return (representation_size(M.manifold)..., size_to_tuple(TSize)...) end +@doc raw""" + Y = riemannian_Hessian(M::AbstractPowerManifold, p, G, H, X) + riemannian_Hessian!(M::AbstractPowerManifold, Y, p, G, H, X) + +Compute the Riemannian Hessian ``\operatorname{Hess} f(p)[X]`` given the +Euclidean gradient ``∇ f(\tilde p)`` in `G` and the Euclidean Hessian ``∇^2 f(\tilde p)[\tilde X]`` in `H`, +where ``\tilde p, \tilde X`` are the representations of ``p,X`` in the embedding,. + +On an abstract power manifold, this decouples and can be computed elementwise. +""" +riemannian_Hessian(M::AbstractPowerManifold, p, G, H, X) + +function riemannian_Hessian!(M::AbstractPowerManifold, Y, p, G, H, X) + rep_size = representation_size(M.manifold) + for i in get_iterator(M) + riemannian_Hessian!( + M.manifold, + _write(M, rep_size, Y, i), + _read(M, rep_size, p, i), + _read(M, rep_size, G, i), + _read(M, rep_size, H, i), + _read(M, rep_size, X, i), + ) + end + return Y +end + @doc raw""" sharp(M::AbstractPowerManifold, p, ξ::RieszRepresenterCotangentVector) @@ -306,7 +333,7 @@ end @doc raw""" volume_density(M::PowerManifold, p, X) -Return volume density on the [`PowerManifold`](@ref) `M`, i.e. product of constituent +Return volume density on the [`PowerManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.PowerManifold) `M`, i.e. product of constituent volume densities. """ function volume_density(M::PowerManifold, p, X) @@ -320,6 +347,29 @@ function volume_density(M::PowerManifold, p, X) return density end +@doc raw""" + Y = Weingarten(M::AbstractPowerManifold, p, X, V) + Weingarten!(M::AbstractPowerManifold, Y, p, X, V) + +Since the metric decouples, also the computation of the Weingarten map +``\mathcal W_p`` can be computed elementwise on the single elements of the [`PowerManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds/#sec-power-manifold) `M`. +""" +Weingarten(::AbstractPowerManifold, p, X, V) + +function Weingarten!(M::AbstractPowerManifold, Y, p, X, V) + rep_size = representation_size(M.manifold) + for i in get_iterator(M) + Weingarten!( + M.manifold, + _write(M, rep_size, Y, i), + _read(M, rep_size, p, i), + _read(M, rep_size, X, i), + _read(M, rep_size, V, i), + ) + end + return Y +end + @inline function _write( ::PowerManifoldMultidimensional, rep_size::Tuple, diff --git a/src/manifolds/ProbabilitySimplex.jl b/src/manifolds/ProbabilitySimplex.jl index 5a1d2c5ec6..83fea8ec58 100644 --- a/src/manifolds/ProbabilitySimplex.jl +++ b/src/manifolds/ProbabilitySimplex.jl @@ -28,17 +28,11 @@ simplex in the $n$-sphere of radius 2. The corresponding diffeomorphism $\varphi: \mathbb Δ^n → \mathcal N$, where $\mathcal N \subset 2𝕊^n$ is given by $\varphi(p) = 2\sqrt{p}$. -This implementation follows the notation in [^ÅströmPetraSchmitzerSchnörr2017]. +This implementation follows the notation in [AastroemPetraSchmitzerSchnoerr:2017](@cite). # Constructor ProbabilitySimplex(n::Int; boundary::Symbol=:open) - -[^ÅströmPetraSchmitzerSchnörr2017]: - > F. Åström, S. Petra, B. Schmitzer, C. Schnörr: “Image Labeling by Assignment”, - > Journal of Mathematical Imaging and Vision, 58(2), pp. 221–238, 2017. - > doi: [10.1007/s10851-016-0702-4](https://doi.org/10.1007/s10851-016-0702-4) - > arxiv: [1603.05285](https://arxiv.org/abs/1603.05285). """ struct ProbabilitySimplex{n,boundary} <: AbstractDecoratorManifold{ℝ} end @@ -74,7 +68,7 @@ the representer of a linear functional on the tangent space is adapted as ``Z = The first part “compensates” for the divsion by ``p`` in the Riemannian metric on the [`ProbabilitySimplex`](@ref) and the second part performs appropriate projection to keep the vector tangent. -For details see Proposition 2.3 in [^ÅströmPetraSchmitzerSchnörr2017]. +For details see Proposition 2.3 in [AastroemPetraSchmitzerSchnoerr:2017](@cite). """ change_representer(::ProbabilitySimplex, ::EuclideanMetric, ::Any, ::Any) @@ -244,13 +238,7 @@ Compute the inner product of two tangent vectors `X`, `Y` from the tangent space g_p(X,Y) = \sum_{i=1}^{n+1}\frac{X_iY_i}{p_i} ```` When `M` includes boundary, we can just skip coordinates where ``p_i`` is equal to 0, see -Proposition 2.1 in [^AyJostLeSchwachhöfer2017]. - -[^AyJostLeSchwachhöfer2017]: - > N. Ay, J. Jost, H. V. Le, and L. Schwachhöfer, Information Geometry. in Ergebnisse der - > Mathematik und ihrer Grenzgebiete. 3. Folge / A Series of Modern Surveys in - > Mathematics. Springer International Publishing, 2017. - > doi: [10.1007/978-3-319-56478-4](https://doi.org/10.1007/978-3-319-56478-4) +Proposition 2.1 in [AyJostLeSchwachhoefer:2017](@cite). """ function inner(::ProbabilitySimplex{n,boundary}, p, X, Y) where {n,boundary} d = zero(Base.promote_eltype(p, X, Y)) @@ -368,18 +356,12 @@ end When `vector_at` is `nothing`, return a random (uniform over the Fisher-Rao metric; that is, uniform with respect to the `n`-sphere whose positive orthant is mapped to the simplex). -point `x` on the [`ProbabilitySimplex`](@ref) manifold `M` according to the isometric embedding into -the `n`-sphere by normalizing the vector length of a sample from a multivariate Gaussian. See [^Marsaglia1972]. +point `x` on the [`ProbabilitySimplex`](@ref) manifold `M` according to the isometric embedding into +the `n`-sphere by normalizing the vector length of a sample from a multivariate Gaussian. See [Marsaglia:1972](@cite). When `vector_at` is not `nothing`, return a (Gaussian) random vector from the tangent space -``T_{p}\mathrm{\Delta}^n``by shifting a multivariate Gaussian with standard deviation `σ` +``T_{p}\mathrm{\Delta}^n``by shifting a multivariate Gaussian with standard deviation `σ` to have a zero component sum. - -[^Marsaglia1972]: - > Marsaglia, G.: - > _Choosing a Point from the Surface of a Sphere_. - > Annals of Mathematical Statistics, 43 (2): 645–646, 1972. - > doi: [10.1214/aoms/1177692644](https://doi.org/10.1214/aoms/1177692644) """ rand(::ProbabilitySimplex; σ::Real=1.0) diff --git a/src/manifolds/ProbabilitySimplexEuclideanMetric.jl b/src/manifolds/ProbabilitySimplexEuclideanMetric.jl index 791237a6c6..b44540a3a7 100644 --- a/src/manifolds/ProbabilitySimplexEuclideanMetric.jl +++ b/src/manifolds/ProbabilitySimplexEuclideanMetric.jl @@ -26,17 +26,12 @@ end When `vector_at` is `nothing`, return a random (uniform) point `x` on the [`ProbabilitySimplex`](@ref) with the Euclidean metric -manifold `M` by normalizing independent exponential draws to unit sum, see [^Devroye1986], Theorems 2.1 and 2.2 on p. 207 and 208, respectively. +manifold `M` by normalizing independent exponential draws to unit sum, see [Devroye:1986](@cite), Theorems 2.1 and 2.2 on p. 207 and 208, respectively. When `vector_at` is not `nothing`, return a (Gaussian) random vector from the tangent space -``T_{p}\mathrm{\Delta}^n``by shifting a multivariate Gaussian with standard deviation `σ` +``T_{p}\mathrm{\Delta}^n``by shifting a multivariate Gaussian with standard deviation `σ` to have a zero component sum. -[^Devroye1986]: - > Devroye, L.: - > _Non-Uniform Random Variate Generation_. - > Springer New York, NY, 1986. - > doi: [10.1007/978-1-4613-8643-8](https://doi.org/10.1007/978-1-4613-8643-8) """ rand(::MetricManifold{ℝ,<:ProbabilitySimplex,<:EuclideanMetric}; σ::Real=1.0) diff --git a/src/manifolds/ProductManifold.jl b/src/manifolds/ProductManifold.jl index 197bf9088e..76b6e90c45 100644 --- a/src/manifolds/ProductManifold.jl +++ b/src/manifolds/ProductManifold.jl @@ -1341,6 +1341,31 @@ function representation_size(M::ProductManifold) return (mapreduce(m -> prod(representation_size(m)), +, M.manifolds),) end +@doc raw""" + Y = riemannian_Hessian(M::ProductManifold, p, G, H, X) + riemannian_Hessian!(M::ProductManifold, Y, p, G, H, X) + +Compute the Riemannian Hessian ``\operatorname{Hess} f(p)[X]`` given the +Euclidean gradient ``∇ f(\tilde p)`` in `G` and the Euclidean Hessian ``∇^2 f(\tilde p)[\tilde X]`` in `H`, +where ``\tilde p, \tilde X`` are the representations of ``p,X`` in the embedding,. + +On a product manifold, this decouples and can be computed elementwise. +""" +riemannian_Hessian(M::ProductManifold, p, G, H, X) + +function riemannian_Hessian!(M::ProductManifold, Y, p, G, H, X) + map( + riemannian_Hessian!, + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, G), + submanifold_components(M, H), + submanifold_components(M, X), + ) + return Y +end + @doc raw""" riemann_tensor(M::ProductManifold, p, X, Y, Z) @@ -1663,6 +1688,27 @@ function volume_density(M::ProductManifold, p, X) return prod(dens) end +@doc raw""" + Y = Weingarten(M::ProductManifold, p, X, V) + Weingarten!(M::ProductManifold, Y, p, X, V) + +Since the metric decouples, also the computation of the Weingarten map +``\mathcal W_p`` can be computed elementwise on the single elements of the [`ProductManifold`](@ref) `M`. +""" +Weingarten(::ProductManifold, p, X, V) + +function Weingarten!(M::ProductManifold, Y, p, X, V) + map( + Weingarten!, + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + submanifold_components(M, V), + ) + return Y +end + function zero_vector!(M::ProductManifold, X, p) map( zero_vector!, diff --git a/src/manifolds/ProjectiveSpace.jl b/src/manifolds/ProjectiveSpace.jl index 3a5bc4e6cd..7ff9180fb3 100644 --- a/src/manifolds/ProjectiveSpace.jl +++ b/src/manifolds/ProjectiveSpace.jl @@ -354,12 +354,7 @@ Volume of the ``n``-dimensional [`AbstractProjectiveSpace`](@ref) `M`. The formu ```` where ``Γ`` denotes the [Gamma function](https://en.wikipedia.org/wiki/Gamma_function). -For details see [^BoyaSudarshanTilma2003]. - -[^BoyaSudarshanTilma2003]: - > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports - > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003, - > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1) +For details see [BoyaSudarshanTilma:2003](@cite). """ function manifold_volume(M::AbstractProjectiveSpace{ℝ}) n = manifold_dimension(M) + 1 diff --git a/src/manifolds/QuotientManifold.jl b/src/manifolds/QuotientManifold.jl index 72a833b4d3..3c806f8b09 100644 --- a/src/manifolds/QuotientManifold.jl +++ b/src/manifolds/QuotientManifold.jl @@ -25,7 +25,7 @@ on ``\mathcal N`` we have \mathcal M = \mathcal N / ∼ = \bigl\{ [p] : p ∈ \mathcal N \bigr\}, ``` where ``[p] ≔ \{ q ∈ \mathcal N : q ∼ p\}`` denotes the equivalence class containing ``p``. -For more details see Subsection 3.4.1[^AbsilMahonySepulchre2008]. +For more details see Subsection 3.4.1 [AbsilMahonySepulchre:2008](@cite). This manifold type models an explicit quotient structure. This should be done if either the default implementation of ``\mathcal M`` @@ -42,13 +42,6 @@ it provides a (default) quotient structure that is different from the one introd QuotientManifold(M,N) Create a manifold where `M` is the quotient manifold and `N`is its total space. - -[^AbsilMahonySepulchre2008]: - > Absil, P.-A., Mahony, R. and Sepulchre R., - > _Optimization Algorithms on Matrix Manifolds_ - > Princeton University Press, 2008, - > doi: [10.1515/9781400830244](https://doi.org/10.1515/9781400830244) - > [open access](http://press.princeton.edu/chapters/absil/) """ struct QuotientManifold{𝔽,MT<:AbstractManifold{𝔽},NT<:AbstractManifold} <: AbstractDecoratorManifold{𝔽} @@ -165,7 +158,7 @@ to ``T_q\mathcal N` in place of `Y`. horizontal_lift!(N::AbstractManifold, Y, q, X) """ - horizontal_component(N::AbstractManifold, p, X) + horizontal_component(N::AbstractManifold, p, X) Compute the horizontal component of tangent vector `X` at point `p` in the total space of quotient manifold `N`. @@ -180,7 +173,7 @@ function Base.show(io::IO, M::QuotientManifold) end """ - vertical_component(N::AbstractManifold, p, X) + vertical_component(N::AbstractManifold, p, X) Compute the vertical component of tangent vector `X` at point `p` in the total space of quotient manifold `N`. diff --git a/src/manifolds/Rotations.jl b/src/manifolds/Rotations.jl index ded758a85f..956394822b 100644 --- a/src/manifolds/Rotations.jl +++ b/src/manifolds/Rotations.jl @@ -340,7 +340,7 @@ end parallel_transport_direction(M::Rotations, p, X, d) Compute parallel transport of vector `X` tangent at `p` on the [`Rotations`](@ref) -manifold in the direction `d`. The formula, provided in [^Rentmeesters2011], reads: +manifold in the direction `d`. The formula, provided in [Rentmeesters:2011](@cite), reads: ```math \mathcal P_{q\gets p}X = q^\mathrm{T}p \operatorname{Exp}(d/2) X \operatorname{Exp}(d/2) @@ -348,11 +348,6 @@ manifold in the direction `d`. The formula, provided in [^Rentmeesters2011], rea where ``q=\exp_p d``. The formula simplifies to identity for 2-D rotations. - -[^Rentmeesters2011]: - > Rentmeesters Q., “A gradient method for geodesic data fitting on some symmetric - > Riemannian manifolds,” in 2011 50th IEEE Conference on Decision and Control and - > European Control Conference, Dec. 2011, pp. 7141–7146. doi: [10.1109/CDC.2011.6161280](https://doi.org/10.1109/CDC.2011.6161280). """ parallel_transport_direction(M::Rotations, p, X, d) @@ -390,8 +385,45 @@ function Base.show(io::IO, ::Rotations{n}) where {n} return print(io, "Rotations($(n))") end +@doc raw""" + riemannian_Hessian(M::Rotations, p, G, H, X) + +The Riemannian Hessian can be computed by adopting Eq. (5.6) [Nguyen:2023](@cite), +so very similar to the Stiefel manifold. +The only difference is, that here the tangent vectors are stored +in the Lie algebra, i.e. the update direction is actually ``pX`` instead of just ``X`` (in Stiefel). +and that means the inverse has to be appliead to the (Euclidean) Hessian +to map it into the Lie algebra. +""" +riemannian_Hessian(M::Rotations, p, G, H, X) +function riemannian_Hessian!(::Rotations{N}, Y, p, G, H, X) where {N} + symmetrize!(Y, G' * p) + project!(SkewSymmetricMatrices(N), Y, p' * H - X * Y) + return Y +end + Distributions.support(d::NormalRotationDistribution) = MPointSupport(d.manifold) +@doc raw""" + Weingarten(M::Rotations, p, X, V) + +Compute the Weingarten map ``\mathcal W_p`` at `p` on the [`Stiefel`](@ref) `M` with respect to the +tangent vector ``X \in T_p\mathcal M`` and the normal vector ``V \in N_p\mathcal M``. + +The formula is due to [AbsilMahonyTrumpf:2013](@cite) given by + +```math +\mathcal W_p(X,V) = -\frac{1}{2}p\bigl(V^{\mathrm{T}}X - X^\mathrm{T}V\bigr) +``` +""" +Weingarten(::Rotations, p, X, V) + +function Weingarten!(::Rotations, Y, p, X, V) + Y .= V' * X + Y .= -p * 1 / 2 * (Y - Y') + return Y +end + @doc raw""" zero_vector(M::Rotations, p) diff --git a/src/manifolds/SkewHermitian.jl b/src/manifolds/SkewHermitian.jl index 56c6485d7d..234e3f295d 100644 --- a/src/manifolds/SkewHermitian.jl +++ b/src/manifolds/SkewHermitian.jl @@ -252,3 +252,16 @@ end function Base.show(io::IO, ::SkewSymmetricMatrices{n}) where {n} return print(io, "SkewSymmetricMatrices($(n))") end + +@doc raw""" + Y = Weingarten(M::SkewSymmetricMatrices, p, X, V) + Weingarten!(M::SkewSymmetricMatrices, Y, p, X, V) + +Compute the Weingarten map ``\mathcal W_p`` at `p` on the [`SkewSymmetricMatrices`](@ref) `M` with respect to the +tangent vector ``X \in T_p\mathcal M`` and the normal vector ``V \in N_p\mathcal M``. + +Since this a flat space by itself, the result is always the zero tangent vector. +""" +Weingarten(::SkewSymmetricMatrices, p, X, V) + +Weingarten!(::SkewSymmetricMatrices, Y, p, X, V) = fill!(Y, 0) diff --git a/src/manifolds/Spectrahedron.jl b/src/manifolds/Spectrahedron.jl index fb3a0e043a..28b883538f 100644 --- a/src/manifolds/Spectrahedron.jl +++ b/src/manifolds/Spectrahedron.jl @@ -34,20 +34,13 @@ endowed with the [`Euclidean`](@ref) metric from the embedding, i.e. from the $ This manifold was for example -investigated in[^JourneeBachAbsilSepulchre2010]. +investigated in [JourneeBachAbsilSepulchre:2010](@cite). # Constructor Spectrahedron(n,k) generates the manifold $\mathcal S(n,k) \subset ℝ^{n × n}$. - -[^JourneeBachAbsilSepulchre2010]: - > Journée, M., Bach, F., Absil, P.-A., and Sepulchre, R.: - > “Low-Rank Optimization on the Cone of Positive Semidefinite Matrices”, - > SIAM Journal on Optimization (20)5, pp. 2327–2351, 2010. - > doi: [10.1137/080731359](https://doi.org/10.1137/080731359), - > arXiv: [0807.4423](http://arxiv.org/abs/0807.4423). """ struct Spectrahedron{N,K} <: AbstractDecoratorManifold{ℝ} end diff --git a/src/manifolds/Sphere.jl b/src/manifolds/Sphere.jl index 7d91c0432b..54e0fd38cb 100644 --- a/src/manifolds/Sphere.jl +++ b/src/manifolds/Sphere.jl @@ -513,17 +513,12 @@ end riemann_tensor(M::AbstractSphere{ℝ}, p, X, Y, Z) Compute the Riemann tensor ``R(X,Y)Z`` at point `p` on [`AbstractSphere`](@ref) `M`. -The formula reads [^MuralidharanFletcher2021] (though note that a different convention is +The formula reads [MuralidharanFlecther:2012](@cite) (though note that a different convention is used in that paper than in Manifolds.jl): ````math R(X,Y)Z = \langle Z, Y \rangle X - \langle Z, X \rangle Y ```` - -[^MuralidharanFletcher2021]: - > P. Muralidharan and P. T. Fletcher, “Sasaki Metrics for Analysis of Longitudinal Data - > on Manifolds,” Proc IEEE Comput Soc Conf Comput Vis Pattern Recognit, vol. 2012, - > pp. 1027–1034, Jun. 2012, doi: [10.1109/CVPR.2012.6247780](https://doi.org/10.1109/CVPR.2012.6247780). """ riemann_tensor(M::AbstractSphere{ℝ}, p, X, Y, Z) @@ -539,12 +534,7 @@ end Compute volume density function of a sphere, i.e. determinant of the differential of exponential map `exp(M, p, X)`. The formula reads ``(\sin(\lVert X\rVert)/\lVert X\rVert)^(n-1)`` -where `n` is the dimension of `M`. It is derived from Eq. (4.1) in [^ChevallierLiLuDunson2022]. - -[^ChevallierLiLuDunson2022]: - > E. Chevallier, D. Li, Y. Lu, and D. B. Dunson, “Exponential-wrapped distributions on - > symmetric spaces.” arXiv, Oct. 09, 2022. - > doi: [10.48550/arXiv.2009.01983](https://doi.org/10.48550/arXiv.2009.01983). +where `n` is the dimension of `M`. It is derived from Eq. (4.1) in [ChevallierLiLuDunson:2022](@cite). """ function volume_density(M::AbstractSphere{ℝ}, p, X) Xnorm = norm(X) @@ -552,6 +542,26 @@ function volume_density(M::AbstractSphere{ℝ}, p, X) return usinc(Xnorm)^n end +@doc raw""" + Y = Weingarten(M::Sphere, p, X, V) + Weingarten!(M::Sphere, Y, p, X, V) + +Compute the Weingarten map ``\mathcal W_p`` at `p` on the [`Sphere`](@ref) `M` with respect to the +tangent vector ``X \in T_p\mathcal M`` and the normal vector ``V \in N_p\mathcal M``. + +The formula is due to [AbsilMahonyTrumpf:2013](@cite) given by + +```math +\mathcal W_p(X,V) = -Xp^{\mathrm{T}}V +``` +""" +Weingarten(::Sphere, p, X, V) + +function Weingarten!(::Sphere, Y, p, X, V) + Y .= -dot(p, V) .* X + return Y +end + """ StereographicAtlas() diff --git a/src/manifolds/Stiefel.jl b/src/manifolds/Stiefel.jl index 4f753edebe..2cbee11ac0 100644 --- a/src/manifolds/Stiefel.jl +++ b/src/manifolds/Stiefel.jl @@ -172,12 +172,7 @@ inverse_retract(::Stiefel, ::Any, ::Any, ::PolarInverseRetraction) Compute the inverse retraction based on a qr decomposition for two points `p`, `q` on the [`Stiefel`](@ref) manifold `M` and return the resulting tangent vector in `X`. The computation follows Algorithm 1 -in [^KanekoFioriTanaka2013]. - -[^KanekoFioriTanaka2013]: - > T. Kaneko, S. Fiori, T. Tanaka: "Empirical Arithmetic Averaging over the - > Compact Stiefel AbstractManifold", IEEE Transactions on Signal Processing, 2013, - > doi: [10.1109/TSP.2012.2226167](https://doi.org/10.1109/TSP.2012.2226167). +in [KanekoFioriTanaka:2013](@cite). """ inverse_retract(::Stiefel, ::Any, ::Any, ::QRInverseRetraction) @@ -342,11 +337,11 @@ end @doc raw""" retract(::Stiefel, p, X, ::CayleyRetraction) -Compute the retraction on the [`Stiefel`](@ref) that is based on the Cayley transform[^Zhu2017]. +Compute the retraction on the [`Stiefel`](@ref) that is based on the Cayley transform[Zhu:2016](@cite). Using ````math W_{p,X} = \operatorname{P}_pXp^{\mathrm{H}} - pX^{\mathrm{H}}\operatorname{P_p} - \quad\text{where}  + \quad\text{where} \operatorname{P}_p = I - \frac{1}{2}pp^{\mathrm{H}} ```` the formula reads @@ -355,46 +350,44 @@ the formula reads ```` It is implemented as the case $m=1$ of the `PadeRetraction`. - -[^Zhu2017]: - > X. Zhu: - > A Riemannian conjugate gradient method for optimizazion on the Stiefel manifold, - > Computational Optimization and Applications 67(1), pp. 73–110, 2017. - > doi [10.1007/s10589-016-9883-4](https://doi.org/10.1007/s10589-016-9883-4). """ retract(::Stiefel, ::Any, ::Any, ::CayleyRetraction) @doc raw""" retract(M::Stiefel, p, X, ::PadeRetraction{m}) -Compute the retraction on the [`Stiefel`](@ref) manifold `M` based on the Padé approximation of order $m$[^ZhuDuan2018]. +Compute the retraction on the [`Stiefel`](@ref) manifold `M` based on the Padé approximation of order $m$ [ZhuDuan:2018](@cite). Let $p_m$ and $q_m$ be defined for any matrix $A ∈ ℝ^{n×x}$ as + ````math p_m(A) = \sum_{k=0}^m \frac{(2m-k)!m!}{(2m)!(m-k)!}\frac{A^k}{k!} ```` + and + ````math q_m(A) = \sum_{k=0}^m \frac{(2m-k)!m!}{(2m)!(m-k)!}\frac{(-A)^k}{k!} ```` + respectively. Then the Padé approximation (of the matrix exponential $\exp(A)$) reads + ````math r_m(A) = q_m(A)^{-1}p_m(A) ```` + Defining further + ````math W_{p,X} = \operatorname{P}_pXp^{\mathrm{H}} - pX^{\mathrm{H}}\operatorname{P_p} \quad\text{where } \operatorname{P}_p = I - \frac{1}{2}pp^{\mathrm{H}} ```` + the retraction reads + ````math \operatorname{retr}_pX = r_m(W_{p,X})p ```` -[^ZhuDuan2018]: - > X. Zhu, C. Duan: - > On matrix exponentials and their approximations related to optimization on the Stiefel manifold, - > Optimizazion Letters 13(5), pp. 1069–1083, 2018. - > doi [10.1007/s11590-018-1341-z](https://doi.org/10.1007/s11590-018-1341-z). """ retract(::Stiefel, ::Any, ::Any, ::PadeRetraction) @@ -489,12 +482,8 @@ Uniform distribution on given (real-valued) [`Stiefel`](@ref) `M`. Specifically, this is the normalized Haar and Hausdorff measure on `M`. Generated points will be of similar type as `p`. -The implementation is based on Section 2.5.1 in [^Chikuse2003]; -see also Theorem 2.2.1(iii) in [^Chikuse2003]. - -[^Chikuse2003]: - > Y. Chikuse: "Statistics on Special Manifolds", Springer New York, 2003, - > doi: [10.1007/978-0-387-21540-2](https://doi.org/10.1007/978-0-387-21540-2). +The implementation is based on Section 2.5.1 in [Chikuse:2003](@cite); +see also Theorem 2.2.1(iii) in [Chikuse:2003](@cite). """ function uniform_distribution(M::Stiefel{n,k,ℝ}, p) where {n,k} μ = Distributions.Zeros(n, k) @@ -509,7 +498,7 @@ end @doc raw""" vector_transport_direction(::Stiefel, p, X, d, ::DifferentiatedRetractionVectorTransport{CayleyRetraction}) -Compute the vector transport given by the differentiated retraction of the [`CayleyRetraction`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.CayleyRetraction), cf. [^Zhu2017] Equation (17). +Compute the vector transport given by the differentiated retraction of the [`CayleyRetraction`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.CayleyRetraction), cf. [Zhu:2016](@cite) Equation (17). The formula reads ````math @@ -538,7 +527,7 @@ vector_transport_direction( vector_transport_direction(M::Stiefel, p, X, d, DifferentiatedRetractionVectorTransport{PolarRetraction}) Compute the vector transport by computing the push forward of -[`retract(::Stiefel, ::Any, ::Any, ::PolarRetraction)`](@ref) Section 3.5 of [^Zhu2017]: +[`retract(::Stiefel, ::Any, ::Any, ::PolarRetraction)`](@ref) Section 3.5 of [Zhu:2016](@cite): ```math T_{p,d}^{\text{Pol}}(X) = q*Λ + (I-qq^{\mathrm{T}})X(1+d^\mathrm{T}d)^{-\frac{1}{2}}, @@ -563,7 +552,7 @@ vector_transport_direction( Compute the vector transport by computing the push forward of the [`retract(::Stiefel, ::Any, ::Any, ::QRRetraction)`](@ref), -See [^AbsilMahonySepulchre2008], p. 173, or Section 3.5 of [^Zhu2017]. +See [AbsilMahonySepulchre:2008](@cite), p. 173, or Section 3.5 of [Zhu:2016](@cite). ```math T_{p,d}^{\text{QR}}(X) = q*\rho_{\mathrm{s}}(q^\mathrm{T}XR^{-1}) + (I-qq^{\mathrm{T}})XR^{-1}, ``` @@ -577,12 +566,6 @@ A_{ij}&\text{ if } i > j\\ -A_{ji} \text{ if } i < j.\\ \end{cases} ``` -[^AbsilMahonySepulchre2008]: - >Absil, P.-A., Mahony, R. and Sepulchre R., - > _Optimization Algorithms on Matrix Manifolds_ - > Princeton University Press, 2008, - > doi: [10.1515/9781400830244](https://doi.org/10.1515/9781400830244) - > [open access](http://press.princeton.edu/chapters/absil/) """ vector_transport_direction( ::Stiefel, @@ -628,7 +611,7 @@ end Compute the vector transport by computing the push forward of the [`retract(M::Stiefel, ::Any, ::Any, ::PolarRetraction)`](@ref), see -Section 4 of [^HuangGallivanAbsil2015] or Section 3.5 of [^Zhu2017]: +Section 4 of [HuangGallivanAbsil:2015](@cite) or Section 3.5 of [Zhu:2016](@cite): ```math T_{q\gets p}^{\text{Pol}}(X) = q*Λ + (I-qq^{\mathrm{T}})X(1+d^\mathrm{T}d)^{-\frac{1}{2}}, @@ -640,12 +623,6 @@ and $Λ$ is the unique solution of the Sylvester equation ```math Λ(I+d^\mathrm{T}d)^{\frac{1}{2}} + (I + d^\mathrm{T}d)^{\frac{1}{2}} = q^\mathrm{T}X - X^\mathrm{T}q ``` -[^HuangGallivanAbsil2015]: - > Huang, W., Gallivan, K. A., and Absil, P.-A.: - > _A Broyden class of quasi-Newton methods for Riemannian optimization_ - > SIAM Journal of Optimization, 2015, Vol. 25, No. 3, pp. 1660–1685 - > doi: [10.1137/140955483](https://doi.org/10.1137/140955483) - > pdf: [tech. report](https://www.math.fsu.edu/~whuang2/pdf/RBroydenBasic_techrep.pdf) """ vector_transport_to( ::Stiefel, @@ -660,7 +637,7 @@ vector_transport_to( Compute the vector transport by computing the push forward of the [`retract(M::Stiefel, ::Any, ::Any, ::QRRetraction)`](@ref), -see [^AbsilMahonySepulchre2008], p. 173, or Section 3.5 of [^Zhu2017]. +see [AbsilMahonySepulchre:2008](@cite), p. 173, or Section 3.5 of [Zhu:2016](@cite). ```math T_{q \gets p}^{\text{QR}}(X) = q*\rho_{\mathrm{s}}(q^\mathrm{T}XR^{-1}) + (I-qq^{\mathrm{T}})XR^{-1}, diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl index d27b0bec51..68d57a6f28 100644 --- a/src/manifolds/StiefelCanonicalMetric.jl +++ b/src/manifolds/StiefelCanonicalMetric.jl @@ -2,14 +2,8 @@ CanonicalMetric <: AbstractMetric The Canonical Metric refers to a metric for the [`Stiefel`](@ref) -manifold, see[^EdelmanAriasSmith1998]. - -[^EdelmanAriasSmith1998]: - > Edelman, A., Ariar, T. A., Smith, S. T.: - > _The Geometry of Algorihthms with Orthogonality Constraints_, - > SIAM Journal on Matrix Analysis and Applications (20(2), pp. 303–353, 1998. - > doi: [10.1137/S0895479895290954](https://doi.org/10.1137/S0895479895290954) - > arxiv: [9806030](https://arxiv.org/abs/physics/9806030) +manifold, see[EdelmanAriasSmith:1998](@cite). + """ struct CanonicalMetric <: RiemannianMetric end @@ -68,7 +62,7 @@ the exponential map reads ```math q = \exp_p X = pC + QB. ``` -For more details, see [^EdelmanAriasSmith1998][^Zimmermann2017]. +For more details, see [EdelmanAriasSmith:1998](@cite)[Zimmermann:2017](@cite). """ exp(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, ::Any...) where {n,k} @@ -111,14 +105,8 @@ Compute an approximation to the logarithmic map on the [`Stiefel`](@ref)`(n,k)` using a matrix-algebraic based approach to an iterative inversion of the formula of the [`exp`](@ref exp(::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, ::Any...) where {n,k}). -The algorithm is derived in[^Zimmermann2017] and it uses the `max_iterations` and the `tolerance` field +The algorithm is derived in [Zimmermann:2017](@cite) and it uses the `max_iterations` and the `tolerance` field from the [`ApproximateLogarithmicMap`](@ref). - -[^Zimmermann2017]: - > Zimmermann, R.: _A matrix-algebraic algorithm for the Riemannian logarithm on the Stiefel manifold under the canoncial metric. - > SIAM Journal on Matrix Analysis and Applications 28(2), pp. 322-342, 2017. - > doi: [10.1137/16M1074485](https://doi.org/10.1137/16M1074485), - > arXiv: [1604.05054](https://arxiv.org/abs/1604.05054). """ inverse_retract( ::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, @@ -186,3 +174,46 @@ function inverse_retract!( @views mul!(X, qfact.U, LV[1:(2k), 1:k]) return X end + +@doc raw""" + Y = riemannian_Hessian(M::MetricManifold{ℝ, Stiefel{n,k}, CanonicalMetric}, p, G, H, X) + riemannian_Hessian!(M::MetricManifold{ℝ, Stiefel{n,k}, CanonicalMetric}, Y, p, G, H, X) + +Compute the Riemannian Hessian ``\operatorname{Hess} f(p)[X]`` given the +Euclidean gradient ``∇ f(\tilde p)`` in `G` and the Euclidean Hessian ``∇^2 f(\tilde p)[\tilde X]`` in `H`, +where ``\tilde p, \tilde X`` are the representations of ``p,X`` in the embedding,. + +Here, we adopt Eq. (5.6) [Nguyen:2023](@cite), for the [`CanonicalMetric`](@ref) +``α_0=1, α_1=\frac{1}{2}`` in their formula. The formula reads + +```math + \operatorname{Hess}f(p)[X] + = + \operatorname{proj}_{T_p\mathcal M}\Bigl( + ∇^2f(p)[X] - \frac{1}{2} X \bigl( (∇f(p))^{\mathrm{H}}p + p^{\mathrm{H}}∇f(p)\bigr) + - \frac{1}{2} \bigl( P ∇f(p) p^{\mathrm{H}} + p ∇f(p))^{\mathrm{H}} P)X + \Bigr), +``` +where ``P = I-pp^{\mathrm{H}}``. +""" +riemannian_Hessian( + M::MetricManifold{𝔽,Stiefel{n,k,𝔽},CanonicalMetric}, + p, + G, + H, + X, +) where {n,k,𝔽} + +function riemannian_Hessian!( + M::MetricManifold{𝔽,Stiefel{n,k,𝔽},CanonicalMetric}, + Y, + p, + G, + H, + X, +) where {n,k,𝔽} + Gp = symmetrize(G' * p) + Z = symmetrize((I - p * p') * G * p') + project!(M, Y, p, H - X * Gp - Z * X) + return Y +end diff --git a/src/manifolds/StiefelEuclideanMetric.jl b/src/manifolds/StiefelEuclideanMetric.jl index eba7b5d2ce..f897f4053b 100644 --- a/src/manifolds/StiefelEuclideanMetric.jl +++ b/src/manifolds/StiefelEuclideanMetric.jl @@ -136,8 +136,8 @@ end project(M::Stiefel,p) Projects `p` from the embedding onto the [`Stiefel`](@ref) `M`, i.e. compute `q` -as the polar decomposition of $p$ such that $q^{\mathrm{H}q$ is the identity, -where $\cdot^{\mathrm{H}}$ denotes the hermitian, i.e. complex conjugate transposed. +as the polar decomposition of $p$ such that ``q^{\mathrm{H}}q`` is the identity, +where ``\cdot^{\mathrm{H}}`` denotes the hermitian, i.e. complex conjugate transposed. """ project(::Stiefel, ::Any, ::Any) @@ -154,7 +154,7 @@ Project `X` onto the tangent space of `p` to the [`Stiefel`](@ref) manifold `M`. The formula reads ````math -\operatorname{proj}_{\mathcal M}(p, X) = X - p \operatorname{Sym}(p^{\mathrm{H}}X), +\operatorname{proj}_{T_p\mathcal M}(X) = X - p \operatorname{Sym}(p^{\mathrm{H}}X), ```` where $\operatorname{Sym}(q)$ is the symmetrization of $q$, e.g. by @@ -185,6 +185,55 @@ function retract_project!(M::Stiefel, q, p, X, t::Number) return q end +@doc raw""" + Y = riemannian_Hessian(M::Stiefel, p, G, H, X) + riemannian_Hessian!(M::Stiefel, Y, p, G, H, X) + +Compute the Riemannian Hessian ``\operatorname{Hess} f(p)[X]`` given the +Euclidean gradient ``∇ f(\tilde p)`` in `G` and the Euclidean Hessian ``∇^2 f(\tilde p)[\tilde X]`` in `H`, +where ``\tilde p, \tilde X`` are the representations of ``p,X`` in the embedding,. + +Here, we adopt Eq. (5.6) [Nguyen:2023](@cite), where we use for the [`EuclideanMetric`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.EuclideanMetric) +``α_0=α_1=1`` in their formula. Then the formula reads + +```math + \operatorname{Hess}f(p)[X] + = + \operatorname{proj}_{T_p\mathcal M}\Bigl( + ∇^2f(p)[X] - \frac{1}{2} X \bigl((∇f(p))^{\mathrm{H}}p + p^{\mathrm{H}}∇f(p)\bigr) + \Bigr). +``` + +Compared to Eq. (5.6) also the metric conversion simplifies to the identity. +""" +riemannian_Hessian(M::Stiefel, p, G, H, X) + +function riemannian_Hessian!(M::Stiefel, Y, p, G, H, X) + Z = symmetrize(G' * p) + project!(M, Y, p, H - X * Z) + return Y +end + function vector_transport_to!(M::Stiefel, Y, ::Any, X, q, ::ProjectionTransport) return project!(M, Y, q, X) end + +@doc raw""" + Weingarten(M::Stiefel, p, X, V) + +Compute the Weingarten map ``\mathcal W_p`` at `p` on the [`Stiefel`](@ref) `M` with respect to the +tangent vector ``X \in T_p\mathcal M`` and the normal vector ``V \in N_p\mathcal M``. + +The formula is due to [AbsilMahonyTrumpf:2013](@cite) given by + +```math +\mathcal W_p(X,V) = -Xp^{\mathrm{H}}V - \frac{1}{2}p\bigl(X^\mathrm{H}V + V^{\mathrm{H}}X\bigr) +``` +""" +Weingarten(::Stiefel, p, X, V) + +function Weingarten!(::Stiefel, Y, p, X, V) + Z = symmetrize(X' * V) + Y .= -X * p' * V .- p * Z + return Y +end diff --git a/src/manifolds/StiefelSubmersionMetric.jl b/src/manifolds/StiefelSubmersionMetric.jl index 1b67c6cf6a..4c0fbbe00e 100644 --- a/src/manifolds/StiefelSubmersionMetric.jl +++ b/src/manifolds/StiefelSubmersionMetric.jl @@ -7,18 +7,8 @@ The family, with a single real parameter ``α>-1``, has two special cases: - ``α = -\frac{1}{2}``: [`EuclideanMetric`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.EuclideanMetric) - ``α = 0``: [`CanonicalMetric`](@ref) -The family was described in [^HüperMarkinaLeite2021]. This implementation follows the -description in [^ZimmermanHüper2022]. - -[^HüperMarkinaLeite2021]: - > Hüper, M., Markina, A., Leite, R. T. (2021) - > "A Lagrangian approach to extremal curves on Stiefel manifolds" - > Journal of Geometric Mechanics, 13(1): 55-72. - > doi: [10.3934/jgm.2020031](http://dx.doi.org/10.3934/jgm.2020031) -[^ZimmermanHüper2022]: - > Ralf Zimmerman and Knut Hüper. (2022). - > "Computing the Riemannian logarithm on the Stiefel manifold: metrics, methods and performance." - > arXiv: [2103.12046](https://arxiv.org/abs/2103.12046) +The family was described in [HueperMarkinaSilvaLeite:2021](@cite). This implementation follows the +description in [ZimmermannHueper:2022](@cite). # Constructor @@ -45,7 +35,7 @@ The exponential map is given by X p^\mathrm{T} - p X^\mathrm{T} \bigr) p \operatorname{Exp}\bigl(\frac{\alpha}{\alpha+1} p^\mathrm{T} X\bigr) ```` -This implementation is based on [^ZimmermanHüper2022]. +This implementation is based on [ZimmermannHueper:2022](@cite). For ``k < \frac{n}{2}`` the exponential is computed more efficiently using [`StiefelFactorization`](@ref). @@ -207,8 +197,8 @@ end Compute the logarithmic map on the [`Stiefel(n,k)`](@ref) manifold with respect to the [`StiefelSubmersionMetric`](@ref). The logarithmic map is computed using [`ShootingInverseRetraction`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.ShootingInverseRetraction). For -``k ≤ \lfloor\frac{n}{2}\rfloor``, this is sped up using the ``k``-shooting method of -[^ZimmermanHüper2022]. Keyword arguments are forwarded to `ShootingInverseRetraction`; see +``k ≤ \lfloor\frac{n}{2}\rfloor``, this is sped up using the ``k``-shooting method of [ZimmermannHueper:2022](@cite). +Keyword arguments are forwarded to `ShootingInverseRetraction`; see that documentation for details. Their defaults are: - `num_transport_points=4` - `tolerance=sqrt(eps())` @@ -257,13 +247,59 @@ function log!( return inverse_retract!(M, X, p, q, inverse_retraction) end +@doc raw""" + Y = riemannian_Hessian(M::MetricManifold{ℝ,Stiefel{n,k,ℝ}, StiefelSubmersionMetric},, p, G, H, X) + riemannian_Hessian!(MetricManifold{ℝ,Stiefel{n,k,ℝ}, StiefelSubmersionMetric},, Y, p, G, H, X) + +Compute the Riemannian Hessian ``\operatorname{Hess} f(p)[X]`` given the +Euclidean gradient ``∇ f(\tilde p)`` in `G` and the Euclidean Hessian ``∇^2 f(\tilde p)[\tilde X]`` in `H`, +where ``\tilde p, \tilde X`` are the representations of ``p,X`` in the embedding,. + +Here, we adopt Eq. (5.6) [Nguyen:2023](@cite), for the [`CanonicalMetric`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.EuclideanMetric) +``α_0=1, α_1=\frac{1}{2}`` in their formula. The formula reads + +```math + \operatorname{Hess}f(p)[X] + = + \operatorname{proj}_{T_p\mathcal M}\Bigl( + ∇^2f(p)[X] - \frac{1}{2} X \bigl( (∇f(p))^{\mathrm{H}}p + p^{\mathrm{H}}∇f(p)\bigr) + - \frac{2α+1}{2(α+1)} \bigl( P ∇f(p) p^{\mathrm{H}} + p ∇f(p))^{\mathrm{H}} P)X + \Bigr), +``` +where ``P = I-pp^{\mathrm{H}}``. + +Compared to Eq. (5.6) we have that their ``α_0 = 1``and ``\alpha_1 = \frac{2α+1}{2(α+1)} + 1``. +""" +riemannian_Hessian( + M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, + p, + G, + H, + X, +) where {n,k} + +function riemannian_Hessian!( + M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, + Y, + p, + G, + H, + X, +) where {n,k} + α = metric(M).α + Gp = symmetrize(G' * p) + Z = symmetrize((I - p * p') * G * p') + project!(M, Y, p, H - X * Gp - (2 * α + 1) / (α + 1) * Z * X) + return Y +end + # StiefelFactorization code # Note: intended only for internal use @doc raw""" StiefelFactorization{UT,XT} <: AbstractManifoldPoint -Represent points (and vectors) on `Stiefel(n, k)` with ``2k × k`` factors.[^ZimmermanHüper2022] +Represent points (and vectors) on `Stiefel(n, k)` with ``2k × k`` factors [ZimmermannHueper:2022](@cite). Given a point ``p ∈ \mathrm{St}(n, k)`` and another matrix ``B ∈ ℝ^{n × k}`` for ``k ≤ \lfloor\frac{n}{2}\rfloor`` the factorization is diff --git a/src/manifolds/Symmetric.jl b/src/manifolds/Symmetric.jl index c64d651c93..6bad4bcb9c 100644 --- a/src/manifolds/Symmetric.jl +++ b/src/manifolds/Symmetric.jl @@ -233,3 +233,16 @@ project!(M::SymmetricMatrices, Y, p, X) = (Y .= (X .+ transpose(X)) ./ 2) function Base.show(io::IO, ::SymmetricMatrices{n,F}) where {n,F} return print(io, "SymmetricMatrices($(n), $(F))") end + +@doc raw""" + Y = Weingarten(M::SymmetricMatrices, p, X, V) + Weingarten!(M::SymmetricMatrices, Y, p, X, V) + +Compute the Weingarten map ``\mathcal W_p`` at `p` on the [`SymmetricMatrices`](@ref) `M` with respect to the +tangent vector ``X \in T_p\mathcal M`` and the normal vector ``V \in N_p\mathcal M``. + +Since this a flat space by itself, the result is always the zero tangent vector. +""" +Weingarten(::SymmetricMatrices, p, X, V) + +Weingarten!(::SymmetricMatrices, Y, p, X, V) = fill!(Y, 0) diff --git a/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl b/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl index c6987940bb..dde2356aef 100644 --- a/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl +++ b/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl @@ -419,19 +419,39 @@ function parallel_transport_to!(M::SymmetricPositiveDefinite{N}, Y, p, X, q) whe return copyto!(Y, vtp) end +@doc raw""" + riemannian_Hessian(M::SymmetricPositiveDefinite, p, G, H, X) + +The Riemannian Hessian can be computed as stated in Eq. (7.3) [Nguyen:2023](@cite). +Let ``\nabla f(p)`` denote the Euclidean gradient `G`, +``\nabla^2 f(p)[X]`` the Euclidean Hessian `H`, and +``\operatorname{sym}(X) = \frac{1}{2}\bigl(X^{\mathrm{T}}+X\bigr)`` +the symmetrization operator. Then the formula reads + +```math + \operatorname{Hess}f(p)[X] + = + p\operatorname{sym}(∇^2 f(p)[X])p + + \operatorname{sym}\bigl( X\operatorname{sym}\bigl(∇ f(p)\bigr)p) +``` +""" +riemannian_Hessian(M::SymmetricPositiveDefinite, p, G, H, X) + +function riemannian_Hessian!(::SymmetricPositiveDefinite, Y, p, G, H, X) + symmetrize!(Y, G) + symmetrize!(Y, X * Y * p) + Z = symmetrize(H) + Y .+= p * Z * p + return Y +end + @doc raw""" riemann_tensor(::SymmetricPositiveDefinite, p, X, Y, Z) Compute the value of Riemann tensor on the [`SymmetricPositiveDefinite`](@ref) manifold. -The formula reads[^Rentmeesters2011] ``R(X,Y)Z=p^{1/2}R(X_I, Y_I)Z_Ip^{1/2}``, where +The formula reads [Rentmeesters:2011](@cite) ``R(X,Y)Z=p^{1/2}R(X_I, Y_I)Z_Ip^{1/2}``, where ``R_I(X_I, Y_I)Z_I=\frac{1}{4}[Z_I, [X_I, Y_I]]``, ``X_I=p^{-1/2}Xp^{-1/2}``, ``Y_I=p^{-1/2}Yp^{-1/2}`` and ``Z_I=p^{-1/2}Zp^{-1/2}``. - -[^Rentmeesters2011]: - > Q. Rentmeesters, “A gradient method for geodesic data fitting on some symmetric - > Riemannian manifolds,” in 2011 50th IEEE Conference on Decision and Control and - > European Control Conference, Dec. 2011, pp. 7141–7146. - > doi: [10.1109/CDC.2011.6161280](https://doi.org/10.1109/CDC.2011.6161280). """ riemann_tensor(::SymmetricPositiveDefinite, p, X, Y, Z) @@ -450,14 +470,8 @@ end volume_density(::SymmetricPositiveDefinite{n}, p, X) where {n} Compute the volume density of the [`SymmetricPositiveDefinite`](@ref) manifold at `p` -in direction `X`. See [^ChevallierKalungaAngulo2017], Section 6.2 for details. +in direction `X`. See [ChevallierKalungaAngulo:2017](@cite), Section 6.2 for details. Note that metric in Manifolds.jl has a different scaling factor than the reference. - -[^ChevallierKalungaAngulo2017]: - > E. Chevallier, E. Kalunga, and J. Angulo, “Kernel Density Estimation on Spaces of - > Gaussian Distributions and Symmetric Positive Definite Matrices,” SIAM J. Imaging Sci., - > vol. 10, no. 1, pp. 191–215, Jan. 2017, - > doi: [10.1137/15M1053566](https://doi.org/10.1137/15M1053566). """ function volume_density(::SymmetricPositiveDefinite{n}, p, X) where {n} eig = eigvals(X) diff --git a/src/manifolds/SymmetricPositiveDefiniteBuresWasserstein.jl b/src/manifolds/SymmetricPositiveDefiniteBuresWasserstein.jl index b41c3531ec..95e06819bb 100644 --- a/src/manifolds/SymmetricPositiveDefiniteBuresWasserstein.jl +++ b/src/manifolds/SymmetricPositiveDefiniteBuresWasserstein.jl @@ -1,13 +1,7 @@ @doc raw""" BurresWassertseinMetric <: AbstractMetric -The Bures Wasserstein metric for symmetric positive definite matrices[^MalagoMontruccioPistone2018]. - -[^MalagoMontruccioPistone2018]: - > Malagò, L., Montrucchio, L., Pistone, G.: - > _Wasserstein Riemannian geometry of Gaussian densities_. - > Information Geometry, 1, pp. 137–179, 2018. - > doi: [10.1007/s41884-018-0014-4](https://doi.org/10.1007/s41884-018-0014-4) +The Bures Wasserstein metric for symmetric positive definite matrices [MalagoMontruccioPistone:2018](@cite) """ struct BuresWassersteinMetric <: RiemannianMetric end diff --git a/src/manifolds/SymmetricPositiveDefiniteGeneralizedBuresWasserstein.jl b/src/manifolds/SymmetricPositiveDefiniteGeneralizedBuresWasserstein.jl index 46fd24e117..b710406fa3 100644 --- a/src/manifolds/SymmetricPositiveDefiniteGeneralizedBuresWasserstein.jl +++ b/src/manifolds/SymmetricPositiveDefiniteGeneralizedBuresWasserstein.jl @@ -1,15 +1,10 @@ @doc raw""" GeneralizedBurresWassertseinMetric{T<:AbstractMatrix} <: AbstractMetric -The generalized Bures Wasserstein metric for symmetric positive definite matrices, see[^HanMishraJawanpuriaGao2021]. +The generalized Bures Wasserstein metric for symmetric positive definite matrices, see [HanMishraJawanpuriaGao:2021](@cite). This metric internally stores the symmetric positive definite matrix ``M`` to generalise the metric, where the name also follows the mentioned preprint. - -[^HanMishraJawanpuriaGao2021]: - > Han, A., Mishra, B., Jawanpuria, P., Gao, J.: - > _Generalized Bures-Wasserstein geometry for positive definite matrices_. - > arXiv: [2110.10464](https://arxiv.org/abs/2110.10464). """ struct GeneralizedBuresWassersteinMetric{T<:AbstractMatrix} <: RiemannianMetric M::T diff --git a/src/manifolds/SymmetricPositiveDefiniteLogCholesky.jl b/src/manifolds/SymmetricPositiveDefiniteLogCholesky.jl index 0ed0d57d6d..703dc0b241 100644 --- a/src/manifolds/SymmetricPositiveDefiniteLogCholesky.jl +++ b/src/manifolds/SymmetricPositiveDefiniteLogCholesky.jl @@ -2,11 +2,7 @@ LogCholeskyMetric <: RiemannianMetric The Log-Cholesky metric imposes a metric based on the Cholesky decomposition as -introduced by [^Lin2019]. - -[^Lin2019]: - > Lin, Zenhua: "Riemannian Geometry of Symmetric Positive Definite Matrices via - > Cholesky Decomposition", arXiv: [1908.09326](https://arxiv.org/abs/1908.09326). +introduced by [Lin:2019](@cite). """ struct LogCholeskyMetric <: RiemannianMetric end diff --git a/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl b/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl index dd3ab4d89c..c1232f4d5b 100644 --- a/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl +++ b/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl @@ -31,7 +31,7 @@ X ∈ 𝔽^{n × n}\,|\,X = qY^{\mathrm{H}} + Yq^{\mathrm{H}} ```` Note that the metric used yields a non-complete manifold. -The metric was used in[^JourneeBachAbsilSepulchre2010][^MassartAbsil2020]. +The metric was used in [JourneeBachAbsilSepulchre:2010](@cite)[MassartAbsil:2020](@cite). # Constructor @@ -39,19 +39,6 @@ The metric was used in[^JourneeBachAbsilSepulchre2010][^MassartAbsil2020]. Generate the manifold of $n × n$ symmetric positive semidefinite matrices of rank $k$ over the `field` of real numbers `ℝ` or complex numbers `ℂ`. - -[^JourneeBachAbsilSepulchre2010]: - > Journée, M., Bach, F., Absil, P.-A., and Sepulchre, R.: - > “Low-Rank Optimization on the Cone of Positive Semidefinite Matrices”, - > SIAM Journal on Optimization (20)5, pp. 2327–2351, 2010. - > doi: [10.1137/080731359](https://doi.org/10.1137/080731359), - > arXiv: [0807.4423](http://arxiv.org/abs/0807.4423). -[^MassartAbsil2020]: - > Massart, E., Absil, P.-A.: - > "Quotient Geometry with Simple Geodesics for the AbstractManifold of Fixed-Rank Positive-Semidefinite Matrices", - > SIAM Journal on Matrix Analysis and Applications (41)1, pp. 171–198, 2020. - > doi: [10.1137/18m1231389](https://doi.org/10.1137/18m1231389), - > preprint: [sites.uclouvain.be/absil/2018.06](https://sites.uclouvain.be/absil/2018.06). """ struct SymmetricPositiveSemidefiniteFixedRank{n,k,𝔽} <: AbstractDecoratorManifold{𝔽} end diff --git a/src/manifolds/Symplectic.jl b/src/manifolds/Symplectic.jl index cb8a11290c..7de358593d 100644 --- a/src/manifolds/Symplectic.jl +++ b/src/manifolds/Symplectic.jl @@ -22,7 +22,7 @@ That is, the symplectic manifold consists of with ``0_n`` and ``I_n`` denoting the ``n × n`` zero-matrix and indentity matrix in ``ℝ^{n \times n}`` respectively. -The tangent space at a point ``p`` is given by [^BendokatZimmermann2021] +The tangent space at a point ``p`` is given by [BendokatZimmermann:2021](@cite) ````math \begin{align*} T_p\operatorname{Sp}(2n) @@ -54,7 +54,7 @@ end RealSymplecticMetric <: RiemannianMetric The canonical Riemannian metric on the symplectic manifold, -defined pointwise for ``p \in \operatorname{Sp}(2n)`` by [^Fiori2011] +defined pointwise for ``p \in \operatorname{Sp}(2n)`` by [Fiori:2011](@cite)] ````math \begin{align*} & g_p \colon T_p\operatorname{Sp}(2n) \times T_p\operatorname{Sp}(2n) \rightarrow ℝ, \\ @@ -264,7 +264,7 @@ ManifoldsBase.default_retraction_method(::Symplectic) = CayleyRetraction() distance(M::Symplectic, p, q) Compute an approximate geodesic distance between two Symplectic matrices -``p, q \in \operatorname{Sp}(2n)``, as done in [^WangSunFiori2018]. +``p, q \in \operatorname{Sp}(2n)``, as done in [WangSunFiori:2018](@cite). ````math \operatorname{dist}(p, q) ≈ ||\operatorname{Log}(p^+q)||_{\operatorname{Fr}}, @@ -317,18 +317,12 @@ The Exponential mapping on the Symplectic manifold with the [`RealSymplecticMetric`](@ref) Riemannian metric. For the point ``p \in \operatorname{Sp}(2n)`` the exponential mapping along the tangent -vector ``X \in T_p\operatorname{Sp}(2n)`` is computed as [^WangSunFiori2018] +vector ``X \in T_p\operatorname{Sp}(2n)`` is computed as [WangSunFiori:2018](@cite) ````math \operatorname{exp}_p(X) = p \operatorname{Exp}((p^{-1}X)^{\mathrm{T}}) \operatorname{Exp}(p^{-1}X - (p^{-1}X)^{\mathrm{T}}), ```` where ``\operatorname{Exp}(\cdot)`` denotes the matrix exponential. - -[^WangSunFiori2018]: - > Wang, Jing and Sun, Huafei and Fiori, Simone: - > A Riemannian-steepest-descent approach for optimization on the real symplectic group, - > Mathematical Methods in the Applied Sciences, 41(11) pp. 4273-4286, 2018 - > doi [10.1002/mma.4890](https://doi.org/10.1002/mma.4890) """ exp(::Symplectic, ::Any...) @@ -352,7 +346,7 @@ Compute the manifold gradient ``\text{grad}f(p)`` of a scalar function The element ``\text{grad}f(p)`` is found as the Riesz representer of the differential ``\text{D}f(p) \colon T_p\operatorname{Sp}(2n) \rightarrow ℝ`` w.r.t. -the Riemannian metric inner product at ``p`` [^Fiori2011]. +the Riemannian metric inner product at ``p`` [Fiori:2011](@cite)]. That is, ``\text{grad}f(p) \in T_p\operatorname{Sp}(2n)`` solves the relation ````math g_p(\text{grad}f(p), X) = \text{D}f(p) \quad\forall\; X \in T_p\operatorname{Sp}(2n). @@ -372,13 +366,6 @@ w.r.t the Riemannian metric ``g_p`` extended to the entire embedding space. If `false`, compute the gradient by first projecting ``∇f(p)`` onto the tangent vector space, before changing the representer in the tangent vector space to comply with the [`RealSymplecticMetric`](@ref). - - -[^Fiori2011]: - > Simone Fiori: - > Solving minimal-distance problems over the manifold of real-symplectic matrices, - > SIAM Journal on Matrix Analysis and Applications 32(3), pp. 938-968, 2011. - > doi [10.1137/100817115](https://doi.org/10.1137/100817115). """ function ManifoldDiff.gradient( M::Symplectic, @@ -503,7 +490,7 @@ end Compute the Cayley Inverse Retraction ``X = \mathcal{L}_p^{\operatorname{Sp}}(q)`` such that the Cayley Retraction from ``p`` along ``X`` lands at ``q``, i.e. -``\mathcal{R}_p(X) = q`` [^BendokatZimmermann2021]. +``\mathcal{R}_p(X) = q`` [BendokatZimmermann:2021](@cite). First, recall the definition the standard symplectic matrix ````math @@ -527,7 +514,7 @@ If that is the case, the inverse cayley retration at ``p`` applied to ``q`` is ∈ T_p\operatorname{Sp}(2n). ```` -[^BendokatZimmermann2021]: +[BendokatZimmermann:2021](@cite): > Bendokat, Thomas and Zimmermann, Ralf: > The real symplectic Stiefel and Grassmann manifolds: metrics, geodesics and applications > arXiv preprint arXiv:[2108.12447](https://arxiv.org/abs/2108.12447), 2021. @@ -603,10 +590,12 @@ end Compute the projection of ``X ∈ R^{2n × 2n}`` onto ``T_p\operatorname{Sp}(2n, ℝ)`` w.r.t. the Riemannian metric ``g`` [`RealSymplecticMetric`](@ref). -The closed form projection mapping is given by [^GaoSonAbsilStykel2021] +The closed form projection mapping is given by [GaoSonAbsilStykel:2021](@cite) + ````math \operatorname{P}^{T_p\operatorname{Sp}(2n)}_{g_p}(X) = pQ\operatorname{sym}(p^{\mathrm{T}}Q^{\mathrm{T}}X), ```` + where ``\operatorname{sym}(A) = \frac{1}{2}(A + A^{\mathrm{T}})``. This function is not exported. """ @@ -632,22 +621,20 @@ Project onto the normal of the tangent space ``(T_p\operatorname{Sp}(2n))^{\perp a point ``p ∈ \operatorname{Sp}(2n)``, relative to the riemannian metric ``g`` [`RealSymplecticMetric`](@ref). That is, + ````math (T_p\operatorname{Sp}(2n))^{\perp_g} = \{Y \in \mathbb{R}^{2n \times 2n} : g_p(Y, X) = 0 \;\forall\; X \in T_p\operatorname{Sp}(2n)\}. ```` -The closed form projection operator onto the normal space is given by [^GaoSonAbsilStykel2021] + +The closed form projection operator onto the normal space is given by [GaoSonAbsilStykel:2021](@cite) + ````math \operatorname{P}^{(T_p\operatorname{Sp}(2n))\perp}_{g_p}(X) = pQ\operatorname{skew}(p^{\mathrm{T}}Q^{\mathrm{T}}X), ```` + where ``\operatorname{skew}(A) = \frac{1}{2}(A - A^{\mathrm{T}})``. This function is not exported. - -[^GaoSonAbsilStykel2021]: - > Gao, Bin and Son, Nguyen Thanh and Absil, P-A and Stykel, Tatjana: - > Riemannian optimization on the symplectic Stiefel manifold, - > SIAM Journal on Optimization 31(2), pp. 1546-1575, 2021. - > doi [10.1137/20M1348522](https://doi.org/10.1137/20M1348522) """ function project_normal!( ::MetricManifold{𝔽,Euclidean{Tuple{m,n},𝔽},ExtendedSymplecticMetric}, @@ -686,7 +673,7 @@ and then symmetrizes it as `S = S + S'`. Then ``S`` is normalized to have Frobenius norm of `hamiltonian_norm` and `X = pQS` is returned, where `Q` is the [`SymplecticMatrix`](@ref). """ -function Base.rand( +function Random.rand( M::Symplectic; vector_at=nothing, hamiltonian_norm=(vector_at === nothing ? 1 / 2 : 1.0), @@ -725,7 +712,7 @@ end Compute the Cayley retraction on ``p ∈ \operatorname{Sp}(2n, ℝ)`` in the direction of tangent vector ``X ∈ T_p\operatorname{Sp}(2n, ℝ)``, -as defined in by Birtea et al in proposition 2 [^BirteaCaşuComănescu2020]. +as defined in by Birtea et al in proposition 2 [BirteaCaşuComănescu:2020](@cite). Using the symplectic inverse of a matrix ``A \in ℝ^{2n \times 2n}``, `` @@ -752,12 +739,6 @@ is defined pointwise as Here ``\operatorname{exp}_{1/1}(z) = (2 - z)^{-1}(2 + z)`` denotes the Padé (1, 1) approximation to ``\operatorname{exp}(z)``. - -[^BirteaCaşuComănescu2020]: - > Birtea, Petre and Caşu, Ioan and Comănescu, Dan: - > Optimization on the real symplectic group, - > Monatshefte f{\"u}r Mathematik, Springer, 2020. - > doi [10.1007/s00605-020-01369-9](https://doi.org/10.1007/s00605-020-01369-9) """ retract(M::Symplectic, p, X) diff --git a/src/manifolds/SymplecticStiefel.jl b/src/manifolds/SymplecticStiefel.jl index 6d5dbebb52..b9e47f487c 100644 --- a/src/manifolds/SymplecticStiefel.jl +++ b/src/manifolds/SymplecticStiefel.jl @@ -16,7 +16,7 @@ Q_{2n} = \end{bmatrix}. ```` -The symplectic Stiefel tangent space at ``p`` can be parametrized as [^BendokatZimmermann2021] +The symplectic Stiefel tangent space at ``p`` can be parametrized as [BendokatZimmermann:2021](@cite) ````math \begin{align*} T_p\operatorname{SpSt}(2n, 2k) @@ -61,7 +61,7 @@ ManifoldsBase.default_retraction_method(::SymplecticStiefel) = CayleyRetraction( Define the canonical projection from ``\operatorname{Sp}(2n, 2n)`` onto ``\operatorname{SpSt}(2n, 2k)``, by projecting onto the first ``k`` columns -and the ``n + 1``'th onto the ``n + k``'th columns [^BendokatZimmermann2021]. +and the ``n + 1``'th onto the ``n + k``'th columns [BendokatZimmermann:2021](@cite). It is assumed that the point ``p`` is on ``\operatorname{Sp}(2n, 2n)``. """ @@ -125,7 +125,7 @@ $A^{+} = Q_{2k}^{\mathrm{T}}A^{\mathrm{T}}Q_{2n}$ is the symplectic inverse, wit ```` The we check that ``H = p^{+}X \in 𝔤_{2k}``, where ``𝔤`` is the Lie Algebra of the symplectic group ``\operatorname{Sp}(2k)``, -characterized as [^BendokatZimmermann2021], +characterized as [BendokatZimmermann:2021](@cite), ````math 𝔤_{2k} = \{H \in ℝ^{2k \times 2k} \;|\; H^+ = -H \}. ```` @@ -164,7 +164,7 @@ at a point ``p \in \operatorname{SpSt}(2n, 2k)`` in the direction of ``X \in T_p\operatorname{SpSt}(2n, 2k)``. The tangent vector ``X`` can be written in the form -``X = \bar{\Omega}p`` [^BendokatZimmermann2021], with +``X = \bar{\Omega}p`` [BendokatZimmermann:2021](@cite), with ````math \bar{\Omega} = X (p^{\mathrm{T}}p)^{-1}p^{\mathrm{T}} + Q_{2n}p(p^{\mathrm{T}}p)^{-1}X^{\mathrm{T}}(I_{2n} - Q_{2n}^{\mathrm{T}}p(p^{\mathrm{T}}p)^{-1}p^{\mathrm{T}}Q_{2n})Q_{2n} @@ -180,7 +180,7 @@ where ``\operatorname{Exp}(\cdot)`` denotes the matrix exponential. Computing the above mapping directly however, requires taking matrix exponentials of two ``2n \times 2n`` matrices, which is computationally expensive when ``n`` -increases. Therefore we instead follow [^BendokatZimmermann2021] who express the above +increases. Therefore we instead follow [BendokatZimmermann:2021](@cite) who express the above exponential mapping in a way which only requires taking matrix exponentials of an ``8k \times 8k`` matrix and a ``4k \times 4k`` matrix. @@ -290,7 +290,7 @@ get_total_space(::SymplecticStiefel{n,k,𝔽}) where {n,k,𝔽} = Symplectic{n, Compute the Riemannian inner product ``g^{\operatorname{SpSt}}`` at ``p \in \operatorname{SpSt}`` between tangent vectors ``X, X \in T_p\operatorname{SpSt}``. -Given by Proposition 3.10 in [^BendokatZimmermann2021]. +Given by Proposition 3.10 in [BendokatZimmermann:2021](@cite). ````math g^{\operatorname{SpSt}}_p(X, Y) = \operatorname{tr}\left(X^{\mathrm{T}}\left(I_{2n} - @@ -375,7 +375,7 @@ end Compute the Cayley Inverse Retraction ``X = \mathcal{L}_p^{\operatorname{SpSt}}(q)`` such that the Cayley Retraction from ``p`` along ``X`` lands at ``q``, i.e. -``\mathcal{R}_p(X) = q`` [^BendokatZimmermann2021]. +``\mathcal{R}_p(X) = q`` [BendokatZimmermann:2021](@cite). First, recall the definition the standard symplectic matrix ````math @@ -400,11 +400,6 @@ If that is the case, the inverse cayley retration at ``p`` applied to ``q`` is \mathcal{L}_p^{\operatorname{Sp}}(q) = 2p\bigl(V - U\bigr) + 2\bigl((p + q)U - p\bigr) ∈ T_p\operatorname{Sp}(2n). ```` - -[^BendokatZimmermann2021]: - > Bendokat, Thomas and Zimmermann, Ralf: - > The real symplectic Stiefel and Grassmann manifolds: metrics, geodesics and applications - > arXiv preprint arXiv:[2108.12447](https://arxiv.org/abs/2108.12447), 2021. """ inverse_retract(::SymplecticStiefel, p, q, ::CayleyInverseRetraction) @@ -427,7 +422,7 @@ is_flat(M::SymplecticStiefel) = false manifold_dimension(::SymplecticStiefel{n, k}) Returns the dimension of the symplectic Stiefel manifold embedded in ``ℝ^{2n \times 2k}``, -i.e. [^BendokatZimmermann2021] +i.e. [BendokatZimmermann:2021](@cite) ````math \operatorname{dim}(\operatorname{SpSt}(2n, 2k)) = (4n - 2k + 1)k. ```` @@ -495,7 +490,7 @@ To generate random tangent vectors at ``p`` then, this function sets ``B_X = 0`` and generates a random Hamiltonian matrix ``Ω_X \in \mathfrak{sp}(2n,F)`` with Frobenius norm of `hamiltonian_norm` before returning ``X = pΩ_X``. """ -function Base.rand( +function Random.rand( M::SymplecticStiefel{n}; vector_at=nothing, hamiltonian_norm=(vector_at === nothing ? 1 / 2 : 1.0), @@ -530,7 +525,7 @@ Given a point ``p \in \operatorname{SpSt}(2n, 2k)``, every tangent vector \tilde{\Omega} = \left(I_{2n} - \frac{1}{2}pp^+\right)Xp^+ - pX^+\left(I_{2n} - \frac{1}{2}pp^+\right) \in ℝ^{2n \times 2n}, ```` -as shown in Proposition 3.5 of [^BendokatZimmermann2021]. +as shown in Proposition 3.5 of [BendokatZimmermann:2021](@cite). Using this representation of ``X``, the Cayley retraction on ``\operatorname{SpSt}(2n, 2k)`` is defined pointwise as ````math @@ -540,7 +535,7 @@ The operator ``\operatorname{cay}(A) = (I - A)^{-1}(I + A)`` is the Cayley trans However, the computation of an ``2n \times 2n`` matrix inverse in the expression above can be reduced down to inverting a ``2k \times 2k`` matrix due to Proposition -5.2 of [^BendokatZimmermann2021]. +5.2 of [BendokatZimmermann:2021](@cite). Let ``A = p^+X`` and ``H = X - pA``. Then an equivalent expression for the Cayley retraction defined pointwise above is diff --git a/src/manifolds/Tucker.jl b/src/manifolds/Tucker.jl index 24c91cb03a..8ae708c871 100644 --- a/src/manifolds/Tucker.jl +++ b/src/manifolds/Tucker.jl @@ -10,7 +10,7 @@ Segre manifold, i.e., the set of rank-1 tensors. Let ``\mathbb{F}`` be the real or complex numbers. Any tensor ``p`` on the Tucker manifold can be represented as a multilinear product in HOSVD -[^DeLathauwer2000] form +[DeLathauwerDeMoorVanderwalle:2000](@cite) form ```math p = (U_1,\dots,U_D) \cdot \mathcal{C} ``` @@ -21,7 +21,7 @@ the matrix ``U_d \in \mathbb{F}^{N_d \times R_d}`` contains the singular vectors # Tangent space The tangent space to the Tucker manifold at -``p = (U_1,\dots,U_D) \cdot \mathcal{C}`` is [^Koch2010] +``p = (U_1,\dots,U_D) \cdot \mathcal{C}`` is [KochLubich:2010](@cite) ```math T_p \mathcal{M} = \bigl\{ @@ -40,16 +40,6 @@ where ``\mathcal{C}^\prime`` is arbitrary, ``U_d^{\mathrm{H}}`` is the Hermitian Generate the manifold of `field`-valued tensors of dimensions `N[1] × … × N[D]` and multilinear rank `R = (R[1], …, R[D])`. - -[^DeLathauwer2000]: - > Lieven De Lathauwer, Bart De Moor, Joos Vandewalle: "A multilinear singular value decomposition" - > SIAM Journal on Matrix Analysis and Applications, 21(4), pp. 1253-1278, 2000 - > doi: [10.1137/S0895479896305696](https://doi.org/10.1137/S0895479896305696) - -[^Koch2010]: - > Othmar Koch, Christian Lubic, "Dynamical Tensor approximation" - > SIAM Journal on Matrix Analysis and Applications, 31(5), pp. 2360-2375, 2010 - > doi: [10.1137/09076578X](https://doi.org/10.1137/09076578X) """ struct Tucker{N,R,D,𝔽} <: AbstractManifold{𝔽} end function Tucker(n⃗::NTuple{D,Int}, r⃗::NTuple{D,Int}, field::AbstractNumbers=ℝ) where {D} @@ -92,13 +82,7 @@ assumptions are made. The low-multilinear rank tensor arising from the sequentially truncated the higher-order singular value decomposition of the `D`-dimensional array `p` of type `T`. The singular values are truncated to get a multilinear rank `mlrank` -[^Vannieuwenhoven2012]. - -[^Vannieuwenhoven2012]: - > Nick Vannieuwenhoven, Raf Vandebril, Karl Meerbergen: "A new truncation strategy for the higher-order singular value decomposition" - > SIAM Journal on Scientific Computing, 34(2), pp. 1027-1052, 2012 - > doi: [10.1137/110836067](https://doi.org/10.1137/110836067) - +[VannieuwenhovenVanderbrilMeerbergen:2012](@cite). """ struct TuckerPoint{T,D} <: AbstractManifoldPoint hosvd::HOSVD{T,D} @@ -322,7 +306,7 @@ end Base.convert(::Type{Matrix{T}}, basis::CachedBasis{𝔽,DefaultOrthonormalBasis{𝔽, TangentSpaceType},HOSVDBasis{T, D}}) where {𝔽, T, D} Base.convert(::Type{Matrix}, basis::CachedBasis{𝔽,DefaultOrthonormalBasis{𝔽, TangentSpaceType},HOSVDBasis{T, D}}) where {𝔽, T, D} -Convert a HOSVD-derived cached basis from [^Dewaele2021] of the `D`th order +Convert a HOSVD-derived cached basis from [DewaeleBreidingVannieuwenhoven:2021](@cite) of the `D`th order [`Tucker`](@ref) manifold with number type `T` to a matrix. The columns of this matrix are the vectorisations of the [`embed`](@ref)dings of the basis vectors. @@ -454,7 +438,7 @@ An implicitly stored basis of the tangent space to the Tucker manifold. Assume ``p = (U_1,\dots,U_D) \cdot \mathcal{C}`` is in HOSVD format and that, for ``d=1,\dots,D``, the singular values of the ``d``'th unfolding are ``\sigma_{dj}``, with ``j = 1,\dots,R_d``. -The basis of the tangent space is as follows: [^Dewaele2021] +The basis of the tangent space is as follows: [DewaeleBreidingVannieuwenhoven:2021](@cite) ````math \bigl\{ @@ -467,10 +451,6 @@ The basis of the tangent space is as follows: [^Dewaele2021] for all ``d = 1,\dots,D`` and all canonical basis vectors ``e_i`` and ``e_j``. Every ``U_d^\perp`` is such that ``[U_d \quad U_d^{\perp}]`` forms an orthonormal basis of ``\mathbb{R}^{N_d}``. - -[^Dewaele2021]: - > Nick Dewaele, Paul Breiding, Nick Vannieuwenhoven, "The condition number of many tensor decompositions is invariant under Tucker compression" - > arxiv: [2106.13034](https://arxiv.org/abs/2106.13034) """ function get_basis( ::Tucker, @@ -681,17 +661,11 @@ end @doc raw""" retract(::Tucker, p::TuckerPoint, X::TuckerTVector, ::PolarRetraction) -The truncated HOSVD-based retraction [^Kressner2014] to the Tucker manifold, i.e. +The truncated HOSVD-based retraction [KressnerSteinlechnerVandereycken:2013](@cite) to the Tucker manifold, i.e. the result is the sequentially tuncated HOSVD approximation of ``p + X``. In the exceptional case that the multilinear rank of ``p + X`` is lower than that of ``p``, this retraction produces a boundary point, which is outside the manifold. - -[^Kressner2014]: - > Daniel Kressner, Michael Steinlechner, Bart Vandereycken: "Low-rank tensor completion by Riemannian optimization" - > BIT Numerical Mathematics, 54(2), pp. 447-468, 2014 - > doi: [10.1007/s10543-013-0455-z](https://doi.org/10.1007/s10543-013-0455-z) - """ retract(::Tucker, ::Any, ::Any, ::PolarRetraction) diff --git a/src/manifolds/Unitary.jl b/src/manifolds/Unitary.jl index 2dc220c724..7345f90143 100644 --- a/src/manifolds/Unitary.jl +++ b/src/manifolds/Unitary.jl @@ -116,4 +116,41 @@ end show(io::IO, ::UnitaryMatrices{n,ℂ}) where {n} = print(io, "UnitaryMatrices($(n))") show(io::IO, ::UnitaryMatrices{n,ℍ}) where {n} = print(io, "UnitaryMatrices($(n), ℍ)") +@doc raw""" + riemannian_Hessian(M::UnitaryMatrices, p, G, H, X) + +The Riemannian Hessian can be computed by adopting Eq. (5.6) [Nguyen:2023](@cite), +so very similar to the complex Stiefel manifold. +The only difference is, that here the tangent vectors are stored +in the Lie algebra, i.e. the update direction is actually ``pX`` instead of just ``X`` (in Stiefel). +and that means the inverse has to be appliead to the (Euclidean) Hessian +to map it into the Lie algebra. +""" +riemannian_Hessian(M::UnitaryMatrices, p, G, H, X) +function riemannian_Hessian!(M::UnitaryMatrices, Y, p, G, H, X) + symmetrize!(Y, G' * p) + project!(M, Y, p, p' * H - X * Y) + return Y +end + +@doc raw""" + Weingarten(M::UnitaryMatrices, p, X, V) + +Compute the Weingarten map ``\mathcal W_p`` at `p` on the [`Stiefel`](@ref) `M` with respect to the +tangent vector ``X \in T_p\mathcal M`` and the normal vector ``V \in N_p\mathcal M``. + +The formula is due to [AbsilMahonyTrumpf:2013](@cite) given by + +```math +\mathcal W_p(X,V) = -\frac{1}{2}p\bigl(V^{\mathrm{H}}X - X^\mathrm{H}V\bigr) +``` +""" +Weingarten(::UnitaryMatrices, p, X, V) + +function Weingarten!(::UnitaryMatrices, Y, p, X, V) + Y .= V' * X + Y .= -p * 1 / 2 * (Y - Y') + return Y +end + Manifolds.zero_vector(::UnitaryMatrices{1,ℍ}, p) = zero(p) diff --git a/src/manifolds/VectorBundle.jl b/src/manifolds/VectorBundle.jl index 774822760c..64a0a993c6 100644 --- a/src/manifolds/VectorBundle.jl +++ b/src/manifolds/VectorBundle.jl @@ -85,7 +85,7 @@ TangentSpaceAtPoint(M::AbstractManifold, p) = VectorSpaceAtPoint(TangentBundleFi """ TangentSpace(M::AbstractManifold, p) -Return a [`TangentSpaceAtPoint`](@ref Manifolds.TangentSpaceAtPoint) representing tangent space at `p` on the [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold) `M`. +Return a [`TangentSpaceAtPoint`](@ref) representing tangent space at `p` on the [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold) `M`. """ TangentSpace(M::AbstractManifold, p) = VectorSpaceAtPoint(TangentBundleFibers(M), p) @@ -106,7 +106,7 @@ end struct SasakiRetraction <: AbstractRetractionMethod end Exponential map on [`TangentBundle`](@ref) computed via Euler integration as described -in [^Muralidharan2012]. The system of equations for $\gamma : ℝ \to T\mathcal M$ such that +in [MuralidharanFlecther:2012](@cite). The system of equations for $\gamma : ℝ \to T\mathcal M$ such that $\gamma(1) = \exp_{p,X}(X_M, X_F)$ and $\gamma(0)=(p, X)$ reads ```math @@ -120,11 +120,6 @@ where $R$ is the Riemann curvature tensor (see [`riemann_tensor`](@ref)). SasakiRetraction(L::Int) In this constructor `L` is the number of integration steps. - -[^Muralidharan2012]: - > P. Muralidharan and P. T. Fletcher, “Sasaki Metrics for Analysis of Longitudinal Data - > on Manifolds,” Proc IEEE Comput Soc Conf Comput Vis Pattern Recognit, vol. 2012, - > pp. 1027–1034, Jun. 2012, doi: [10.1109/CVPR.2012.6247780](https://doi.org/10.1109/CVPR.2012.6247780). """ struct SasakiRetraction <: AbstractRetractionMethod L::Int @@ -209,16 +204,12 @@ end """ TangentBundle{𝔽,M} = VectorBundle{𝔽,TangentSpaceType,M} where {𝔽,M<:AbstractManifold{𝔽}} -Tangent bundle for manifold of type `M`, as a manifold with the Sasaki metric [^Sasaki1958]. +Tangent bundle for manifold of type `M`, as a manifold with the Sasaki metric [Sasaki:1958](@cite). Exact retraction and inverse retraction can be approximated using [`VectorBundleProductRetraction`](@ref), [`VectorBundleInverseProductRetraction`](@ref) and [`SasakiRetraction`](@ref). [`VectorBundleProductVectorTransport`](@ref) can be used as a vector transport. -[^Sasaki1958]: - > S. Sasaki, “On the differential geometry of tangent bundles of Riemannian manifolds,” - > Tohoku Math. J. (2), vol. 10, no. 3, pp. 338–354, 1958, doi: [10.2748/tmj/1178244668](https://doi.org/10.2748/tmj/1178244668). - # Constructors TangentBundle(M::AbstractManifold) @@ -603,7 +594,7 @@ end @doc raw""" injectivity_radius(M::TangentSpaceAtPoint) -Return the injectivity radius on the [`TangentSpaceAtPoint`](@ref Manifolds.TangentSpaceAtPoint) `M`, which is $∞$. +Return the injectivity radius on the [`TangentSpaceAtPoint`](@ref) `M`, which is $∞$. """ injectivity_radius(::TangentSpaceAtPoint) = Inf @@ -726,7 +717,7 @@ end """ is_flat(::TangentSpaceAtPoint) -Return true. [`TangentSpaceAtPoint`](@ref Manifolds.TangentSpaceAtPoint) is a flat manifold. +Return true. [`TangentSpaceAtPoint`](@ref) is a flat manifold. """ is_flat(::TangentSpaceAtPoint) = true """ @@ -1266,6 +1257,19 @@ function zero_vector!(B::TangentBundleFibers, X, p) return zero_vector!(B.manifold, X, p) end +@doc raw""" + Y = Weingarten(M::VectorSpaceAtPoint, p, X, V) + Weingarten!(M::VectorSpaceAtPoint, Y, p, X, V) + +Compute the Weingarten map ``\mathcal W_p`` at `p` on the [`VectorSpaceAtPoint`](@ref) `M` with respect to the +tangent vector ``X \in T_p\mathcal M`` and the normal vector ``V \in N_p\mathcal M``. + +Since this a flat space by itself, the result is always the zero tangent vector. +""" +Weingarten(::VectorSpaceAtPoint, p, X, V) + +Weingarten!(::VectorSpaceAtPoint, Y, p, X, V) = fill!(Y, 0) + @doc raw""" zero_vector(B::VectorBundle, p) diff --git a/src/statistics.jl b/src/statistics.jl index 717498ecc4..7a1030e60f 100644 --- a/src/statistics.jl +++ b/src/statistics.jl @@ -67,52 +67,16 @@ points are in an open geodesic ball about the mean with corresponding radius (see [`GeodesicInterpolationWithinRadius`](@ref)): * All simply connected complete Riemannian manifolds with non-positive sectional - curvature at radius $∞$ [^Cheng2016], in particular: + curvature at radius $∞$ [ChengHoSalehianVemuri:2016](@cite), in particular: + [`Euclidean`](@ref) - + [`SymmetricPositiveDefinite`](@ref) [^Ho2013] + + [`SymmetricPositiveDefinite`](@ref) [HoChengSalehianVemuri:2013](@cite) * Other manifolds: - + [`Sphere`](@ref): $\frac{π}{2}$ [^Salehian2015] - + [`Grassmann`](@ref): $\frac{π}{4}$ [^Chakraborty2015] - + [`Stiefel`](@ref)/[`Rotations`](@ref): $\frac{π}{2 \sqrt 2}$ [^Chakraborty2019] + + [`Sphere`](@ref): $\frac{π}{2}$ [SalehianEtAl:2015](@cite) + + [`Grassmann`](@ref): $\frac{π}{4}$ [ChakrabortyVemuri:2015](@cite) + + [`Stiefel`](@ref)/[`Rotations`](@ref): $\frac{π}{2 \sqrt 2}$ [ChakrabortyVemuri:2019](@cite) For online variance computation, the algorithm additionally uses an analogous -recursion to the weighted Welford algorithm [^West1979]. - -[^Ho2013]: - > Ho J.; Cheng G.; Salehian H.; Vemuri B. C.; Recursive Karcher expectation - > estimators and geometric law of large numbers. - > Proceedings of the 16th International Conference on Artificial Intelligence - > and Statistics (2013), pp. 325–332. - > [pdf](http://proceedings.mlr.press/v31/ho13a.pdf). -[^Salehian2015]: - > Salehian H.; Chakraborty R.; Ofori E.; Vaillancourt D.; An efficient - > recursive estimator of the Fréchet mean on a hypersphere with applications - > to Medical Image Analysis. - > Mathematical Foundations of Computational Anatomy (2015). - > [pdf](https://www-sop.inria.fr/asclepios/events/MFCA15/Papers/MFCA15_4_2.pdf). -[^Chakraborty2015]: - > Chakraborty R.; Vemuri B. C.; Recursive Fréchet Mean Computation on the - > Grassmannian and Its Applications to Computer Vision. - > Proceedings of the IEEE International Conference on Computer Vision (ICCV) (2015), - > pp. 4229-4237. - > doi: [10.1109/ICCV.2015.481](https://doi.org/10.1109/ICCV.2015.481), - > [link](http://openaccess.thecvf.com/content_iccv_2015/html/Chakraborty_Recursive_Frechet_Mean_ICCV_2015_paper.html). -[^Cheng2016]: - > Cheng G.; Ho J.; Salehian H.; Vemuri B. C.; Recursive Computation of the - > Fréchet Mean on Non-positively Curved Riemannian Manifolds with Applications. - > Riemannian Computing in Computer Vision. Springer, Cham (2016), pp. 21-43. - > doi: [10.1007/978-3-319-22957-7_2](https://doi.org/10.1007/978-3-319-22957-7_2), - > [pdf](https://www.cise.ufl.edu/~vemuri/paperphp/article.php?y=2016&i=5). -[^Chakraborty2019]: - > Chakraborty R.; Vemuri B. C.; Statistics on the (compact) Stiefel manifold: - > Theory and Applications. - > The Annals of Statistics (2019), 47(1), pp. 415-438. - > doi: [10.1214/18-AOS1692](https://doi.org/10.1214/18-AOS1692), - > arxiv: [1708.00045](https://arxiv.org/abs/1708.00045). -[^West1979]: - > West D. H. D.; Updating Mean and Variance Estimates: An Improved Method. - > Communications of the ACM (1979), 22(9), pp. 532–535. - > doi: [10.1145/359146.359153](https://doi.org/10.1145/359146.359153). +recursion to the weighted Welford algorithm [West:1979](@cite). """ struct GeodesicInterpolation <: AbstractEstimationMethod end @@ -189,7 +153,7 @@ end Estimate the covariance matrix of a set of points `x` on manifold `M`. Since the covariance matrix on a manifold is a rank 2 tensor, the function returns its coefficients in basis induced by -the given tangent space basis. See Section 5 of [^Pennec2006] for details. +the given tangent space basis. See Section 5 of [Pennec:2006](@cite) for details. The mean is calculated using the specified `mean_estimation_method` using [mean](@ref Statistics.mean(::AbstractManifold, ::AbstractVector, ::AbstractEstimationMethod), @@ -199,11 +163,6 @@ Finally, the covariance matrix in the tangent plane is estimated using the Eucli in [`StatsBase.jl`](https://juliastats.org/StatsBase.jl/stable/cov/#StatsBase.CovarianceEstimator) and examples of covariance estimation methods can be found in [`CovarianceEstimation.jl`](https://github.com/mateuszbaran/CovarianceEstimation.jl/). - -[^Pennec2006]: - > X. Pennec, “Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric - > Measurements,” J Math Imaging Vis, vol. 25, no. 1, p. 127, Jul. 2006, - > doi: [10.1007/s10851-006-6228-4](https:///doi.org/10.1007/s10851-006-6228-4). """ function Statistics.cov( M::AbstractManifold, @@ -279,26 +238,9 @@ between two iterates is small. For more stopping criteria check the Optionally, pass `retraction` and `inverse_retraction` method types to specify the (inverse) retraction. -The Theory stems from[^Karcher1977] and is also described in[^PennecArsigny2013] +The Theory stems from [Karcher:1977](@cite) and is also described in [PennecArsigny:2012](@cite) as the exponential barycenter. -The algorithm is further described in[^Afsari2013]. - -[^Afsari2013]: - > Afsari, B; Tron, R.; Vidal, R.: On the Convergence of Gradient - > Descent for Finding the Riemannian Center of Mass, - > SIAM Journal on Control and Optimization (2013), 51(3), pp. 2230–2260, - > doi: [10.1137/12086282X](https://doi.org/10.1137/12086282X), - > arxiv: [1201.0925](https://arxiv.org/abs/1201.0925) -[^PennecArsigny2013]: - > Pennec X., Arsigny V.: Exponential Barycenters of the Canonical Cartan Connection and - > Invariant Means on Lie Groups. - > In: Nielsen F., Bhatia R. (eds) Matrix Information Geometry, (2013), pp. 123-166. - > doi: [10.1007/978-3-642-30232-9_7](https://doi.org/10.1007/978-3-642-30232-9_7), - > hal: [https://hal.inria.fr/hal-00699361/document](https://hal.inria.fr/hal-00699361/document) -[^Karcher1977]: - > Karcher, H.: Riemannian center of mass and mollifier smoothing. - > Communications on Pure Applied Mathematics (1977), 30, pp. 509–541. - > doi [10.1002/cpa.3160300502](https://doi.org/10.1002/cpa.3160300502) +The algorithm is further described in[AfsariTronVidal:2013](@cite). """ mean(::AbstractManifold, ::Any...) function Statistics.mean( @@ -646,13 +588,9 @@ change between two iterates is small. For more stopping criteria check the Optionally, pass `retraction` and `inverse_retraction` method types to specify the (inverse) retraction. -The algorithm is further described in [^Bačák2014]. +The algorithm is further described in [Bacak:2014](@cite). + -[^Bačák2014]: - > Bačák, M: Computing Medians and Means in Hadamard Spaces. - > SIAM Journal on Optimization (2014), 24(3), pp. 1542–1566, - > doi: [10.1137/140953393](https://doi.org/10.1137/140953393), - > arxiv: [1210.2145](https://arxiv.org/abs/1210.2145) """ Statistics.median( ::AbstractManifold, @@ -710,7 +648,7 @@ change between two iterates is small. For more stopping criteria check the The parameter ``α\in (0,2]`` is a step size. -The algorithm is further described in [^FletcherVenkatasubramanianJoshi2008], +The algorithm is further described in [FletcherVenkatasubramanianJoshi:2008](@cite), especially the update rule in Eq. (6), i.e. Let ``q_{k}`` denote the current iterate, $n$ the number of points ``x_1,\ldots,x_n``, and @@ -734,12 +672,6 @@ and where $\mathrm{d}_{\mathcal M}$ denotes the Riemannian [`distance`](@ref). Optionally, pass `retraction` and `inverse_retraction` method types to specify the (inverse) retraction, which by default use the exponential and logarithmic map, respectively. - -[^FletcherVenkatasubramanianJoshi2008]: - > Fletcher, T., Venkatasubramanian, S., Joshi, S: - > Robust statistics on Riemannian manifolds via the geometric median - > 2008 IEEE Conference on Computer Vision and Pattern Recognition, - > doi: [10.1109/CVPR.2008.4587747](https://doi.org/10.1109/CVPR.2008.4587747), """ Statistics.median( ::AbstractManifold, diff --git a/src/utils.jl b/src/utils.jl index 2956778da4..d934ce8c9e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -221,6 +221,25 @@ For example `select_from_tuple(("a", "b", "c"), Val((3, 1, 1)))` returns return Expr(:tuple, [Expr(:ref, :t, k) for k in P]...) end +@doc raw""" + symmetrize!(Y, X) + +Given a quare matrix `X` compute `1/2 .* (X' + X)` in place of `Y` +""" +function symmetrize!(Y, X) + Y .= (X' .+ X) ./ 2 + return Y +end + +@doc raw""" + symmetrize(X) + +Given a quare matrix `X` compute `1/2 .* (X' + X)`. +""" +function symmetrize(X) + return (X' .+ X) ./ 2 +end + @doc raw""" vec2skew!(X, v, k) diff --git a/test/manifolds/euclidean.jl b/test/manifolds/euclidean.jl index eb071d4403..c4024d6850 100644 --- a/test/manifolds/euclidean.jl +++ b/test/manifolds/euclidean.jl @@ -366,6 +366,15 @@ using Manifolds: induced_basis ) === 2.0 end + @testset "Weingarten & Hessian" begin + M = Euclidean(2) + p = [1.0, 2.0] + G = [3.0, 4.0] + H = [5.0, 6.0] + X = [7.0, 8.0] + rH = riemannian_Hessian(M, p, G, H, X) + @test rH == H + end @testset "Volume" begin @test manifold_volume(Euclidean(2)) == Inf @test volume_density(E, p, X) == 1.0 diff --git a/test/manifolds/fixed_rank.jl b/test/manifolds/fixed_rank.jl index c3424830aa..e381bb4da8 100644 --- a/test/manifolds/fixed_rank.jl +++ b/test/manifolds/fixed_rank.jl @@ -230,4 +230,13 @@ include("../utils.jl") ) end end + @testset "Riemannian Hessian" begin + M = FixedRankMatrices(3, 2, 2) + pE = [1.0 0.0; 0.0 1.0; 0.0 0.0] + p = SVDMPoint(pE) + X = UMVTVector([0.0 0.0; 0.0 0.0; 1.0 1.0], [1.0 0.0; 0.0 1.0], [1.0 0.0; 0.0 1.0]) + G = [1.0 0.0; 0.0 2.0; 0.0 0.0] + H = [0.0 3.0; 0.0 4.0; 0.0 1.0] + @test is_vector(M, p, riemannian_Hessian(M, p, G, H, X)) + end end diff --git a/test/manifolds/grassmann.jl b/test/manifolds/grassmann.jl index ef91560650..bbab887d7b 100644 --- a/test/manifolds/grassmann.jl +++ b/test/manifolds/grassmann.jl @@ -380,4 +380,12 @@ include("../utils.jl") end end end + @testset "Riemannian Hessian" begin + M = Grassmann(3, 2) + p = [1.0 0.0; 0.0 1.0; 0.0 0.0] + X = [0.0 0.0; 0.0 0.0; 0.0 1.0] + Y = [0.0 1.0; -1.0 0.0; 1.0 0.0] + Z = [0.0 -1.0; 1.0 0.0; 1.0 0.0] + @test riemannian_Hessian(M, p, Y, Z, X) == [0.0 0.0; 0.0 0.0; 2.0 0.0] + end end diff --git a/test/manifolds/hyperbolic.jl b/test/manifolds/hyperbolic.jl index 734e091ccd..0fede44efe 100644 --- a/test/manifolds/hyperbolic.jl +++ b/test/manifolds/hyperbolic.jl @@ -339,6 +339,16 @@ include("../utils.jl") @test dpX[2][1] == -1.0 @test dpX[2][2].X == 2 * normalize(X) end + @testset "Riemannian Hessian" begin + M = Hyperbolic(2) + p = [0.0, 0.0, 1.0] + G = [1.0, 0.2, 0.3] + H = [2.0, 0.3, 0.4] + X = [0.3, 0.4, 0.0] + D = diagm([1.0, 1.0, -1.0]) + rH = project(M, p, D * H + dot(p, D * G) .* X) + @test riemannian_Hessian(M, p, G, H, X) == rH + end @testset "Manifold volume" begin M = Hyperbolic(2) @test manifold_volume(M) == Inf diff --git a/test/manifolds/positive_numbers.jl b/test/manifolds/positive_numbers.jl index c86215d834..6670b03f14 100644 --- a/test/manifolds/positive_numbers.jl +++ b/test/manifolds/positive_numbers.jl @@ -77,6 +77,15 @@ include("../utils.jl") ) end end + @testset "Weingarten & Hessian" begin + M = PositiveNumbers() + p = 1.0 + G = 2.0 + H = 3.0 + X = 4.0 + rH = riemannian_Hessian(M, p, G, H, X) + @test rH == p * H * p + X * G * p + end @testset "Manifold volume" begin M5 = PositiveVectors(3) diff --git a/test/manifolds/power_manifold.jl b/test/manifolds/power_manifold.jl index 541f04830e..dbce456621 100644 --- a/test/manifolds/power_manifold.jl +++ b/test/manifolds/power_manifold.jl @@ -423,6 +423,28 @@ end @test norm(N, P, Z .- Zc) ≈ 0 end + @testset "Hessian conversion" begin + M = Sphere(2) + N = PowerManifold(M, 2) + p = [1.0 0.0; 0.0 0.0; 1.0 0.0] + q = 1 / sqrt(2) * [1.0 0.0; 1.0 1.0; 0.0 1.0] + q = 1 / sqrt(2) * [0.0 1.0; 1.0 1.0; 1.0 0.0] + r = 1 / sqrt(3) * [1.0 1.0; 1.0 1.0; 1.0 1.0] + X = log(M, p, q) + Y = log(M, p, r) + Z = -X + H1 = riemannian_Hessian(N, p, Y, Z, X) + H2 = cat( + [riemannian_Hessian(M, p[:, i], Y[:, i], Z[:, i], X[:, i]) for i in 1:2]..., + dims=2, + ) + @test H1 == H2 + V = [0.2 0.0; 0.0 0.0; 0.0 0.3] + W1 = Weingarten(N, p, X, V) + W2 = cat([Weingarten(M, p[:, i], X[:, i], V[:, i]) for i in 1:2]..., dims=2) + @test W1 == W2 + end + @testset "Nested replacing RNG" begin M = PowerManifold(Ms, NestedReplacingPowerRepresentation(), 2) @test is_point(M, rand(M)) diff --git a/test/manifolds/product_manifold.jl b/test/manifolds/product_manifold.jl index 237239250b..de0e96ebf7 100644 --- a/test/manifolds/product_manifold.jl +++ b/test/manifolds/product_manifold.jl @@ -756,6 +756,26 @@ using RecursiveArrayTools: ArrayPartition ArrayPartition([0.3535533905932738, -0.35355339059327373, 0.0], [1.0, 1.5]) end + @testset "Hessian conversion" begin + M = Sphere(2) + N = M × M + p = ArrayPartition([1.0, 0.0, 0.0], [0.0, 0.0, 1.0]) + q = 1 / sqrt(2) * ArrayPartition([1.0, 1.0, 0.0], [0.0, 1.0, 1.0]) + q = 1 / sqrt(2) * ArrayPartition([0.0, 1.0, 1.0], [1.0, 1.0, 0.0]) + r = 1 / sqrt(3) * ArrayPartition([1.0, 1.0, 1.0], [1.0, 1.0, 1.0]) + X = log(M, p, q) + Y = log(M, p, r) + Z = -X + H1 = riemannian_Hessian(N, p, Y, Z, X) + H2 = ArrayPartition( + [riemannian_Hessian(M, p[i, :], Y[i, :], Z[i, :], X[i, :]) for i in 1:2]..., + ) + @test H1 == H2 + V = ArrayPartition([0.2, 0.0, 0.0], [0.0, 0.0, 0.3]) + W1 = Weingarten(N, p, X, V) + W2 = ArrayPartition([Weingarten(M, p[i, :], X[i, :], V[i, :]) for i in 1:2]...) + @test W1 == W2 + end @testset "Manifold volume" begin MS2 = Sphere(2) MS3 = Sphere(3) diff --git a/test/manifolds/rotations.jl b/test/manifolds/rotations.jl index 5a6ca9f38b..a729cfa62a 100644 --- a/test/manifolds/rotations.jl +++ b/test/manifolds/rotations.jl @@ -235,6 +235,16 @@ include("../utils.jl") @test injectivity_radius(M, PolarRetraction()) == 0.0 @test injectivity_radius(M, p, PolarRetraction()) == 0.0 end + @testset "Riemannian Hessian" begin + M = Rotations(2) + p = Matrix{Float64}(I, 2, 2) + X = [0.0 3.0; -3.0 0.0] + V = [1.0 0.0; 1.0 0.0] + @test Weingarten(M, p, X, V) == -1 / 2 * p * (V' * X - X' * V) + G = [0.0 1.0; 0.0 0.0] + H = [0.0 0.0; 2.0 0.0] + @test riemannian_Hessian(M, p, G, H, X) == [0.0 -1.0; 1.0 0.0] + end @testset "riemann_tensor" begin M = Rotations(3) diff --git a/test/manifolds/sphere.jl b/test/manifolds/sphere.jl index 61cbe0d013..98e2b955a8 100644 --- a/test/manifolds/sphere.jl +++ b/test/manifolds/sphere.jl @@ -264,6 +264,14 @@ using ManifoldsBase: TFVector @test dpX[2][1] == 1.0 @test dpX[2][2].X == normalize(X) end + @testset "Weingarten" begin + M = Sphere(2) + p = [1.0, 0.0, 0.0] + X = [0.0, 2.0, 0.0] + V = [0.3, 0.0, 0.0] + Y = -X * (p'V) + @test Weingarten(M, p, X, V) == Y + end @testset "Volume density" begin M = Sphere(2) diff --git a/test/manifolds/stiefel.jl b/test/manifolds/stiefel.jl index edd765dc5e..2daa966b0b 100644 --- a/test/manifolds/stiefel.jl +++ b/test/manifolds/stiefel.jl @@ -60,6 +60,10 @@ include("../utils.jl") @test Y == X Z = change_representer(M, EuclideanMetric(), p, X) @test Z == X + # In this case it stays as is + @test riemannian_Hessian(M, p, Y, Z, X) == Z + V = [1.0 1.0; 0.0 0.0; 0.0 0.0] # From T\bot_pM. + @test Weingarten(M, p, X, V) == [0.0 0.0; 0.0 0.0; 1.0 1.0] end types = [Matrix{Float64}] TEST_FLOAT32 && push!(types, Matrix{Float32}) @@ -331,6 +335,8 @@ include("../utils.jl") @test inner(M3, p, X, Y) == 0 @test inner(M3, p, X, 2 * X + 3 * Y) == 2 * inner(M3, p, X, X) @test norm(M3, p, X) ≈ distance(M3, p, q) + Z = [0.0 0.0; 0.0 0.0; -1.0 -1.0] + @test riemannian_Hessian(M3, p, Y, Z, X) == [0.0 0.5; -0.5 0.0; -1.0 -1.0] # check on a higher dimensional manifold, that the iterations are actually used M4 = MetricManifold(Stiefel(10, 2), CanonicalMetric()) p = Matrix{Float64}(I, 10, 2) @@ -511,5 +517,28 @@ include("../utils.jl") @test isapprox(MM, p, log(MM, p, p), zero_vector(MM, p); atol=1e-6) end end + + @testset "Hessian Conversion" begin + M1 = MetricManifold(Stiefel(3, 2), StiefelSubmersionMetric(-0.5)) + M2 = Stiefel(3, 2) + M2b = MetricManifold(Stiefel(3, 2), EuclideanMetric()) + M3 = MetricManifold(Stiefel(3, 2), StiefelSubmersionMetric(0.0)) + M4 = MetricManifold(Stiefel(3, 2), CanonicalMetric()) + p = [1.0 0.0; 0.0 1.0; 0.0 0.0] + X = [0.0 0.0; 0.0 0.0; 1.0 1.0] + Y = [0.0 0.0; 0.0 0.0; -1.0 1.0] + Z = [0.0 0.0; 0.0 0.0; -1.0 -1.0] + rH = riemannian_Hessian(M2, p, Y, Z, X) + @test riemannian_Hessian(M1, p, Y, Z, X) == rH #Special case of submersion metric + @test riemannian_Hessian(M2b, p, Y, Z, X) == rH # metric is default + @test riemannian_Hessian(M3, p, Y, Z, X) == riemannian_Hessian(M4, p, Y, Z, X) + + V = [0.0 -1.0; 1.0 0.0; 0.0 0.0] + W = zero_vector(M2, p) + Weingarten!(M2, W, p, X, V) + Wb = zero_vector(M2b, p) + Weingarten!(M2b, Wb, p, X, V) + @test W == Wb + end end end diff --git a/test/manifolds/symmetric_positive_definite.jl b/test/manifolds/symmetric_positive_definite.jl index 27cb3f7d3d..0dc0a1854f 100644 --- a/test/manifolds/symmetric_positive_definite.jl +++ b/test/manifolds/symmetric_positive_definite.jl @@ -277,6 +277,17 @@ include("../utils.jl") @test is_point(M, p1) end + @testset "Riemannian Hessian Conversion" begin + M = SymmetricPositiveDefinite(2) + p = [2.0 1.0; 1.0 1.0] + G = [1.0 0.0; 1.0 0.1] + H = [2.0 1.0; 0.0 0.0] + X = [0.0 3.0; 3.0 0.0] + Y = X * 1 / 2 * (G' + G) * p + Y = 1 / 2 * p * (H' + H) * p + 1 / 2 * (Y' + Y) + @test riemannian_Hessian(M, p, G, H, X) == Y + end + @testset "Volume density" begin M = SymmetricPositiveDefinite(3) @test manifold_volume(M) == Inf diff --git a/test/manifolds/unitary_matrices.jl b/test/manifolds/unitary_matrices.jl index 05888dfd26..3dee3b3d3e 100644 --- a/test/manifolds/unitary_matrices.jl +++ b/test/manifolds/unitary_matrices.jl @@ -54,6 +54,16 @@ end @test isapprox(M, p, X, X2) r1 = exp(M, p, X, 1.0) @test isapprox(M, r, r1; atol=1e-10) + + @testset "Riemannian Hessian" begin + p = Matrix{Float64}(I, 2, 2) + X = [0.0 3.0; -3.0 0.0] + V = [1.0 0.0; 1.0 0.0] + @test Weingarten(M, p, X, V) == -1 / 2 * p * (V' * X - X' * V) + G = [0.0 1.0; 0.0 0.0] + H = [0.0 0.0; 2.0 0.0] + @test riemannian_Hessian(M, p, G, H, X) == [0.0 -1.0; 1.0 0.0] + end end @testset "Special unitary matrices" begin diff --git a/test/manifolds/vector_bundle.jl b/test/manifolds/vector_bundle.jl index b8d8018fd2..89472ee238 100644 --- a/test/manifolds/vector_bundle.jl +++ b/test/manifolds/vector_bundle.jl @@ -271,4 +271,13 @@ struct TestVectorSpaceType <: VectorSpaceType end @test is_flat(M) @test injectivity_radius(M) == Inf end + + @testset "Weingarten Map" begin + p0 = [1.0, 0.0, 0.0] + M = TangentSpaceAtPoint(Sphere(2), p0) + p = [0.0, 1.0, 1.0] + X = [0.0, 1.0, 0.0] + V = [1.0, 0.0, 0.0] + @test Weingarten(M, p, X, V) == zero_vector(M, p) + end end diff --git a/tutorials/Project.toml b/tutorials/Project.toml index ffcf36e7be..51df6bcbef 100644 --- a/tutorials/Project.toml +++ b/tutorials/Project.toml @@ -19,6 +19,6 @@ CSV = "0.10" DataFrames = "1" Distributions = "0.22.6, 0.23, 0.24, 0.25" IJulia = "1" -Manifolds = "0.8.46" -ManifoldsBase = "0.14.5" +Manifolds = "0.8.71" +ManifoldsBase = "0.14.10" MultivariateStats = "0.10" diff --git a/tutorials/getstarted.qmd b/tutorials/getstarted.qmd index 3bab7d6005..2e70ea10a1 100644 --- a/tutorials/getstarted.qmd +++ b/tutorials/getstarted.qmd @@ -250,7 +250,7 @@ In this section we take a look how to implement generic functions on manifolds. For our example here, we want to implement the so-called [📖 Bézier curve](https://en.wikipedia.org /wiki/Bézier_curve) using the so-called [📖 de-Casteljau algorithm](https://en.wikipedia.org/wiki/De_Casteljau%27s_algorithm). -The linked algorithm can easily be generalised to manifolds by replacing lines with geodesics. This was for example used in [^BergmannGousenbourger2018] and the following example is an extended version of an example from [^AxenBaranBergmannRzecki2022]. +The linked algorithm can easily be generalised to manifolds by replacing lines with geodesics. This was for example used in [BergmannGousenbourger:2018](@cite) and the following example is an extended version of an example from [AxenBaranBergmannRzecki:2023](@cite). The algorithm works recursively. For the case that we have a Bézier curve with just two points, the algorithm just evaluates the geodesic connecting both at some time point $t∈[0,1]$. The function to evaluate a shortest geodesic (it might not be unique, but then a deterministic choice is taken) between two points `p` and `q` on a manifold `M` [🔗 `shortest_geodesic(M, p, q, t)`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions.html#ManifoldsBase.shortest_geodesic-Tuple{AbstractManifold,%20Any,%20Any}). @@ -271,7 +271,7 @@ end Which can now be used on any manifold where the shortest geodesic is implemented -Now on several manifolds the [📖 exponential map](https://en.wikipedia.org/wiki/Exponential_map_(Riemannian_geometry)) and its (locally defined) inverse, the logarithmic map might not be available in an implementation. So one way to generalise this, is the use of a retraction (see [^AbsilMahonySepulchre2008], Def. 4.1.1 for details) and its (local) inverse. +Now on several manifolds the [📖 exponential map](https://en.wikipedia.org/wiki/Exponential_map_(Riemannian_geometry)) and its (locally defined) inverse, the logarithmic map might not be available in an implementation. So one way to generalise this, is the use of a retraction (see [AbsilMahonySepulchre:2008](@cite), Def. 4.1.1 for details) and its (local) inverse. The function itself is quite similar to the expponential map, just that [🔗 `retract(M, p, X, m)`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.retract) has one further parameter, the type of retraction to take, so `m` is a subtype of [`AbstractRetractionMethod`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.AbstractRetractionMethod) `m`, the same for the [🔗 `inverse_retract(M, p, q, n)`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.inverse_retract) with an [`AbstractInverseRetractionMethod`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/retractions.html#ManifoldsBase.AbstractInverseRetractionMethod) `n`. @@ -424,22 +424,11 @@ The same approach is used for Again, for all of these, the concrete types only have to be used if you want to do a second, different from the details, property, for example a second way to embed a manfiold. If a manifold is (in its usual representation) an embedded manifold, this works with the default manifold type already, since then it is again set as the reasonable default. -## References - -[^AbsilMahonySepulchre2008]: - > Absil, P.-A., Mahony, R. and Sepulchre R., - > _Optimization Algorithms on Matrix Manifolds_ - > Princeton University Press, 2008, - > doi: [10.1515/9781400830244](https://doi.org/10.1515/9781400830244) - > [open access](http://press.princeton.edu/chapters/absil/) -[^AxenBaranBergmannRzecki2022]: - >Axen, S. D., Baran, M., Bergmann, R. and Rzecki, K: - > _Manifolds.jl: An Extensible Julia Framework for Data Analysis on Manifolds_, - > arXiv preprint, 2022, [2106.08777](https://arxiv.org/abs/2106.08777) -[^BergmannGousenbourger2018]: - > Bergmann, R. and Gousenbourger, P.-Y.: - > _A variational model for data fitting on manifolds - > by minimizing the acceleration of a Bézier curve_. - > Frontiers in Applied Mathematics and Statistics, 2018. - > doi: [10.3389/fams.2018.00059](https://dx.doi.org/10.3389/fams.2018.00059), - > arXiv: [1807.10090](https://arxiv.org/abs/1807.10090) \ No newline at end of file +## Literature + +````{=commonmark} +```@bibliography +Pages = ["tutorials/getstarted.md"] +Canonical=false +``` +```` \ No newline at end of file diff --git a/tutorials/integration.qmd b/tutorials/integration.qmd index 4f51f6631d..019ac52199 100644 --- a/tutorials/integration.qmd +++ b/tutorials/integration.qmd @@ -24,7 +24,7 @@ In principle, there are many ways in which a manifold can be equipped with a mea One of the most popular ways is based on pushing the Lebesgue measure on a tangent space through the exponential map. Any other suitable atlas could be used, not just the one defined by normal coordinates, though each one requires different volume density corrections due to the Jacobian determinant of the pushforward. `Manifolds.jl` provides the function [`volume_density`](@ref) that calculates that quantity, denoted $\theta_p(X)$. -See for example [^leBrigantPuechmorel2019], Definition 11, for a precise description using Jacobi fields. +See for example [LeBrignantPuechmorel:2019](@cite), Definition 11, for a precise description using Jacobi fields. While many sources define volume density as a function of two points, `Manifolds.jl` decided to use the more general point-tangent vector formulation. The two-points variant can be implemented as ```{julia} @@ -91,7 +91,7 @@ Note that our `pdf_manifold` implements a simplified version of $f_{X_M}$ which In such case there is only one non-zero summand in the formula for $f_{X_M}(q)$, namely $X=\log_p(q)$. Otherwise we would have to consider other vectors $Y\in T_p \mathcal{M}$ such that $\exp_p(Y) = q$ in that sum. -Remarkably, exponential-wrapped distributions possess three important qualities [^ChevallierLiLuDunson2022]: +Remarkably, exponential-wrapped distributions possess three important qualities [ChevallierLiLuDunson:2022](@cite): * Densities of $X_M$ are explicit. There is no normalization constant that needs to be computed like in truncated distributions. * Sampling from $X_M$ is easy. It suffices to get a sample from $X_T$ and pass it to the exponential map. @@ -107,7 +107,7 @@ f(q) = \frac{1}{n h^d} \sum_{i=1}^n \frac{1}{\theta_q(\log_q(p_i))}K\left( \frac where $h$ is the bandwidth, a small positive number less than the injectivity radius of $\mathcal M$ and $K\colon\mathbb{R}\to\mathbb{R}$ is a kernel function. Note that Pelletier's estimator can only use radially-symmetric kernels. -The radially symmetric multivariate Epanechnikov kernel used in the example below is described in [^LangrenéWarin2019]. +The radially symmetric multivariate Epanechnikov kernel used in the example below is described in [LangreneMarin:2019](@cite). ```{julia} struct PelletierKDE{TM<:AbstractManifold,TPts<:AbstractVector} @@ -159,7 +159,7 @@ However, once a constant is picked (and it must be picked before any useful comp ### Haar measures On Lie groups the situation regarding integration is more complicated. -Invariance under left or right group action is a desired property that leads one to consider Haar measures [^Tornier2020]. +Invariance under left or right group action is a desired property that leads one to consider Haar measures [Tornier:2020](@cite). It is, however, unclear what are the practical benefits of considering Haar measures over the Lebesgue measure of the underlying manifold, which often turns out to be invariant anyway. ### Integration in charts @@ -167,32 +167,14 @@ It is, however, unclear what are the practical benefits of considering Haar meas Integration through charts is an approach currently not supported by `Manifolds.jl`. One has to define a suitable set of disjoint charts covering the entire manifold and use a method for multivariate Euclidean integration. Note that ranges of parameters have to be adjusted for each manifold and scaling based on the metric needs to be applied. -See [^BoyaSudarshanTilma2003] for some considerations on symmetric spaces. +See [BoyaSudarshanTilma:2003](@cite) for some considerations on symmetric spaces. ## References -[^Tornier2020]: - > S. Tornier, “Haar Measures.” arXiv, Jun. 19, 2020. - > doi: [10.48550/arXiv.2006.10956](https://doi.org/10.48550/arXiv.2006.10956). - -[^LangrenéWarin2019]: - > N. Langrené and X. Warin, “Fast and Stable Multivariate Kernel Density Estimation by - > Fast Sum Updating,” Journal of Computational and Graphical Statistics, vol. 28, no. 3, - > pp. 596–608, Jul. 2019, - > doi: [10.1080/10618600.2018.1549052](https://doi.org/10.1080/10618600.2018.1549052). - -[^leBrigantPuechmorel2019]: - > A. le Brigant and S. Puechmorel, “Approximation of Densities on Riemannian Manifolds,” - > Entropy, vol. 21, no. 1, Art. no. 1, Jan. 2019, - > doi: [10.3390/e21010043](https://doi.org/10.3390/e21010043). - -[^ChevallierLiLuDunson2022]: - > E. Chevallier, D. Li, Y. Lu, and D. B. Dunson, “Exponential-wrapped distributions on - > symmetric spaces.” arXiv, Oct. 09, 2022. - > doi: [10.48550/arXiv.2009.01983](https://doi.org/10.48550/arXiv.2009.01983). - -[^BoyaSudarshanTilma2003]: - > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports - > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003, - > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1) +## Literature + +```@bibliography +Pages = ["tutorials/integration.md"] +Canonical=false +``` \ No newline at end of file