From b0ab892c7f911ad044925f13a3d98f3b812722dc Mon Sep 17 00:00:00 2001 From: Klingbiel <104769091+Klingbiel@users.noreply.github.com> Date: Wed, 25 Oct 2023 10:50:26 +0200 Subject: [PATCH] Introduction of the orthographic retraction for the fixed rank manifold (#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 Co-authored-by: Mateusz Baran --- .zenodo.json | 5 + NEWS.md | 12 ++- Project.toml | 2 +- docs/src/references.bib | 23 ++++- ext/ManifoldsTestExt/tests_general.jl | 17 +++- src/Manifolds.jl | 4 + src/manifolds/FixedRankMatrices.jl | 134 ++++++++++++++++++++++++++ test/manifolds/fixed_rank.jl | 3 +- 8 files changed, 193 insertions(+), 7 deletions(-) diff --git a/.zenodo.json b/.zenodo.json index d4a28f3ea3..4076fbb329 100644 --- 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 index eedd76e525..84ca00e2f7 100644 --- 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 index 8c5a5e0c7f..8f2d4c8e94 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Manifolds" uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann ", "Antoine Levitt "] -version = "0.9.0" +version = "0.9.1" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" diff --git a/docs/src/references.bib b/docs/src/references.bib index 46e3f3d19b..f313545904 100644 --- 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 index 435f1a7484..54b37d7cd4 100644 --- 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 index e56361fcac..2c9a5315cb 100644 --- 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 index 9b09232b1c..97bc66fd18 100644 --- 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 index f779c06444..ef128f8e85 100644 --- 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()],