Introduction of the orthographic retraction for the fixed rank manifo…

…ld (#662) * Add orthographic retraction and inverse retraction. * Add boilerplate code for orthographic retractions and their inverses and already 2 new references. * add test to new retraction method and fixes typo. * add docs to retract_orthographic and inverse_retract_orthographic * Adds Lukas to Project members. * runs formatter * extend tests to in-place inverse retractions. * Bump version, adapt and fix News.md --------- Co-authored-by: Ronny Bergmann <[email protected]> Co-authored-by: Mateusz Baran <[email protected]>
JuliaManifolds · Oct 25, 2023 · b0ab892 · b0ab892 · kellertuer · Oct 25, 2023
1 parent 3536123
commit b0ab892
Show file tree

Hide file tree

Showing 8 changed files with 193 additions and 7 deletions.
diff --git a/.zenodo.json b/.zenodo.json
@@ -31,6 +31,11 @@
             "name": "Weiß, Manuel",
             "type": "ProjectMember"
         },
+        {
+            "affiliation": "Georg-August-University Göttingen",
+            "name": "Klingbiel, Lukas",
+            "type": "ProjectMember"
+        },
         {
             "affiliation": "NTNU Trondheim",
             "name": "Kolstø, Johannes Voll",

diff --git a/NEWS.md b/NEWS.md
@@ -5,7 +5,13 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [0.9.0] - 2023-mm-dd
+## [0.9.1] - 2023-10-25
+
+### Added
+
+- a new retraction and its inverse for the fixed Rank Manifolds, the orthographic retraction.
+
+## [0.9.0] - 2023-10-24
 
 ### Added
 
@@ -57,7 +63,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   - `SymplecticStiefel{n,k}`,
   - `TranslationGroup`,
   - `Tucker`.
-  
+
   For example
 
   ```{julia}
@@ -84,7 +90,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   ```
 
   for groups with size stored in field. Alternatively, you can use a single generic method like this:
-  
+
   ```{julia}
   function Base.show(io::IO, M::CenteredMatrices{T}) where {T}
       m, n = get_parameter(M)

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>", "Antoine Levitt <[email protected]>"]
-version = "0.9.0"
+version = "0.9.1"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

diff --git a/docs/src/references.bib b/docs/src/references.bib
@@ -25,7 +25,28 @@ @incollection{AbsilMahonyTrumpf:2013
     PUBLISHER = {Springer Berlin Heidelberg},
     PAGES     = {361--368},
     TITLE     = {An Extrinsic Look at the Riemannian Hessian},
-    BOOKTITLE = {Geometric Science of Information},
+    BOOKTITLE = {Geometric Science of Information}
+}
+@article{AbsilMalick:2012,
+    AUTHOR   = {Absil, P.-A. and Malick, Jérôme},
+    TITLE    = {Projection-like Retractions on Matrix Manifolds},
+    JOURNAL  = {SIAM Journal on Optimization},
+    VOLUME   = {22},
+    NUMBER   = {1},
+    PAGES    = {135-158},
+    YEAR     = {2012},
+    DOI      = {10.1137/100802529}
+}
+@article{AbsilOseledets:2014,
+    DOI       = {10.1007/s10589-014-9714-4},
+    YEAR      = {2014},
+    PUBLISHER = {Springer Science and Business Media LLC},
+    VOLUME    = {62},
+    NUMBER    = {1},
+    PAGES     = {5--29},
+    AUTHOR    = {P.-A. Absil and I. V. Oseledets},
+    TITLE     = {Low-rank retractions: a survey and new results},
+    JOURNAL   = {Computational Optimization and Applications}
 }
 @article{AfsariTronVidal:2013,
     DOI        = {10.1137/12086282x},

diff --git a/ext/ManifoldsTestExt/tests_general.jl b/ext/ManifoldsTestExt/tests_general.jl
@@ -333,21 +333,36 @@ function test_manifold(
                         Test.@test isapprox(M, p2, q; atol=point_atol)
                         # This test is not reasonable for `inverse_retract!(M, X, p, q, m)`,
                         # since X is of different type/concept than p,q
+
                     end
             end
         end
         for p in pts
             epsx = find_eps(p)
             for inv_retr_method in inverse_retraction_methods
+                X = inverse_retract(M, p, p, inv_retr_method)
                 Test.@test isapprox(
                     M,
                     p,
                     zero_vector(M, p),
-                    inverse_retract(M, p, p, inv_retr_method);
+                    X;
                     atol=epsx * retraction_atol_multiplier,
                     rtol=retraction_atol_multiplier == 0 ?
                          sqrt(epsx) * retraction_rtol_multiplier : 0,
                 )
+                if (test_inplace && is_mutating)
+                    Y = copy(M, p, X)
+                    inverse_retract!(M, Y, p, p, inv_retr_method)
+                    Test.@test isapprox(
+                        M,
+                        p,
+                        zero_vector(M, p),
+                        Y;
+                        atol=epsx * retraction_atol_multiplier,
+                        rtol=retraction_atol_multiplier == 0 ?
+                             sqrt(epsx) * retraction_rtol_multiplier : 0,
+                    )
+                end
             end
         end
     end

diff --git a/src/Manifolds.jl b/src/Manifolds.jl
@@ -132,6 +132,8 @@ import ManifoldsBase:
     representation_size,
     retract,
     retract!,
+    _retract,
+    retract!,
     retract_cayley!,
     retract_exp_ode!,
     retract_pade!,
@@ -729,6 +731,7 @@ export AbstractRetractionMethod,
     ProjectionRetraction,
     SoftmaxRetraction,
     ODEExponentialRetraction,
+    OrthographicRetraction,
     PadeRetraction,
     ProductRetraction,
     SasakiRetraction
@@ -739,6 +742,7 @@ export AbstractInverseRetractionMethod,
     CayleyInverseRetraction,
     LogarithmicInverseRetraction,
     QRInverseRetraction,
+    OrthographicInverseRetraction,
     PolarInverseRetraction,
     ProjectionInverseRetraction,
     ShootingInverseRetraction,

diff --git a/src/manifolds/FixedRankMatrices.jl b/src/manifolds/FixedRankMatrices.jl
@@ -118,6 +118,67 @@ Base.:-(v::UMVTVector) = UMVTVector(-v.U, -v.M, -v.Vt)
 Base.:+(v::UMVTVector) = UMVTVector(v.U, v.M, v.Vt)
 Base.:(==)(v::UMVTVector, w::UMVTVector) = (v.U == w.U) && (v.M == w.M) && (v.Vt == w.Vt)
 
+# Move to Base when name is established – i.e. used in more than one manifold
+# |/---
+"""
+    OrthographicRetraction <: AbstractRetractionMethod
+
+Retractions that are related to orthographic projections, which was first
+used in [AbsilMalick:2012](@cite).
+"""
+struct OrthographicRetraction <: AbstractRetractionMethod end
+
+"""
+    OrthographicInverseRetraction <: AbstractInverseRetractionMethod
+
+Retractions that are related to orthographic projections, which was first
+used in [AbsilMalick:2012](@cite).
+"""
+struct OrthographicInverseRetraction <: AbstractInverseRetractionMethod end
+
+# Layer II
+function _inverse_retract!(
+    M::AbstractManifold,
+    X,
+    p,
+    q,
+    ::OrthographicInverseRetraction;
+    kwargs...,
+)
+    return inverse_retract_orthographic!(M, X, p, q; kwargs...)
+end
+
+# Layer III
+"""
+    inverse_retract_orthographic!(M::AbstractManifold, X, p, q)
+
+Compute the in-place variant of the [`OrthographicInverseRetraction`](@ref).
+"""
+inverse_retract_orthographic!(M::AbstractManifold, X, p, q)
+
+## Layer II
+function _retract!(
+    M::AbstractManifold,
+    q,
+    p,
+    X,
+    t::Number,
+    ::OrthographicRetraction;
+    kwargs...,
+)
+    return retract_orthographic!(M, q, p, X, t; kwargs...)
+end
+
+## Layer III
+"""
+    retract_orthographic!(M::AbstractManifold, q, p, X, t::Number)
+
+Compute the in-place variant of the [`OrthographicRetraction`](@ref).
+"""
+retract_orthographic!(M::AbstractManifold, q, p, X, t::Number)
+
+# \|---
+
 allocate(p::SVDMPoint) = SVDMPoint(allocate(p.U), allocate(p.S), allocate(p.Vt))
 function allocate(p::SVDMPoint, ::Type{T}) where {T}
     return SVDMPoint(allocate(p.U, T), allocate(p.S, T), allocate(p.Vt, T))
@@ -127,6 +188,9 @@ function allocate(X::UMVTVector, ::Type{T}) where {T}
     return UMVTVector(allocate(X.U, T), allocate(X.M, T), allocate(X.Vt, T))
 end
 
+function allocate_result(M::FixedRankMatrices, ::typeof(inverse_retract), p, q)
+    return zero_vector(M, p)
+end
 function allocate_result(M::FixedRankMatrices, ::typeof(project), X, p, vals...)
     m, n, k = get_parameter(M.size)
     # vals are p and X, so we can use their fields to set up those of the UMVTVector
@@ -378,6 +442,31 @@ function inner(::FixedRankMatrices, x::SVDMPoint, v::UMVTVector, w::UMVTVector)
     return dot(v.U, w.U) + dot(v.M, w.M) + dot(v.Vt, w.Vt)
 end
 
+@doc raw"""
+    inverse_retract(M, p, q, ::OrthographicInverseRetraction)
+
+Compute the orthographic inverse retraction [`FixedRankMatrices`](@ref) `M` by computing
+
+```math
+    X = P_{T_{p}M}(q - p) = qVV^\mathrm{T} + UU^{\mathrm{T}}q - UU^{\mathrm{T}}qVV^{\mathrm{T}} - p,
+```
+where ``p`` is a [`SVDMPoint`](@ref)`(U,S,Vt)` and ``P_{T_{p}M}`` is the [`project`](@ref)ion
+onto the tangent space at ``p``.
+
+For more details, see [AbsilOseledets:2014](@cite).
+"""
+inverse_retract(::FixedRankMatrices, ::Any, ::Any, ::OrthographicInverseRetraction)
+
+function inverse_retract_orthographic!(
+    M::FixedRankMatrices,
+    X::UMVTVector,
+    p::SVDMPoint,
+    q::SVDMPoint,
+)
+    project!(M, X, p, embed(M, q) - embed(M, p))
+    return X
+end
+
 function _isapprox(::FixedRankMatrices, p::SVDMPoint, q::SVDMPoint; kwargs...)
     return isapprox(p.U * Diagonal(p.S) * p.Vt, q.U * Diagonal(q.S) * q.Vt; kwargs...)
 end
@@ -523,6 +612,51 @@ function representation_size(M::FixedRankMatrices)
     return (m, n)
 end
 
+@doc raw"""
+    retract(M::FixedRankMatrices, p, X, ::OrthographicRetraction)
+
+Compute the OrthographicRetraction on the [`FixedRankMatrices`](@ref) `M` by finding
+the nearest point to ``p + X`` in
+
+```math
+    p + X + N_{p}\mathcal M \cap \mathcal M
+```
+
+where ``N_{p}\mathcal M `` is the Normal Space of ``T_{p}\mathcal M ``.
+
+If ``X`` is sufficiently small, then the nearest such point is unique and can be expressed by
+
+```math
+    q = (U(S + M) + U_{p})(S + M)^{-1}((S + M)V^{\mathrm{T}} + V^{\mathrm{T}}_{p}),
+```
+
+where ``p`` is a [`SVDMPoint`](@ref)`(U,S,Vt)` and ``X`` is an [`UMVTVector`](@ref)`(Up,M,Vtp)`.
+
+For more details, see [AbsilOseledets:2014](@cite).
+"""
+retract(::FixedRankMatrices, ::Any, ::Any, ::OrthographicRetraction)
+
+function retract_orthographic!(
+    M::FixedRankMatrices,
+    q::SVDMPoint,
+    p::SVDMPoint,
+    X::UMVTVector,
+    t::Number,
+)
+    m, n, k = get_parameter(M.size)
+    tX = t * X
+    QU, RU = qr(p.U * (diagm(p.S) + tX.M) + tX.U)
+    QV, RV = qr(p.Vt' * (diagm(p.S) + tX.M') + tX.Vt')
+
+    Uk, Sk, Vtk = svd(RU * inv(diagm(p.S) + tX.M) * RV')
+
+    mul!(q.U, QU[:, 1:k], Uk)
+    q.S .= Sk[1:k]
+    mul!(q.Vt, Vtk, QV[:, 1:k]')
+
+    return q
+end
+
 @doc raw"""
     retract(M, p, X, ::PolarRetraction)
 

diff --git a/test/manifolds/fixed_rank.jl b/test/manifolds/fixed_rank.jl
@@ -237,7 +237,8 @@ include("../utils.jl")
                 test_vee_hat=false,
                 test_tangent_vector_broadcasting=true,
                 projection_atol_multiplier=15,
-                retraction_methods=[PolarRetraction()],
+                retraction_methods=[PolarRetraction(), OrthographicRetraction()],
+                inverse_retraction_methods=[OrthographicInverseRetraction()],
                 vector_transport_methods=[ProjectionTransport()],
                 vector_transport_retractions=[PolarRetraction()],
                 vector_transport_inverse_retractions=[PolarInverseRetraction()],