From b3995393379399ac7fc502098ec3d8d7e0f5978a Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 6 Sep 2023 14:02:16 +0200 Subject: [PATCH] Stiefel, SPD, minor fixes --- NEWS.md | 31 ++++- src/manifolds/CholeskySpace.jl | 34 +++-- src/manifolds/Grassmann.jl | 2 +- src/manifolds/GrassmannStiefel.jl | 2 +- src/manifolds/Rotations.jl | 3 +- src/manifolds/Stiefel.jl | 119 +++++++++++++----- src/manifolds/StiefelCanonicalMetric.jl | 61 +++++---- src/manifolds/StiefelEuclideanMetric.jl | 30 ++--- src/manifolds/StiefelSubmersionMetric.jl | 79 ++++++------ src/manifolds/SymmetricPositiveDefinite.jl | 45 ++++--- ...ymmetricPositiveDefiniteAffineInvariant.jl | 54 ++++---- ...mmetricPositiveDefiniteBuresWasserstein.jl | 8 +- ...tiveDefiniteGeneralizedBuresWasserstein.jl | 24 ++-- .../SymmetricPositiveDefiniteLogCholesky.jl | 58 ++++----- .../SymmetricPositiveDefiniteLogEuclidean.jl | 8 +- test/groups/special_orthogonal.jl | 6 +- test/groups/translation_group.jl | 4 +- test/manifolds/stiefel.jl | 8 +- 18 files changed, 325 insertions(+), 251 deletions(-) diff --git a/NEWS.md b/NEWS.md index a2829afd23..e7ba48b346 100644 --- a/NEWS.md +++ b/NEWS.md @@ -18,7 +18,20 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Sizes of all manifolds can now be either encoded in type or stored in a field to avoid over-specialization. The default is set to store the size in a field. To obtain the old behavior, pass the `parameter=:type` keyword argument to manifold constructor. Related changes: - - Statically sized `SpecialEuclidean{N}` is now `SpecialEuclidean{TypeParameter{Tuple{N}}}`, whereas the type of special Euclidean group with field-stored size is `SpecialEuclidean{Tuple{Int}}`. Similar change applies to `GeneralUnitaryMultiplicationGroup{n}`, `Orthogonal{n}`, `SpecialOrthogonal{n}`, `SpecialUnitary{n}`, `SpecialEuclideanManifold{n}`, `TranslationGroup`. For example + - Statically sized `SpecialEuclidean{N}` is now `SpecialEuclidean{TypeParameter{Tuple{N}}}`, whereas the type of special Euclidean group with field-stored size is `SpecialEuclidean{Tuple{Int}}`. Similar change applies to: + - `CholeskySpace{N}`, + - `Euclidean`, + - `GeneralUnitaryMultiplicationGroup{n}`, + - `Grassmann{n,k}`, + - `Orthogonal{n}`, + - `SpecialOrthogonal{n}`, + - `SpecialUnitary{n}`, + - `SpecialEuclideanManifold{n}`, + - `Stiefel{n,k}`, + - `SymmetricPositiveDefinite{n}`, + - `TranslationGroup`. + + For example ```{julia} function Base.show(io::IO, ::SpecialEuclidean{n}) where {n} @@ -43,10 +56,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 end ``` - for groups with size stored in field. + for groups with size stored in field. Alternatively, you can use a single generic method like this: + + ```{julia} + function Base.show(io::IO, G::SpecialEuclidean{T}) where {T} + n = get_n(G) + if T <: TypeParameter + return print(io, "SpecialEuclidean($(n); parameter=:type)") + else + return print(io, "SpecialEuclidean($(n))") + end + end + ``` + - Argument order for type alias `RotationActionOnVector`: most often dispatched on argument is now first. ### Removed - `ProductRepr` is removed; please use `ArrayPartition` instead. -- Default methods throwing "not implemented" `ErrorException` for some group-related operations. \ No newline at end of file +- Default methods throwing "not implemented" `ErrorException` for some group-related operations. diff --git a/src/manifolds/CholeskySpace.jl b/src/manifolds/CholeskySpace.jl index 40f8cbd56f..1c713aceaf 100644 --- a/src/manifolds/CholeskySpace.jl +++ b/src/manifolds/CholeskySpace.jl @@ -1,5 +1,5 @@ @doc raw""" - CholeskySpace{N} <: AbstractManifold{ℝ} + CholeskySpace{T} <: AbstractManifold{ℝ} The manifold of lower triangular matrices with positive diagonal and a metric based on the cholesky decomposition. The formulae for this manifold @@ -7,13 +7,18 @@ are for example summarized in Table 1 of [Lin:2019](@cite). # Constructor - CholeskySpace(n) + CholeskySpace(n; parameter::Symbol=:field) Generate the manifold of $n× n$ lower triangular matrices with positive diagonal. """ -struct CholeskySpace{N} <: AbstractManifold{ℝ} end +struct CholeskySpace{T} <: AbstractManifold{ℝ} + size::T +end -CholeskySpace(n::Int) = CholeskySpace{n}() +function CholeskySpace(n::Int; parameter::Symbol=:field) + size = wrap_type_parameter(parameter, (n,)) + return CholeskySpace{typeof(size)}(size) +end @doc raw""" check_point(M::CholeskySpace, p; kwargs...) @@ -105,6 +110,9 @@ function exp!(::CholeskySpace, q, p, X) return q end +get_n(::CholeskySpace{TypeParameter{N}}) where {N} = N +get_n(M::CholeskySpace{Tuple{Int}}) = get_parameter(M.size)[1] + @doc raw""" inner(M::CholeskySpace, p, X, Y) @@ -164,16 +172,28 @@ Return the manifold dimension for the [`CholeskySpace`](@ref) `M`, i.e. \dim(\mathcal M) = \frac{N(N+1)}{2}. ```` """ -@generated manifold_dimension(::CholeskySpace{N}) where {N} = div(N * (N + 1), 2) +function manifold_dimension(M::CholeskySpace) + N = get_n(M) + return div(N * (N + 1), 2) +end @doc raw""" representation_size(M::CholeskySpace) Return the representation size for the [`CholeskySpace`](@ref)`{N}` `M`, i.e. `(N,N)`. """ -@generated representation_size(::CholeskySpace{N}) where {N} = (N, N) +function representation_size(M::CholeskySpace) + N = get_n(M) + return (N, N) +end -Base.show(io::IO, ::CholeskySpace{N}) where {N} = print(io, "CholeskySpace($(N))") +function Base.show(io::IO, ::CholeskySpace{TypeParameter{Tuple{n}}}) where {n} + return print(io, "CholeskySpace($(n); parameter=:type)") +end +function Base.show(io::IO, M::CholeskySpace{Tuple{Int}}) + n = get_n(M) + return print(io, "CholeskySpace($(n))") +end # two small helpers for strictly lower and upper triangulars strictlyLowerTriangular(p) = LowerTriangular(p) - Diagonal(diag(p)) diff --git a/src/manifolds/Grassmann.jl b/src/manifolds/Grassmann.jl index ed315f3e30..f989907758 100644 --- a/src/manifolds/Grassmann.jl +++ b/src/manifolds/Grassmann.jl @@ -65,7 +65,7 @@ A good overview can be found in[BendokatZimmermannAbsil:2020](@cite). # Constructor - Grassmann(n,k,field=ℝ) + Grassmann(n, k, field=ℝ, parameter::Symbol=:field) Generate the Grassmann manifold $\operatorname{Gr}(n,k)$, where the real-valued case `field = ℝ` is the default. diff --git a/src/manifolds/GrassmannStiefel.jl b/src/manifolds/GrassmannStiefel.jl index 6877de3f86..035ad99cd3 100644 --- a/src/manifolds/GrassmannStiefel.jl +++ b/src/manifolds/GrassmannStiefel.jl @@ -28,7 +28,7 @@ end ManifoldsBase.@manifold_element_forwards StiefelPoint value ManifoldsBase.@manifold_vector_forwards StiefelTVector value ManifoldsBase.@default_manifold_fallbacks Stiefel StiefelPoint StiefelTVector value value -ManifoldsBase.@default_manifold_fallbacks (Stiefel{n,k,ℝ} where {n,k}) StiefelPoint StiefelTVector value value +ManifoldsBase.@default_manifold_fallbacks (Stiefel{<:Any,ℝ}) StiefelPoint StiefelTVector value value ManifoldsBase.@default_manifold_fallbacks Grassmann StiefelPoint StiefelTVector value value function default_vector_transport_method(::Grassmann, ::Type{<:AbstractArray}) diff --git a/src/manifolds/Rotations.jl b/src/manifolds/Rotations.jl index 4006f8080a..4dee6fb28e 100644 --- a/src/manifolds/Rotations.jl +++ b/src/manifolds/Rotations.jl @@ -408,7 +408,8 @@ 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} +function riemannian_Hessian!(M::Rotations, Y, p, G, H, X) + N = get_n(M) symmetrize!(Y, G' * p) project!(SkewSymmetricMatrices(N), Y, p' * H - X * Y) return Y diff --git a/src/manifolds/Stiefel.jl b/src/manifolds/Stiefel.jl index 2cbee11ac0..81db0a3ce0 100644 --- a/src/manifolds/Stiefel.jl +++ b/src/manifolds/Stiefel.jl @@ -1,5 +1,5 @@ @doc raw""" - Stiefel{n,k,𝔽} <: AbstractDecoratorManifold{𝔽} + Stiefel{T,𝔽} <: AbstractDecoratorManifold{𝔽} The Stiefel manifold consists of all $n × k$, $n ≥ k$ unitary matrices, i.e. @@ -27,19 +27,24 @@ The manifold is named after [Eduard L. Stiefel](https://en.wikipedia.org/wiki/Eduard_Stiefel) (1909–1978). # Constructor - Stiefel(n, k, field = ℝ) + Stiefel(n, k, field = ℝ; parameter::Symbol=:field) Generate the (real-valued) Stiefel manifold of $n × k$ dimensional orthonormal matrices. """ -struct Stiefel{n,k,𝔽} <: AbstractDecoratorManifold{𝔽} end +struct Stiefel{T,𝔽} <: AbstractDecoratorManifold{𝔽} + size::T +end -Stiefel(n::Int, k::Int, field::AbstractNumbers=ℝ) = Stiefel{n,k,field}() +function Stiefel(n::Int, k::Int, field::AbstractNumbers=ℝ; parameter::Symbol=:field) + size = wrap_type_parameter(parameter, (n, k)) + return Stiefel{typeof(size),field}(size) +end function active_traits(f, ::Stiefel, args...) return merge_traits(IsIsometricEmbeddedManifold(), IsDefaultMetric(EuclideanMetric())) end -function allocation_promotion_function(::Stiefel{n,k,ℂ}, ::Any, ::Tuple) where {n,k} +function allocation_promotion_function(::Stiefel{<:Any,ℂ}, ::Any, ::Tuple) return complex end @@ -76,7 +81,7 @@ Check whether `p` is a valid point on the [`Stiefel`](@ref) `M`=$\operatorname{S [`AbstractNumbers`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#number-system) type and $p^{\mathrm{H}}p$ is (approximately) the identity, where $\cdot^{\mathrm{H}}$ is the complex conjugate transpose. The settings for approximately can be set with `kwargs...`. """ -function check_point(M::Stiefel{n,k,𝔽}, p; kwargs...) where {n,k,𝔽} +function check_point(M::Stiefel, p; kwargs...) cks = check_size(M, p) (cks === nothing) || return cks c = p' * p @@ -98,7 +103,8 @@ it (approximately) holds that $p^{\mathrm{H}}X + \overline{X^{\mathrm{H}}p} = 0$ where $\cdot^{\mathrm{H}}$ denotes the Hermitian and $\overline{\cdot}$ the (elementwise) complex conjugate. The settings for approximately can be set with `kwargs...`. """ -function check_vector(M::Stiefel{n,k,𝔽}, p, X; kwargs...) where {n,k,𝔽} +function check_vector(M::Stiefel, p, X; kwargs...) + n, k = get_nk(M) cks = check_size(M, p, X) cks === nothing || return cks if !isapprox(p' * X, -conj(X' * p); kwargs...) @@ -138,9 +144,16 @@ end embed(::Stiefel, p) = p embed(::Stiefel, p, X) = X -function get_embedding(::Stiefel{N,K,𝔽}) where {N,K,𝔽} - return Euclidean(N, K; field=𝔽) +function get_embedding(::Stiefel{TypeParameter{Tuple{n,k}},𝔽}) where {n,k,𝔽} + return Euclidean(n, k; field=𝔽, parameter=:type) end +function get_embedding(M::Stiefel{Tuple{Int,Int},𝔽}) where {𝔽} + n, k = get_nk(M) + return Euclidean(n, k; field=𝔽) +end + +get_nk(::Stiefel{TypeParameter{Tuple{n,k}}}) where {n,k} = (n, k) +get_nk(M::Stiefel{Tuple{Int,Int}}) = get_parameter(M.size) @doc raw""" inverse_retract(M::Stiefel, p, q, ::PolarInverseRetraction) @@ -176,7 +189,8 @@ in [KanekoFioriTanaka:2013](@cite). """ inverse_retract(::Stiefel, ::Any, ::Any, ::QRInverseRetraction) -function _stiefel_inv_retr_qr_mul_by_r_generic!(::Stiefel{n,k}, X, q, R, A) where {n,k} +function _stiefel_inv_retr_qr_mul_by_r_generic!(M::Stiefel, X, q, R, A) + n, k = get_nk(M) @inbounds for i in 1:k b = zeros(eltype(R), i) b[i] = 1 @@ -188,12 +202,18 @@ function _stiefel_inv_retr_qr_mul_by_r_generic!(::Stiefel{n,k}, X, q, R, A) wher return mul!(X, q, R) end -function _stiefel_inv_retr_qr_mul_by_r!(::Stiefel{n,1}, X, q, A, ::Type) where {n} +function _stiefel_inv_retr_qr_mul_by_r!( + ::Stiefel{TypeParameter{Tuple{n,1}}}, + X, + q, + A, + ::Type, +) where {n} @inbounds R = SMatrix{1,1}(inv(A[1, 1])) return mul!(X, q, R) end function _stiefel_inv_retr_qr_mul_by_r!( - M::Stiefel{n,1}, + M::Stiefel{TypeParameter{Tuple{n,1}}}, X, q, A::StaticArray, @@ -201,7 +221,13 @@ function _stiefel_inv_retr_qr_mul_by_r!( ) where {n,ElT} return invoke( _stiefel_inv_retr_qr_mul_by_r!, - Tuple{Stiefel{n,1},typeof(X),typeof(q),AbstractArray,typeof(ElT)}, + Tuple{ + Stiefel{TypeParameter{Tuple{n,1}}}, + typeof(X), + typeof(q), + AbstractArray, + typeof(ElT), + }, M, X, q, @@ -209,7 +235,13 @@ function _stiefel_inv_retr_qr_mul_by_r!( ElT, ) end -function _stiefel_inv_retr_qr_mul_by_r!(::Stiefel{n,2}, X, q, A, ::Type{ElT}) where {n,ElT} +function _stiefel_inv_retr_qr_mul_by_r!( + ::Stiefel{TypeParameter{Tuple{n,2}}}, + X, + q, + A, + ::Type{ElT}, +) where {n,ElT} R11 = inv(A[1, 1]) @inbounds R = hcat(SA[R11, zero(ElT)], A[SOneTo(2), SOneTo(2)] \ SA[-R11 * A[2, 1], one(ElT)]) @@ -219,7 +251,7 @@ function _stiefel_inv_retr_qr_mul_by_r!(::Stiefel{n,2}, X, q, A, ::Type{ElT}) wh return mul!(X, q, R) end function _stiefel_inv_retr_qr_mul_by_r!( - M::Stiefel{n,2}, + M::Stiefel{TypeParameter{Tuple{n,2}}}, X, q, A::StaticArray, @@ -227,7 +259,13 @@ function _stiefel_inv_retr_qr_mul_by_r!( ) where {n,ElT} return invoke( _stiefel_inv_retr_qr_mul_by_r!, - Tuple{Stiefel{n,2},typeof(X),typeof(q),AbstractArray,typeof(ElT)}, + Tuple{ + Stiefel{TypeParameter{Tuple{n,2}}}, + typeof(X), + typeof(q), + AbstractArray, + typeof(ElT), + }, M, X, q, @@ -236,7 +274,7 @@ function _stiefel_inv_retr_qr_mul_by_r!( ) end function _stiefel_inv_retr_qr_mul_by_r!( - M::Stiefel{n,k}, + M::Stiefel{TypeParameter{Tuple{n,k}}}, X, q, A::StaticArray, @@ -245,13 +283,8 @@ function _stiefel_inv_retr_qr_mul_by_r!( R = zeros(MMatrix{k,k,ElT}) return _stiefel_inv_retr_qr_mul_by_r_generic!(M, X, q, R, A) end -function _stiefel_inv_retr_qr_mul_by_r!( - M::Stiefel{n,k}, - X, - q, - A, - ::Type{ElT}, -) where {n,k,ElT} +function _stiefel_inv_retr_qr_mul_by_r!(M::Stiefel, X, q, A, ::Type{ElT}) where {ElT} + n, k = get_nk(M) R = zeros(ElT, k, k) return _stiefel_inv_retr_qr_mul_by_r_generic!(M, X, q, R, A) end @@ -264,7 +297,8 @@ function inverse_retract_polar!(::Stiefel, X, p, q) X .-= p return X end -function inverse_retract_qr!(M::Stiefel{n,k}, X, p, q) where {n,k} +function inverse_retract_qr!(M::Stiefel, X, p, q) + n, k = get_nk(M) A = p' * q @boundscheck size(A) === (k, k) ElT = typeof(one(eltype(p)) * one(eltype(q))) @@ -298,9 +332,18 @@ The dimension is given by \end{aligned} ```` """ -manifold_dimension(::Stiefel{n,k,ℝ}) where {n,k} = n * k - div(k * (k + 1), 2) -manifold_dimension(::Stiefel{n,k,ℂ}) where {n,k} = 2 * n * k - k * k -manifold_dimension(::Stiefel{n,k,ℍ}) where {n,k} = 4 * n * k - k * (2k - 1) +function manifold_dimension(M::Stiefel{<:Any,ℝ}) + n, k = get_nk(M) + return n * k - div(k * (k + 1), 2) +end +function manifold_dimension(M::Stiefel{<:Any,ℂ}) + n, k = get_nk(M) + return 2 * n * k - k * k +end +function manifold_dimension(M::Stiefel{<:Any,ℍ}) + n, k = get_nk(M) + return 4 * n * k - k * (2k - 1) +end @doc raw""" rand(::Stiefel; vector_at=nothing, σ::Real=1.0) @@ -318,11 +361,12 @@ rand(::Stiefel; σ::Real=1.0) function Random.rand!( rng::AbstractRNG, - M::Stiefel{n,k,𝔽}, + M::Stiefel{<:Any,𝔽}, pX; vector_at=nothing, σ::Real=one(real(eltype(pX))), -) where {n,k,𝔽} +) where {𝔽} + n, k = get_nk(M) if vector_at === nothing A = σ * randn(rng, 𝔽 === ℝ ? Float64 : ComplexF64, n, k) pX .= Matrix(qr(A).Q) @@ -471,12 +515,18 @@ end Returns the representation size of the [`Stiefel`](@ref) `M`=$\operatorname{St}(n,k)$, i.e. `(n,k)`, which is the matrix dimensions. """ -@generated representation_size(::Stiefel{n,k}) where {n,k} = (n, k) +representation_size(M::Stiefel) = get_nk(M) -Base.show(io::IO, ::Stiefel{n,k,F}) where {n,k,F} = print(io, "Stiefel($(n), $(k), $(F))") +function Base.show(io::IO, ::Stiefel{TypeParameter{Tuple{n,k}},𝔽}) where {n,k,𝔽} + return print(io, "Stiefel($(n), $(k), $(𝔽); parameter=:type)") +end +function Base.show(io::IO, M::Stiefel{Tuple{Int,Int},𝔽}) where {𝔽} + n, k = get_nk(M) + return print(io, "Stiefel($(n), $(k), $(𝔽))") +end """ - uniform_distribution(M::Stiefel{n,k,ℝ}, p) + uniform_distribution(M::Stiefel{<:Any,ℝ}, p) Uniform distribution on given (real-valued) [`Stiefel`](@ref) `M`. Specifically, this is the normalized Haar and Hausdorff measure on `M`. @@ -485,7 +535,8 @@ Generated points will be of similar type as `p`. 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} +function uniform_distribution(M::Stiefel{<:Any,ℝ}, p) + n, k = get_nk(M) μ = Distributions.Zeros(n, k) σ = one(eltype(p)) Σ1 = Distributions.PDMats.ScalMat(n, σ) diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl index 68d57a6f28..c0e05069db 100644 --- a/src/manifolds/StiefelCanonicalMetric.jl +++ b/src/manifolds/StiefelCanonicalMetric.jl @@ -11,7 +11,7 @@ struct CanonicalMetric <: RiemannianMetric end ApproximateLogarithmicMap <: ApproximateInverseRetraction An approximate implementation of the logarithmic map, which is an [`inverse_retract`](@ref)ion. -See [`inverse_retract(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, ::Any, ::Any, ::ApproximateLogarithmicMap) where {n,k}`](@ref) for a use case. +See [`inverse_retract(::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},CanonicalMetric}, ::Any, ::Any, ::ApproximateLogarithmicMap)`](@ref) for a use case. # Fields @@ -24,15 +24,15 @@ struct ApproximateLogarithmicMap{T} <: ApproximateInverseRetraction tolerance::T end -function distance(M::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, q, p) where {n,k} +function distance(M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},CanonicalMetric}, q, p) return norm(M, p, log(M, p, q)) end @doc raw""" - q = exp(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, X) - exp!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, q, CanonicalMetric}, p, X) + q = exp(M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},CanonicalMetric}, p, X) + exp!(M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},CanonicalMetric}, q, p, X) -Compute the exponential map on the [`Stiefel`](@ref)`(n,k)` manifold with respect to the [`CanonicalMetric`](@ref). +Compute the exponential map on the [`Stiefel`](@ref)`(n, k)` manifold with respect to the [`CanonicalMetric`](@ref). First, decompose The tangent vector ``X`` into its horizontal and vertical component with respect to ``p``, i.e. @@ -64,9 +64,10 @@ q = \exp_p X = pC + QB. ``` For more details, see [EdelmanAriasSmith:1998](@cite)[Zimmermann:2017](@cite). """ -exp(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, ::Any...) where {n,k} +exp(::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},CanonicalMetric}, ::Any...) -function exp!(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, q, p, X) where {n,k} +function exp!(M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},CanonicalMetric}, q, p, X) + n, k = get_nk(M.manifold) A = p' * X n == k && return mul!(q, p, exp(A)) QR = qr(X - p * A) @@ -79,7 +80,7 @@ function exp!(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, q, p, X) w end @doc raw""" - inner(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p, X, Y) + inner(M::MetricManifold{ℝ, Stiefel{<:Any,ℝ}, X, CanonicalMetric}, p, X, Y) Compute the inner product on the [`Stiefel`](@ref) manifold with respect to the [`CanonicalMetric`](@ref). The formula reads @@ -88,7 +89,8 @@ Compute the inner product on the [`Stiefel`](@ref) manifold with respect to the g_p(X,Y) = \operatorname{tr}\bigl( X^{\mathrm{T}}(I_n - \frac{1}{2}pp^{\mathrm{T}})Y \bigr). ``` """ -function inner(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, p, X, Y) where {n,k} +function inner(M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},CanonicalMetric}, p, X, Y) + n, k = get_nk(M.manifold) T = Base.promote_eltype(p, X, Y) if n == k return T(dot(X, Y)) / 2 @@ -98,54 +100,55 @@ function inner(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, p, X, Y) end @doc raw""" - X = inverse_retract(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, q, a::ApproximateLogarithmicMap) - inverse_retract!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p, q, a::ApproximateLogarithmicMap) + X = inverse_retract(M::MetricManifold{ℝ, Stiefel{<:Any,ℝ}, CanonicalMetric}, p, q, a::ApproximateLogarithmicMap) + inverse_retract!(M::MetricManifold{ℝ, Stiefel{<:Any,ℝ}, X, CanonicalMetric}, p, q, a::ApproximateLogarithmicMap) -Compute an approximation to the logarithmic map on the [`Stiefel`](@ref)`(n,k)` manifold with respect to the [`CanonicalMetric`](@ref) +Compute an approximation to the logarithmic map on the [`Stiefel`](@ref)`(n, k)` manifold with respect to the [`CanonicalMetric`](@ref) 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}). +[`exp`](@ref exp(::MetricManifold{ℝ, Stiefel{<:Any,ℝ}, CanonicalMetric}, ::Any...)). The algorithm is derived in [Zimmermann:2017](@cite) and it uses the `max_iterations` and the `tolerance` field from the [`ApproximateLogarithmicMap`](@ref). """ inverse_retract( - ::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, + ::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},CanonicalMetric}, ::Any, ::Any, ::ApproximateLogarithmicMap, -) where {n,k} +) function log( - M::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, + M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},CanonicalMetric}, p, q; maxiter::Int=10000, tolerance=1e-9, -) where {n,k} +) X = allocate_result(M, log, p, q) inverse_retract!(M, X, p, q, ApproximateLogarithmicMap(maxiter, tolerance)) return X end function log!( - M::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, + M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},CanonicalMetric}, X, p, q; maxiter::Int=10000, tolerance=1e-9, -) where {n,k} +) inverse_retract!(M, X, p, q, ApproximateLogarithmicMap(maxiter, tolerance)) return X end function inverse_retract!( - ::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, + M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},CanonicalMetric}, X, p, q, a::ApproximateLogarithmicMap, -) where {n,k} +) + n, k = get_nk(M.manifold) qfact = stiefel_factorization(p, q) V = allocate(qfact.Z, Size(2k, 2k)) LV = allocate(V) @@ -176,8 +179,8 @@ function inverse_retract!( 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) + Y = riemannian_Hessian(M::MetricManifold{ℝ, Stiefel, CanonicalMetric}, p, G, H, X) + riemannian_Hessian!(M::MetricManifold{ℝ, Stiefel, 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`, @@ -196,22 +199,16 @@ Here, we adopt Eq. (5.6) [Nguyen:2023](@cite), for the [`CanonicalMetric`](@ref) ``` where ``P = I-pp^{\mathrm{H}}``. """ -riemannian_Hessian( - M::MetricManifold{𝔽,Stiefel{n,k,𝔽},CanonicalMetric}, - p, - G, - H, - X, -) where {n,k,𝔽} +riemannian_Hessian(M::MetricManifold{𝔽,Stiefel,CanonicalMetric}, p, G, H, X) where {𝔽} function riemannian_Hessian!( - M::MetricManifold{𝔽,Stiefel{n,k,𝔽},CanonicalMetric}, + M::MetricManifold{𝔽,<:Stiefel{<:Any,𝔽},CanonicalMetric}, Y, p, G, H, X, -) where {n,k,𝔽} +) where {𝔽} Gp = symmetrize(G' * p) Z = symmetrize((I - p * p') * G * p') project!(M, Y, p, H - X * Gp - Z * X) diff --git a/src/manifolds/StiefelEuclideanMetric.jl b/src/manifolds/StiefelEuclideanMetric.jl index f897f4053b..50a9c451a4 100644 --- a/src/manifolds/StiefelEuclideanMetric.jl +++ b/src/manifolds/StiefelEuclideanMetric.jl @@ -23,7 +23,8 @@ $0_k$ are the identity matrix and the zero matrix of dimension $k × k$, respect """ exp(::Stiefel, ::Any...) -function exp!(::Stiefel{n,k}, q, p, X) where {n,k} +function exp!(M::Stiefel, q, p, X) + n, k = get_nk(M) A = p' * X B = exp([A -X'*X; I A]) @views begin @@ -35,7 +36,7 @@ function exp!(::Stiefel{n,k}, q, p, X) where {n,k} end @doc raw""" - get_basis(M::Stiefel{n,k,ℝ}, p, B::DefaultOrthonormalBasis) where {n,k} + get_basis(M::Stiefel{<:Any,ℝ}, p, B::DefaultOrthonormalBasis) Create the default basis using the parametrization for any $X ∈ T_p\mathcal M$. Set $p_\bot \in ℝ^{n\times(n-k)}$ the matrix such that the $n\times n$ matrix of the common @@ -58,30 +59,24 @@ trangular entries of $a$ is set to $1$ its symmetric entry to $-1$ and we normal the factor $\frac{1}{\sqrt{2}}$ and for $b$ one can just use unit vectors reshaped to a matrix to obtain orthonormal set of parameters. """ -get_basis(M::Stiefel{n,k,ℝ}, p, B::DefaultOrthonormalBasis{ℝ,TangentSpaceType}) where {n,k} +get_basis(M::Stiefel{<:Any,ℝ}, p, B::DefaultOrthonormalBasis{ℝ,TangentSpaceType}) function _get_basis( - M::Stiefel{n,k,ℝ}, + M::Stiefel{<:Any,ℝ}, p, B::DefaultOrthonormalBasis{ℝ,TangentSpaceType}; kwargs..., -) where {n,k} +) return CachedBasis(B, get_vectors(M, p, B)) end -function get_coordinates_orthonormal!( - M::Stiefel{n,k,ℝ}, - c, - p, - X, - N::RealNumbers, -) where {n,k} +function get_coordinates_orthonormal!(M::Stiefel{<:Any,ℝ}, c, p, X, N::RealNumbers) V = get_vectors(M, p, DefaultOrthonormalBasis(N)) c .= inner.(Ref(M), Ref(p), V, Ref(X)) return c end -function get_vector_orthonormal!(M::Stiefel{n,k,ℝ}, X, p, c, N::RealNumbers) where {n,k} +function get_vector_orthonormal!(M::Stiefel{<:Any,ℝ}, X, p, c, N::RealNumbers) V = get_vectors(M, p, DefaultOrthonormalBasis(N)) zero_vector!(M, X, p) length(c) < length(V) && error( @@ -93,11 +88,8 @@ function get_vector_orthonormal!(M::Stiefel{n,k,ℝ}, X, p, c, N::RealNumbers) w return X end -function get_vectors( - ::Stiefel{n,k,ℝ}, - p, - ::DefaultOrthonormalBasis{ℝ,TangentSpaceType}, -) where {n,k} +function get_vectors(M::Stiefel{<:Any,ℝ}, p, ::DefaultOrthonormalBasis{ℝ,TangentSpaceType}) + n, k = get_nk(M) p⊥ = nullspace([p zeros(n, n - k)]) an = div(k * (k - 1), 2) bn = (n - k) * k @@ -126,7 +118,7 @@ function inverse_retract_project!(M::Stiefel, X, p, q) return X end -function log!(M::Stiefel{n,k,ℝ}, X, p, q) where {n,k} +function log!(M::Stiefel{<:Any,ℝ}, X, p, q) MM = MetricManifold(M, StiefelSubmersionMetric(-1 // 2)) log!(MM, X, p, q) return X diff --git a/src/manifolds/StiefelSubmersionMetric.jl b/src/manifolds/StiefelSubmersionMetric.jl index 4c0fbbe00e..982c7f385f 100644 --- a/src/manifolds/StiefelSubmersionMetric.jl +++ b/src/manifolds/StiefelSubmersionMetric.jl @@ -22,8 +22,8 @@ struct StiefelSubmersionMetric{T<:Real} <: RiemannianMetric end @doc raw""" - q = exp(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, <:StiefelSubmersionMetric}, p, X) - exp!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, q, <:StiefelSubmersionMetric}, p, X) + q = exp(M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, p, X) + exp!(M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, q, p, X) Compute the exponential map on the [`Stiefel(n,k)`](@ref) manifold with respect to the [`StiefelSubmersionMetric`](@ref). @@ -40,14 +40,10 @@ This implementation is based on [ZimmermannHueper:2022](@cite). For ``k < \frac{n}{2}`` the exponential is computed more efficiently using [`StiefelFactorization`](@ref). """ -exp(::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, ::Any...) where {n,k} +exp(::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, ::Any...) -function exp!( - M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, - q, - p, - X, -) where {n,k} +function exp!(M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, q, p, X) + n, k = get_nk(M.manifold) α = metric(M).α T = Base.promote_eltype(q, p, X) if k ≤ div(n, 2) @@ -82,7 +78,7 @@ function exp!( end @doc raw""" - inner(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, <:StiefelSubmersionMetric}, p, X, Y) + inner(M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, p, X, Y) Compute the inner product on the [`Stiefel`](@ref) manifold with respect to the [`StiefelSubmersionMetric`](@ref). The formula reads @@ -91,12 +87,8 @@ g_p(X,Y) = \operatorname{tr}\bigl( X^{\mathrm{T}}(I_n - \frac{2α+1}{2(α+1)}pp^ ``` where ``α`` is the parameter of the metric. """ -function inner( - M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, - p, - X, - Y, -) where {n,k} +function inner(M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, p, X, Y) + n, k = get_nk(M.manifold) α = metric(M).α T = typeof(one(Base.promote_eltype(p, X, Y, α))) if n == k @@ -119,7 +111,7 @@ end @doc doc""" inverse_retract( - M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, + M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, p, q, method::ShootingInverseRetraction, @@ -130,7 +122,7 @@ Compute the inverse retraction using [`ShootingInverseRetraction`](https://julia In general the retraction is computed using the generic shooting method. inverse_retract( - M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, + M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, p, q, method::ShootingInverseRetraction{ @@ -153,7 +145,7 @@ inverse_retract( ) function inverse_retract_shooting!( - M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, + M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, X::AbstractMatrix, p::AbstractMatrix, q::AbstractMatrix, @@ -162,13 +154,14 @@ function inverse_retract_shooting!( ProjectionInverseRetraction, <:Union{ProjectionTransport,ScaledVectorTransport{ProjectionTransport}}, }, -) where {n,k} +) + n, k = get_nk(M.manifold) if k > div(n, 2) # fall back to default method invoke( inverse_retract_shooting!, Tuple{ - MetricManifold{ℝ,Stiefel{n,k,ℝ}}, + MetricManifold{ℝ,typeof(M.manifold)}, typeof(X), typeof(p), typeof(q), @@ -192,7 +185,7 @@ function inverse_retract_shooting!( end @doc raw""" - log(M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, p, q; kwargs...) + log(M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, p, q; kwargs...) Compute the logarithmic map on the [`Stiefel(n,k)`](@ref) manifold with respect to the [`StiefelSubmersionMetric`](@ref). @@ -205,13 +198,13 @@ that documentation for details. Their defaults are: - `max_iterations=1_000` """ function log( - M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, + M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, p, q; tolerance=sqrt(eps(float(real(Base.promote_eltype(p, q))))), max_iterations=1_000, num_transport_points=4, -) where {n,k} +) X = allocate_result(M, log, p, q) log!( M, @@ -225,14 +218,14 @@ function log( return X end function log!( - M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, + M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, X, p, q; tolerance=sqrt(eps(float(real(eltype(X))))), max_iterations=1_000, num_transport_points=4, -) where {n,k} +) retraction = ExponentialRetraction() initial_inverse_retraction = ProjectionInverseRetraction() vector_transport = ScaledVectorTransport(ProjectionTransport()) @@ -248,8 +241,8 @@ function log!( 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) + Y = riemannian_Hessian(M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},StiefelSubmersionMetric}, p, G, H, X) + riemannian_Hessian!(MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},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`, @@ -271,21 +264,21 @@ 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}, + M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, p, G, H, X, -) where {n,k} +) function riemannian_Hessian!( - M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, + M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, Y, p, G, H, X, -) where {n,k} +) α = metric(M).α Gp = symmetrize(G' * p) Z = symmetrize((I - p * p') * G * p') @@ -416,38 +409,38 @@ function Base.copyto!( broadcast!(bc.f, dest.Z, Zargs...) return dest end -function project!( - ::Stiefel{n,k,ℝ}, - q::StiefelFactorization, - p::StiefelFactorization, -) where {n,k} +function project!(M::Stiefel{<:Any,ℝ}, q::StiefelFactorization, p::StiefelFactorization) + n, k = get_nk(M) project!(Stiefel(2k, k), q.Z, p.Z) return q end function project!( - ::Stiefel{n,k,ℝ}, + M::Stiefel{<:Any,ℝ}, Y::StiefelFactorization, p::StiefelFactorization, X::StiefelFactorization, -) where {n,k} +) + n, k = get_nk(M) project!(Stiefel(2k, k), Y.Z, p.Z, X.Z) return Y end function inner( - M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, + M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, p::StiefelFactorization, X::StiefelFactorization, Y::StiefelFactorization, -) where {n,k} +) + n, k = get_nk(M.manifold) Msub = MetricManifold(Stiefel(2k, k), metric(M)) return inner(Msub, p.Z, X.Z, Y.Z) end function exp!( - M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, + M::MetricManifold{ℝ,<:Stiefel{<:Any,ℝ},<:StiefelSubmersionMetric}, q::StiefelFactorization, p::StiefelFactorization, X::StiefelFactorization, -) where {n,k} +) + n, k = get_nk(M.manifold) α = metric(M).α @views begin ZM = X.Z[1:k, 1:k] diff --git a/src/manifolds/SymmetricPositiveDefinite.jl b/src/manifolds/SymmetricPositiveDefinite.jl index 44e974741d..d70b1d9282 100644 --- a/src/manifolds/SymmetricPositiveDefinite.jl +++ b/src/manifolds/SymmetricPositiveDefinite.jl @@ -1,5 +1,5 @@ @doc raw""" - SymmetricPositiveDefinite{N} <: AbstractDecoratorManifold{ℝ} + SymmetricPositiveDefinite{T} <: AbstractDecoratorManifold{ℝ} The manifold of symmetric positive definite matrices, i.e. @@ -22,13 +22,18 @@ i.e. the set of symmetric matrices, # Constructor - SymmetricPositiveDefinite(n) + SymmetricPositiveDefinite(n; parameter::Symbol=:field) generates the manifold $\mathcal P(n) \subset ℝ^{n × n}$ """ -struct SymmetricPositiveDefinite{N} <: AbstractDecoratorManifold{ℝ} end +struct SymmetricPositiveDefinite{T} <: AbstractDecoratorManifold{ℝ} + size::T +end -SymmetricPositiveDefinite(n::Int) = SymmetricPositiveDefinite{n}() +function SymmetricPositiveDefinite(n::Int; parameter::Symbol=:field) + size = wrap_type_parameter(parameter, (n,)) + return SymmetricPositiveDefinite{typeof(size)}(size) +end @doc raw""" SPDPoint <: AbstractManifoldsPoint @@ -131,7 +136,7 @@ checks, whether `p` is a valid point on the [`SymmetricPositiveDefinite`](@ref) of size `(N,N)`, symmetric and positive definite. The tolerance for the second to last test can be set using the `kwargs...`. """ -function check_point(M::SymmetricPositiveDefinite{N}, p; kwargs...) where {N} +function check_point(M::SymmetricPositiveDefinite, p; kwargs...) if !isapprox(norm(p - transpose(p)), 0.0; kwargs...) return DomainError( norm(p - transpose(p)), @@ -159,7 +164,7 @@ and a symmetric matrix, i.e. this stores tangent vetors as elements of the corre Lie group. The tolerance for the last test can be set using the `kwargs...`. """ -function check_vector(M::SymmetricPositiveDefinite{N}, p, X; kwargs...) where {N} +function check_vector(M::SymmetricPositiveDefinite, p, X; kwargs...) if !isapprox(norm(X - transpose(X)), 0.0; kwargs...) return DomainError( X, @@ -227,6 +232,9 @@ function get_embedding(M::SymmetricPositiveDefinite) return Euclidean(representation_size(M)...; field=ℝ) end +get_n(::SymmetricPositiveDefinite{TypeParameter{N}}) where {N} = N +get_n(M::SymmetricPositiveDefinite{Tuple{Int}}) = get_parameter(M.size)[1] + @doc raw""" injectivity_radius(M::SymmetricPositiveDefinite[, p]) injectivity_radius(M::MetricManifold{SymmetricPositiveDefinite,AffineInvariantMetric}[, p]) @@ -261,8 +269,9 @@ returns the dimension of \dim \mathcal P(n) = \frac{n(n+1)}{2}. ```` """ -@generated function manifold_dimension(::SymmetricPositiveDefinite{N}) where {N} - return div(N * (N + 1), 2) +function manifold_dimension(M::SymmetricPositiveDefinite) + n = get_n(M) + return div(n * (n + 1), 2) end """ @@ -405,13 +414,14 @@ rand(M::SymmetricPositiveDefinite; σ::Real=1) function Random.rand!( rng::AbstractRNG, - M::SymmetricPositiveDefinite{N}, + M::SymmetricPositiveDefinite, pX; vector_at=nothing, σ::Real=one(eltype(pX)) / (vector_at === nothing ? 1 : norm(convert(AbstractMatrix, vector_at))), tangent_distr=:Gaussian, -) where {N} +) + N = get_n(M) if vector_at === nothing D = Diagonal(1 .+ rand(rng, N)) # random diagonal matrix s = qr(σ * randn(rng, N, N)) # random q @@ -454,10 +464,17 @@ Return the size of an array representing an element on the [`SymmetricPositiveDefinite`](@ref) manifold `M`, i.e. $n × n$, the size of such a symmetric positive definite matrix on $\mathcal M = \mathcal P(n)$. """ -@generated representation_size(::SymmetricPositiveDefinite{N}) where {N} = (N, N) +function representation_size(M::SymmetricPositiveDefinite) + N = get_n(M) + return (N, N) +end -function Base.show(io::IO, ::SymmetricPositiveDefinite{N}) where {N} - return print(io, "SymmetricPositiveDefinite($(N))") +function Base.show(io::IO, ::SymmetricPositiveDefinite{TypeParameter{Tuple{n}}}) where {n} + return print(io, "SymmetricPositiveDefinite($(n); parameter=:type)") +end +function Base.show(io::IO, M::SymmetricPositiveDefinite{Tuple{Int}}) + n = get_n(M) + return print(io, "SymmetricPositiveDefinite($(n))") end function Base.show(io::IO, ::MIME"text/plain", p::SPDPoint) @@ -489,4 +506,4 @@ function zero_vector(M::SymmetricPositiveDefinite, p::SPDPoint) return zero_vector(M, convert(AbstractMatrix, p)) end -zero_vector!(::SymmetricPositiveDefinite{N}, X, ::Any) where {N} = fill!(X, 0) +zero_vector!(::SymmetricPositiveDefinite, X, ::Any) = fill!(X, 0) diff --git a/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl b/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl index dde2356aef..c41acc107c 100644 --- a/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl +++ b/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl @@ -31,7 +31,7 @@ function change_representer!(::SymmetricPositiveDefinite, Y, ::EuclideanMetric, end @doc raw""" - change_metric(M::SymmetricPositiveDefinite{n}, E::EuclideanMetric, p, X) + change_metric(M::SymmetricPositiveDefinite, E::EuclideanMetric, p, X) Given a tangent vector ``X ∈ T_p\mathcal P(n)`` with respect to the [`EuclideanMetric`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.EuclideanMetric) `g_E`, this function changes into the [`AffineInvariantMetric`](@ref) (default) metric on the @@ -67,7 +67,7 @@ d_{\mathcal P(n)}(p,q) where $\operatorname{Log}$ denotes the matrix logarithm and $\lVert\cdot\rVert_{\mathrm{F}}$ denotes the matrix Frobenius norm. """ -function distance(::SymmetricPositiveDefinite{N}, p, q) where {N} +function distance(::SymmetricPositiveDefinite, p, q) # avoid numerical instabilities in cholesky norm(p - q) < eps(eltype(p + q)) && return zero(eltype(p + q)) cq = cholesky(Symmetric(q)) # to avoid numerical inaccuracies @@ -80,7 +80,7 @@ end @doc raw""" exp(M::SymmetricPositiveDefinite, p, X) - exp(M::MetricManifold{SymmetricPositiveDefinite{N},AffineInvariantMetric}, p, X) + exp(M::MetricManifold{<:SymmetricPositiveDefinite,AffineInvariantMetric}, p, X) Compute the exponential map from `p` with tangent vector `X` on the [`SymmetricPositiveDefinite`](@ref) `M` with its default [`MetricManifold`](@ref) having the @@ -95,7 +95,7 @@ where $\operatorname{Exp}$ denotes to the matrix exponential. exp(::SymmetricPositiveDefinite, ::Any...) exp(M::SymmetricPositiveDefinite, p::SPDPoint, X, t::Number) = exp(M, p, t * X) -function exp(::SymmetricPositiveDefinite{N}, p::SPDPoint, X) where {N} +function exp(::SymmetricPositiveDefinite, p::SPDPoint, X) (p_sqrt, p_sqrt_inv) = spd_sqrt_and_sqrt_inv(p) T = Symmetric(p_sqrt_inv * X * p_sqrt_inv) eig1 = eigen(T) # numerical stabilization @@ -114,7 +114,7 @@ end function exp!(M::SymmetricPositiveDefinite, q, p, X, t::Number) return exp!(M, q, p, t * X) end -function exp!(::SymmetricPositiveDefinite{N}, q, p, X) where {N} +function exp!(::SymmetricPositiveDefinite, q, p, X) (p_sqrt, p_sqrt_inv) = spd_sqrt_and_sqrt_inv(p) T = Symmetric(p_sqrt_inv * X * p_sqrt_inv) eig1 = eigen(T) # numerical stabilization @@ -123,7 +123,7 @@ function exp!(::SymmetricPositiveDefinite{N}, q, p, X) where {N} pUe = p_sqrt * Ue return copyto!(q, pUe * Se * transpose(pUe)) end -function exp!(::SymmetricPositiveDefinite{N}, q::SPDPoint, p, X) where {N} +function exp!(::SymmetricPositiveDefinite, q::SPDPoint, p, X) (p_sqrt, p_sqrt_inv) = spd_sqrt_and_sqrt_inv(p) T = Symmetric(p_sqrt_inv * X * p_sqrt_inv) eig1 = eigen(T) # numerical stabilization @@ -146,7 +146,7 @@ end @doc raw""" [Ξ,κ] = get_basis_diagonalizing(M::SymmetricPositiveDefinite, p, B::DiagonalizingOrthonormalBasis) - [Ξ,κ] = get_basis_diagonalizing(M::MetricManifold{SymmetricPositiveDefinite{N},AffineInvariantMetric}, p, B::DiagonalizingOrthonormalBasis) + [Ξ,κ] = get_basis_diagonalizing(M::MetricManifold{<:SymmetricPositiveDefinite,AffineInvariantMetric}, p, B::DiagonalizingOrthonormalBasis) Return a orthonormal basis `Ξ` as a vector of tangent vectors (of length [`manifold_dimension`](@ref) of `M`) in the tangent space of `p` on the @@ -158,10 +158,11 @@ The construction is based on an ONB for the symmetric matrices similar to [`get_ just that the ONB here is build from the eigen vectors of ``p^{\frac{1}{2}}Vp^{\frac{1}{2}}``. """ function get_basis_diagonalizing( - ::SymmetricPositiveDefinite{N}, + M::SymmetricPositiveDefinite, p, B::DiagonalizingOrthonormalBasis, -) where {N} +) + N = get_n(M) (p_sqrt, p_sqrt_inv) = spd_sqrt_and_sqrt_inv(p) eigv = eigen(Symmetric(p_sqrt_inv * B.frame_direction * p_sqrt_inv)) V = eigv.vectors @@ -178,7 +179,7 @@ end @doc raw""" [Ξ,κ] = get_basis(M::SymmetricPositiveDefinite, p, B::DefaultOrthonormalBasis) - [Ξ,κ] = get_basis(M::MetricManifold{SymmetricPositiveDefinite{N},AffineInvariantMetric}, p, B::DefaultOrthonormalBasis) + [Ξ,κ] = get_basis(M::MetricManifold{<:SymmetricPositiveDefinite,AffineInvariantMetric}, p, B::DefaultOrthonormalBasis) Return a default ONB for the tangent space ``T_p\mathcal P(n)`` of the [`SymmetricPositiveDefinite`](@ref) with respect to the [`AffineInvariantMetric`](@ref). @@ -208,11 +209,8 @@ We then form the ONB by """ get_basis(::SymmetricPositiveDefinite, p, B::DefaultOrthonormalBasis) -function get_basis_orthonormal( - M::SymmetricPositiveDefinite{N}, - p, - Ns::RealNumbers, -) where {N} +function get_basis_orthonormal(M::SymmetricPositiveDefinite, p, Ns::RealNumbers) + N = get_n(M) p_sqrt = spd_sqrt(p) Ξ = [similar(convert(AbstractMatrix, p)) for _ in 1:manifold_dimension(M)] k = 1 @@ -240,13 +238,8 @@ where $k$ is trhe linearized index of the $i=1,\ldots,n, j=i,\ldots,n$. """ get_coordinates(::SymmetricPositiveDefinite, c, p, X, ::DefaultOrthonormalBasis) -function get_coordinates_orthonormal!( - M::SymmetricPositiveDefinite{N}, - c, - p, - X, - ::RealNumbers, -) where {N} +function get_coordinates_orthonormal!(M::SymmetricPositiveDefinite, c, p, X, ::RealNumbers) + N = get_n(M) dim = manifold_dimension(M) @assert size(c) == (dim,) @assert size(X) == (N, N) @@ -283,13 +276,8 @@ where $k$ is the linearized index of the $i=1,\ldots,n, j=i,\ldots,n$. """ get_vector(::SymmetricPositiveDefinite, X, p, c, ::DefaultOrthonormalBasis) -function get_vector_orthonormal!( - ::SymmetricPositiveDefinite{N}, - X, - p, - c, - ::RealNumbers, -) where {N} +function get_vector_orthonormal!(M::SymmetricPositiveDefinite, X, p, c, ::RealNumbers) + N = get_n(M) @assert size(c) == (div(N * (N + 1), 2),) @assert size(X) == (N, N) p_sqrt = spd_sqrt(p) @@ -359,7 +347,7 @@ function allocate_result( return allocate_result(M, log, convert(AbstractMatrix, q), convert(AbstractMatrix, p)) end -function log!(::SymmetricPositiveDefinite{N}, X, p, q) where {N} +function log!(::SymmetricPositiveDefinite, X, p, q) (p_sqrt, p_sqrt_inv) = spd_sqrt_and_sqrt_inv(p) T = Symmetric(p_sqrt_inv * convert(AbstractMatrix, q) * p_sqrt_inv) e2 = eigen(T) @@ -402,7 +390,7 @@ and `log` the logarithmic map on [`SymmetricPositiveDefinite`](@ref) """ parallel_transport_to(::SymmetricPositiveDefinite, ::Any, ::Any, ::Any) -function parallel_transport_to!(M::SymmetricPositiveDefinite{N}, Y, p, X, q) where {N} +function parallel_transport_to!(M::SymmetricPositiveDefinite, Y, p, X, q) distance(M, p, q) < 2 * eps(eltype(p)) && copyto!(Y, X) (p_sqrt, p_sqrt_inv) = spd_sqrt_and_sqrt_inv(p) tv = Symmetric(p_sqrt_inv * X * p_sqrt_inv) # p^(-1/2)Xp^{-1/2} @@ -467,13 +455,13 @@ function riemann_tensor!(::SymmetricPositiveDefinite, Xresult, p, X, Y, Z) end """ - volume_density(::SymmetricPositiveDefinite{n}, p, X) where {n} + volume_density(::SymmetricPositiveDefinite, p, X) Compute the volume density of the [`SymmetricPositiveDefinite`](@ref) manifold at `p` 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. """ -function volume_density(::SymmetricPositiveDefinite{n}, p, X) where {n} +function volume_density(::SymmetricPositiveDefinite, p, X) eig = eigvals(X) dens = 1.0 for i in 1:length(eig) diff --git a/src/manifolds/SymmetricPositiveDefiniteBuresWasserstein.jl b/src/manifolds/SymmetricPositiveDefiniteBuresWasserstein.jl index 95e06819bb..84732f071f 100644 --- a/src/manifolds/SymmetricPositiveDefiniteBuresWasserstein.jl +++ b/src/manifolds/SymmetricPositiveDefiniteBuresWasserstein.jl @@ -6,7 +6,7 @@ The Bures Wasserstein metric for symmetric positive definite matrices [MalagoMon struct BuresWassersteinMetric <: RiemannianMetric end @doc raw""" - change_representer(M::MetricManifold{ℝ,SymmetricPositiveDefinite,BuresWassersteinMetric}, E::EuclideanMetric, p, X) + change_representer(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,BuresWassersteinMetric}, E::EuclideanMetric, p, X) Given a tangent vector ``X ∈ T_p\mathcal M`` representing a linear function on the tangent space at `p` with respect to the [`EuclideanMetric`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.EuclideanMetric) @@ -23,7 +23,7 @@ for all ``Y`` and hence we get ``Z``= 2(A+A^{\mathrm{T}})`` with ``A=Xp``. """ change_representer( - ::MetricManifold{ℝ,SymmetricPositiveDefinite,BuresWassersteinMetric}, + ::MetricManifold{ℝ,<:SymmetricPositiveDefinite,BuresWassersteinMetric}, ::EuclideanMetric, p, X, @@ -73,7 +73,7 @@ the [`BuresWassersteinMetric`](@ref) given by where ``q=L_p(X)`` denotes the Lyapunov operator, i.e. it solves ``pq + qp = X``. """ -exp(::MetricManifold{ℝ,SymmetricPositiveDefinite,BuresWassersteinMetric}, p, X) +exp(::MetricManifold{ℝ,<:SymmetricPositiveDefinite,BuresWassersteinMetric}, p, X) function exp!( ::MetricManifold{ℝ,<:SymmetricPositiveDefinite,BuresWassersteinMetric}, @@ -127,7 +127,7 @@ the [`BuresWassersteinMetric`](@ref) given by where ``q=L_p(X)`` denotes the Lyapunov operator, i.e. it solves ``pq + qp = X``. """ -log(::MetricManifold{ℝ,SymmetricPositiveDefinite,BuresWassersteinMetric}, p, q) +log(::MetricManifold{ℝ,<:SymmetricPositiveDefinite,BuresWassersteinMetric}, p, q) function log!( ::MetricManifold{ℝ,<:SymmetricPositiveDefinite,BuresWassersteinMetric}, diff --git a/src/manifolds/SymmetricPositiveDefiniteGeneralizedBuresWasserstein.jl b/src/manifolds/SymmetricPositiveDefiniteGeneralizedBuresWasserstein.jl index b710406fa3..262654a864 100644 --- a/src/manifolds/SymmetricPositiveDefiniteGeneralizedBuresWasserstein.jl +++ b/src/manifolds/SymmetricPositiveDefiniteGeneralizedBuresWasserstein.jl @@ -12,7 +12,7 @@ struct GeneralizedBuresWassersteinMetric{T<:AbstractMatrix} <: RiemannianMetric end @doc raw""" - change_representer(M::MetricManifold{ℝ,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, E::EuclideanMetric, p, X) + change_representer(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric}, E::EuclideanMetric, p, X) Given a tangent vector ``X ∈ T_p\mathcal M`` representing a linear function on the tangent space at `p` with respect to the [`EuclideanMetric`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.EuclideanMetric) @@ -28,7 +28,7 @@ it holds for all ``Y`` and hence we get ``Z = 2pXM + 2MXp``. """ change_representer( - ::MetricManifold{ℝ,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, + ::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric}, ::EuclideanMetric, p, X, @@ -67,7 +67,7 @@ function distance( end @doc raw""" - exp(::MatricManifold{ℝ,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, X) + exp(::MatricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric}, p, X) Compute the exponential map on [`SymmetricPositiveDefinite`](@ref) with respect to the [`GeneralizedBuresWassersteinMetric`](@ref) given by @@ -78,7 +78,11 @@ the [`GeneralizedBuresWassersteinMetric`](@ref) given by where ``q=L_{M,p}(X)`` denotes the generalized Lyapunov operator, i.e. it solves ``pqM + Mqp = X``. """ -exp(::MetricManifold{ℝ,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, X) +exp( + ::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric}, + p, + X, +) function exp!( M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric}, @@ -95,7 +99,7 @@ function exp!( end @doc raw""" - inner(::MetricManifold{ℝ,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, X, Y) + inner(::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric}, p, X, Y) Compute the inner product [`SymmetricPositiveDefinite`](@ref) with respect to the [`GeneralizedBuresWassersteinMetric`](@ref) given by @@ -122,13 +126,13 @@ Return false. [`SymmetricPositiveDefinite`](@ref) with [`GeneralizedBuresWassers is not a flat manifold. """ function is_flat( - M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric}, + ::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric}, ) return false end @doc raw""" - log(::MatricManifold{SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, q) + log(::MatricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric}, p, q) Compute the logarithmic map on [`SymmetricPositiveDefinite`](@ref) with respect to the [`BuresWassersteinMetric`](@ref) given by @@ -137,7 +141,11 @@ the [`BuresWassersteinMetric`](@ref) given by \log_p(q) = M(M^{-1}pM^{-1}q)^{\frac{1}{2}} + (qM^{-1}pM^{-1})^{\frac{1}{2}}M - 2 p. ``` """ -log(::MetricManifold{ℝ,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, q) +log( + ::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric}, + p, + q, +) function log!( M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric}, diff --git a/src/manifolds/SymmetricPositiveDefiniteLogCholesky.jl b/src/manifolds/SymmetricPositiveDefiniteLogCholesky.jl index 703dc0b241..6ac02ede73 100644 --- a/src/manifolds/SymmetricPositiveDefiniteLogCholesky.jl +++ b/src/manifolds/SymmetricPositiveDefiniteLogCholesky.jl @@ -35,12 +35,9 @@ where $x$ and $y$ are the cholesky factors of $p$ and $q$, respectively, $⌊\cdot⌋$ denbotes the strictly lower triangular matrix of its argument, and $\lVert\cdot\rVert_{\mathrm{F}}$ the Frobenius norm. """ -function distance( - ::MetricManifold{ℝ,SymmetricPositiveDefinite{N},LogCholeskyMetric}, - p, - q, -) where {N} - return distance(CholeskySpace{N}(), cholesky(p).L, cholesky(q).L) +function distance(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, p, q) + N = get_n(M.manifold) + return distance(CholeskySpace(N), cholesky(p).L, cholesky(q).L) end @doc raw""" @@ -60,28 +57,24 @@ denotes the lower triangular matrix with the diagonal multiplied by $\frac{1}{2} """ exp(::MetricManifold{ℝ,SymmetricPositiveDefinite,LogCholeskyMetric}, ::Any...) -function exp!( - ::MetricManifold{ℝ,SymmetricPositiveDefinite{N},LogCholeskyMetric}, - q, - p, - X, -) where {N} +function exp!(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, q, p, X) + N = get_n(M.manifold) (y, W) = spd_to_cholesky(p, X) - z = exp(CholeskySpace{N}(), y, W) + z = exp(CholeskySpace(N), y, W) return copyto!(q, z * z') end function exp!( - M::MetricManifold{ℝ,SymmetricPositiveDefinite{N},LogCholeskyMetric}, + M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, q, p, X, t::Number, -) where {N} +) return exp!(M, q, p, t * X) end @doc raw""" - inner(M::MetricManifold{LogCholeskyMetric,ℝ,SymmetricPositiveDefinite}, p, X, Y) + inner(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, p, X, Y) Compute the inner product of two matrices `X`, `Y` in the tangent space of `p` on the [`SymmetricPositiveDefinite`](@ref) manifold `M`, as @@ -96,15 +89,11 @@ $z$ is the cholesky factor of $p$, $a_z(W) = z (z^{-1}Wz^{-\mathrm{T}})_{\frac{1}{2}}$, and $(\cdot)_\frac{1}{2}$ denotes the lower triangular matrix with the diagonal multiplied by $\frac{1}{2}$ """ -function inner( - ::MetricManifold{ℝ,SymmetricPositiveDefinite{N},LogCholeskyMetric}, - p, - X, - Y, -) where {N} +function inner(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, p, X, Y) + N = get_n(M.manifold) (z, Xz) = spd_to_cholesky(p, X) (z, Yz) = spd_to_cholesky(p, z, Y) - return inner(CholeskySpace{N}(), z, Xz, Yz) + return inner(CholeskySpace(N), z, Xz, Yz) end """ @@ -116,7 +105,7 @@ is not a flat manifold. is_flat(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}) = false @doc raw""" - log(M::MetricManifold{SymmetricPositiveDefinite,LogCholeskyMetric}, p, q) + log(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, p, q) Compute the logarithmic map on [`SymmetricPositiveDefinite`](@ref) `M` with respect to the [`LogCholeskyMetric`](@ref) emanating from `p` to `q`. @@ -129,21 +118,17 @@ of $q$ and the just mentioned logarithmic map is the one on [`CholeskySpace`](@r """ log(::MetricManifold{ℝ,SymmetricPositiveDefinite,LogCholeskyMetric}, ::Any...) -function log!( - ::MetricManifold{ℝ,SymmetricPositiveDefinite{N},LogCholeskyMetric}, - X, - p, - q, -) where {N} +function log!(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, X, p, q) + N = get_n(M.manifold) x = cholesky(p).L y = cholesky(q).L - log!(CholeskySpace{N}(), X, x, y) + log!(CholeskySpace(N), X, x, y) return tangent_cholesky_to_tangent_spd!(x, X) end @doc raw""" vector_transport_to( - M::MetricManifold{SymmetricPositiveDefinite,LogCholeskyMetric}, + M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, p, X, q, @@ -163,21 +148,22 @@ transport on [`CholeskySpace`](@ref) from $x$ to $y$. The formula hear reads ```` """ parallel_transport_to( - ::MetricManifold{ℝ,SymmetricPositiveDefinite,LogCholeskyMetric}, + ::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, ::Any, ::Any, ::Any, ) function parallel_transport_to!( - ::MetricManifold{ℝ,SymmetricPositiveDefinite{N},LogCholeskyMetric}, + M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, Y, p, X, q, -) where {N} +) + N = get_n(M.manifold) y = cholesky(q).L (x, W) = spd_to_cholesky(p, X) - parallel_transport_to!(CholeskySpace{N}(), Y, x, W, y) + parallel_transport_to!(CholeskySpace(N), Y, x, W, y) return tangent_cholesky_to_tangent_spd!(y, Y) end diff --git a/src/manifolds/SymmetricPositiveDefiniteLogEuclidean.jl b/src/manifolds/SymmetricPositiveDefiniteLogEuclidean.jl index 579facd4cd..ad579dffc1 100644 --- a/src/manifolds/SymmetricPositiveDefiniteLogEuclidean.jl +++ b/src/manifolds/SymmetricPositiveDefiniteLogEuclidean.jl @@ -7,7 +7,7 @@ into the Lie Algebra, i.e. performing a matrix logarithm beforehand. struct LogEuclideanMetric <: RiemannianMetric end @doc raw""" - distance(M::MetricManifold{SymmetricPositiveDefinite{N},LogEuclideanMetric}, p, q) + distance(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogEuclideanMetric}, p, q) Compute the distance on the [`SymmetricPositiveDefinite`](@ref) manifold between `p` and `q` as a [`MetricManifold`](@ref) with [`LogEuclideanMetric`](@ref). @@ -20,11 +20,7 @@ The formula reads where $\operatorname{Log}$ denotes the matrix logarithm and $\lVert\cdot\rVert_{\mathrm{F}}$ denotes the matrix Frobenius norm. """ -function distance( - ::MetricManifold{ℝ,SymmetricPositiveDefinite{N},LogEuclideanMetric}, - p, - q, -) where {N} +function distance(::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogEuclideanMetric}, p, q) return norm(log(Symmetric(p)) - log(Symmetric(q))) end diff --git a/test/groups/special_orthogonal.jl b/test/groups/special_orthogonal.jl index 9f3f0faa06..e8a862d4b3 100644 --- a/test/groups/special_orthogonal.jl +++ b/test/groups/special_orthogonal.jl @@ -3,10 +3,10 @@ include("group_utils.jl") @testset "Special Orthogonal group" begin for n in [2, 3] - G = SpecialOrthogonal(n) - @test repr(G) == "SpecialOrthogonal($n)" + G = SpecialOrthogonal(n; parameter=:type) + @test repr(G) == "SpecialOrthogonal($n; parameter=:type)" M = base_manifold(G) - @test M === Rotations(n) + @test M === Rotations(n; parameter=:type) p = Matrix(I, n, n) @test is_default_metric(MetricManifold(G, EuclideanMetric())) @test is_flat(G) == (n == 2) diff --git a/test/groups/translation_group.jl b/test/groups/translation_group.jl index f3ea6b0285..58c3520705 100644 --- a/test/groups/translation_group.jl +++ b/test/groups/translation_group.jl @@ -6,8 +6,8 @@ include("group_utils.jl") G = TranslationGroup(2, 3) @test repr(G) == "TranslationGroup(2, 3; field = ℝ)" @test repr(TranslationGroup(2, 3; field=ℂ)) == "TranslationGroup(2, 3; field = ℂ)" - @test repr(TranslationGroup(2, 3; field=ℂ, parameter=:field)) == - "TranslationGroup(2, 3; field = ℂ, parameter = :field)" + @test repr(TranslationGroup(2, 3; field=ℂ, parameter=:type)) == + "TranslationGroup(2, 3; field = ℂ, parameter = :type)" @test has_invariant_metric(G, LeftForwardAction()) @test has_invariant_metric(G, RightBackwardAction()) diff --git a/test/manifolds/stiefel.jl b/test/manifolds/stiefel.jl index 2daa966b0b..3f458b328c 100644 --- a/test/manifolds/stiefel.jl +++ b/test/manifolds/stiefel.jl @@ -2,10 +2,10 @@ include("../utils.jl") @testset "Stiefel" begin @testset "Real" begin - M = Stiefel(3, 2) + M = Stiefel(3, 2; parameter=:type) M2 = MetricManifold(M, EuclideanMetric()) @testset "Basics" begin - @test repr(M) == "Stiefel(3, 2, ℝ)" + @test repr(M) == "Stiefel(3, 2, ℝ; parameter=:type)" p = [1.0 0.0; 0.0 1.0; 0.0 0.0] @test is_default_metric(M, EuclideanMetric()) @test representation_size(M) == (3, 2) @@ -206,9 +206,9 @@ include("../utils.jl") end @testset "Complex" begin - M = Stiefel(3, 2, ℂ) + M = Stiefel(3, 2, ℂ; parameter=:type) @testset "Basics" begin - @test repr(M) == "Stiefel(3, 2, ℂ)" + @test repr(M) == "Stiefel(3, 2, ℂ; parameter=:type)" @test representation_size(M) == (3, 2) @test manifold_dimension(M) == 8 @test !is_flat(M)