From 9fc772a00e256db7f0ad51dfe3f063849128942c Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Wed, 22 Sep 2021 19:46:01 +0200 Subject: [PATCH] Sketch change metric. (#423) * Adds a change_metric and a change_representer * Apply suggestions from code review * bump version. Co-authored-by: Mateusz Baran --- Project.toml | 2 +- src/Manifolds.jl | 10 +- src/groups/group.jl | 2 +- src/manifolds/CenteredMatrices.jl | 4 +- src/manifolds/GeneralizedGrassmann.jl | 33 +++ src/manifolds/Hyperbolic.jl | 5 +- src/manifolds/HyperbolicHyperboloid.jl | 28 ++- src/manifolds/HyperbolicPoincareBall.jl | 45 +++++ src/manifolds/MetricManifold.jl | 190 ++++++++++++++++-- src/manifolds/Multinomial.jl | 1 + src/manifolds/PositiveNumbers.jl | 49 +++++ src/manifolds/PowerManifold.jl | 42 ++++ src/manifolds/ProbabilitySimplex.jl | 38 ++++ src/manifolds/ProductManifold.jl | 36 ++++ src/manifolds/Spectrahedron.jl | 2 +- src/manifolds/Sphere.jl | 6 +- .../SymmetricPositiveDefiniteLinearAffine.jl | 47 +++++ src/riemannian_diff.jl | 5 +- test/manifolds/generalized_grassmann.jl | 75 +++---- test/manifolds/hyperbolic.jl | 29 ++- test/manifolds/positive_numbers.jl | 6 + test/manifolds/power_manifold.jl | 19 ++ test/manifolds/probability_simplex.jl | 6 + test/manifolds/product_manifold.jl | 22 ++ test/manifolds/sphere.jl | 9 + test/manifolds/symmetric_positive_definite.jl | 33 +-- test/metric.jl | 69 ++++++- 27 files changed, 736 insertions(+), 77 deletions(-) diff --git a/Project.toml b/Project.toml index b1cca42bea..71105fc8df 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.6.7" +version = "0.6.8" [deps] Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" diff --git a/src/Manifolds.jl b/src/Manifolds.jl index ddf474c2ef..fb29fdee60 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -108,7 +108,6 @@ using ManifoldsBase: AbstractLinearVectorTransportMethod, ApproximateInverseRetraction, ApproximateRetraction, - DifferentiatedRetractionVectorTransport, ComponentManifoldError, CompositeManifoldError, CotangentSpaceType, @@ -460,6 +459,10 @@ export ×, allocate_result, base_manifold, bundle_projection, + change_metric, + change_metric!, + change_representer, + change_representer!, check_point, check_vector, christoffel_symbols_first, @@ -648,7 +651,10 @@ export get_basis, get_coordinates, get_coordinates!, get_vector, get_vector!, get_vectors, number_system # differentiation export AbstractDiffBackend, - AbstractRiemannianDiffBackend, FiniteDifferencesBackend, RiemannianONBDiffBackend + AbstractRiemannianDiffBackend, + FiniteDifferencesBackend, + RiemannianONBDiffBackend, + RiemannianProjectionGradientBackend export diff_backend, diff_backend!, diff_backends # atlases and charts export get_point, get_point!, get_parameters, get_parameters! diff --git a/src/groups/group.jl b/src/groups/group.jl index 7b74c0edc4..75b6dc6770 100644 --- a/src/groups/group.jl +++ b/src/groups/group.jl @@ -404,7 +404,7 @@ end Compose elements ``p,q ∈ \mathcal{G}`` using the group operation ``p \circ q``. -For implementing composition on a new group manifold, please overload [`_compose`](@ref) +For implementing composition on a new group manifold, please overload `_compose` instead so that methods with [`Identity`](@ref) arguments are not ambiguous. """ compose(::AbstractGroupManifold, ::Any...) diff --git a/src/manifolds/CenteredMatrices.jl b/src/manifolds/CenteredMatrices.jl index 1f6757c2e3..33d922b727 100644 --- a/src/manifolds/CenteredMatrices.jl +++ b/src/manifolds/CenteredMatrices.jl @@ -101,7 +101,7 @@ where $c_i = \frac{1}{m}\sum_{j=1}^m p_{j,i}$ for $i = 1, \dots, n$. """ project(::CenteredMatrices, ::Any) -project!(M::CenteredMatrices, q, p) = copyto!(q, p .- mean(p, dims=1)) +project!(::CenteredMatrices, q, p) = copyto!(q, p .- mean(p, dims=1)) @doc raw""" project(M::CenteredMatrices, p, X) @@ -119,7 +119,7 @@ where $c_i = \frac{1}{m}\sum_{j=1}^m x_{j,i}$ for $i = 1, \dots, n$. """ project(::CenteredMatrices, ::Any, ::Any) -project!(M::CenteredMatrices, Y, p, X) = (Y .= X .- mean(X, dims=1)) +project!(::CenteredMatrices, Y, p, X) = (Y .= X .- mean(X, dims=1)) @generated representation_size(::CenteredMatrices{m,n,𝔽}) where {m,n,𝔽} = (m, n) diff --git a/src/manifolds/GeneralizedGrassmann.jl b/src/manifolds/GeneralizedGrassmann.jl index cd46d6b43b..41112e5dda 100644 --- a/src/manifolds/GeneralizedGrassmann.jl +++ b/src/manifolds/GeneralizedGrassmann.jl @@ -57,6 +57,39 @@ function GeneralizedGrassmann( return GeneralizedGrassmann{n,k,field,typeof(B)}(B) end +@doc raw""" + change_representer(M::GeneralizedGrassmann, ::EuclideanMetric, p, X) + +Change `X` to the corresponding representer of a cotangent vector at `p` with respect to the scaled metric +of the [`GeneralizedGrassmann`](@ref) `M`, i.e, since + +```math +g_p(X,Y) = \operatorname{tr}(Y^{\mathrm{H}}BZ) = \operatorname{tr}(X^{\mathrm{H}}Z) = ⟨X,Z⟩ +``` + +has to hold for all ``Z``, where the repreenter `X` is given, the resulting representer with +respect to the metric on the [`GeneralizedGrassmann`](@ref) is given by ``Y = B^{-1}X``. +""" +change_representer(::GeneralizedGrassmann, ::EuclideanMetric, ::Any, ::Any) + +function change_representer!(M::GeneralizedGrassmann, Y, ::EuclideanMetric, p, X) + return copyto!(M, Y, p, M.B \ X) +end + +@doc raw""" + change_metric(M::GeneralizedGrassmann, ::EuclideanMetric, p X) + +Change `X` to the corresponding vector with respect to the metric of the [`GeneralizedGrassmann`](@ref) `M`, +i.e. let ``B=LL'`` be the Cholesky decomposition of the matrix `M.B`, then the corresponding vector is ``L\X``. +""" +change_metric(M::GeneralizedGrassmann, ::EuclideanMetric, ::Any, ::Any) + +function change_metric!(M::GeneralizedGrassmann, Y, ::EuclideanMetric, p, X) + C2 = cholesky(M.B).L + Y .= C2 \ X + return Y +end + @doc raw""" check_point(M::GeneralizedGrassmann{n,k,𝔽}, p) diff --git a/src/manifolds/Hyperbolic.jl b/src/manifolds/Hyperbolic.jl index 8e44b43ad8..a9e5af25f3 100644 --- a/src/manifolds/Hyperbolic.jl +++ b/src/manifolds/Hyperbolic.jl @@ -133,7 +133,10 @@ for T in _HyperbolicTypes allocate(p::$T, ::Type{P}, dims::Tuple) where {P} = $T(allocate(p.value, P, dims)) @inline Base.copy(p::$T) = $T(copy(p.value)) - Base.copyto!(q::$T, p::$T) = copyto!(q.value, p.value) + function Base.copyto!(q::$T, p::$T) + copyto!(q.value, p.value) + return q + end Base.similar(p::$T) = $T(similar(p.value)) diff --git a/src/manifolds/HyperbolicHyperboloid.jl b/src/manifolds/HyperbolicHyperboloid.jl index e1411c99a6..798545c43f 100644 --- a/src/manifolds/HyperbolicHyperboloid.jl +++ b/src/manifolds/HyperbolicHyperboloid.jl @@ -1,3 +1,29 @@ +@doc raw""" + change_representer(M::Hyperbolic{n}, ::EuclideanMetric, p, X) + +Change the Eucliden representer `X` of a cotangent vector at point `p`. +We only have to correct for the metric, which means that the sign of the last entry changes, since +for the result ``Y`` we are looking for a tangent vector such that + +```math + g_p(Y,Z) = -y_{n+1}z_{n+1} + \sum_{i=1}^n y_iz_i = \sum_{i=1}^{n+1} z_ix_i +``` + +holds, which directly yields ``y_i=x_i`` for ``i=1,\ldots,n`` and ``y_{n+1}=-x_{n+1}``. +""" +change_representer(::Hyperbolic, ::EuclideanMetric, ::Any, ::Any) + +function change_representer!(M::Hyperbolic, Y, ::EuclideanMetric, p, X) + copyto!(M, Y, p, X) + Y[end] *= -1 + return Y +end + +function change_metric!(::Hyperbolic, ::Any, ::EuclideanMetric, ::Any, ::Any) + return error( + "Changing metric from Euclidean to Hyperbolic is not possible (see Sylvester's law of inertia).", + ) +end function check_point(M::Hyperbolic, p; kwargs...) mpv = invoke(check_point, Tuple{supertype(typeof(M)),typeof(p)}, M, p; kwargs...) @@ -325,7 +351,7 @@ component such that for the resulting `p` we have that its [`minkowski_metric`](@ref) is $⟨p,X⟩_{\mathrm{M}} = 0$, i.e. $X_{n+1} = \frac{⟨\tilde p, Y⟩}{p_{n+1}}$, where $\tilde p = (p_1,\ldots,p_n)$. """ -_hyperbolize(M::Hyperbolic, p, Y) = vcat(Y, dot(p[1:(end - 1)], Y) / p[end]) +_hyperbolize(::Hyperbolic, p, Y) = vcat(Y, dot(p[1:(end - 1)], Y) / p[end]) @doc raw""" inner(M::Hyperbolic{n}, p, X, Y) diff --git a/src/manifolds/HyperbolicPoincareBall.jl b/src/manifolds/HyperbolicPoincareBall.jl index 651e208648..0499b02f2a 100644 --- a/src/manifolds/HyperbolicPoincareBall.jl +++ b/src/manifolds/HyperbolicPoincareBall.jl @@ -1,3 +1,48 @@ +@doc raw""" + change_representer(M::Hyperbolic{n}, ::EuclideanMetric, p::PoincareBallPoint, X::PoincareBallTVector) + +Since in the metric we have the term `` α = \frac{2}{1-\sum_{i=1}^n p_i^2}`` per element, +the correction for the gradient reads `` Y = \frac{1}{α^2}X``. +""" +change_representer( + ::Hyperbolic, + ::EuclideanMetric, + ::PoincareBallPoint, + ::PoincareBallTVector, +) + +function change_representer!( + ::Hyperbolic, + Y::PoincareBallTVector, + ::EuclideanMetric, + p::PoincareBallPoint, + X::PoincareBallTVector, +) + α = 2 / (1 - norm(p.value)^2) + Y.value .= X.value ./ α^2 + return Y +end + +@doc raw""" + change_metric(M::Hyperbolic{n}, ::EuclideanMetric, p::PoincareBallPoint, X::PoincareBallTVector) + +Since in the metric we always have the term `` α = \frac{2}{1-\sum_{i=1}^n p_i^2}`` per element, +the correction for the metric reads ``Z = \frac{1}{α}X``. +""" +change_metric(::Hyperbolic, ::EuclideanMetric, ::PoincareBallPoint, ::PoincareBallTVector) + +function change_metric!( + ::Hyperbolic, + Y::PoincareBallTVector, + ::EuclideanMetric, + p::PoincareBallPoint, + X::PoincareBallTVector, +) + α = 2 / (1 - norm(p.value)^2) + Y.value .= X.value ./ α + return Y +end + function check_point(M::Hyperbolic{N}, p::PoincareBallPoint; kwargs...) where {N} mpv = check_point(Euclidean(N), p.value; kwargs...) mpv === nothing || return mpv diff --git a/src/manifolds/MetricManifold.jl b/src/manifolds/MetricManifold.jl index 066b00760c..2eb7db887b 100644 --- a/src/manifolds/MetricManifold.jl +++ b/src/manifolds/MetricManifold.jl @@ -6,11 +6,13 @@ varying inner products on the tangent space. See [`inner`](@ref). # Functor - (metric::Metric)(M::Manifold) + (metric::Metric)(M::AbstractManifold) + (metric::Metric)(M::MetricManifold) Generate the `MetricManifold` that wraps the manifold `M` with given `metric`. This works for both a variable containing the metric as well as a subtype `T<:AbstractMetric`, where a zero parameter constructor `T()` is availabe. +If `M` is already a metric manifold, the inner manifold with the new `metric` is returned. """ abstract type AbstractMetric end @@ -41,6 +43,9 @@ struct MetricManifold{𝔽,M<:AbstractManifold{𝔽},G<:AbstractMetric} <: manifold::M metric::G end +# remetricise instead of double-decorating +(metric::AbstractMetric)(M::MetricManifold) = MetricManifold(M.manifold, metric) +(::Type{T})(M::MetricManifold) where {T<:AbstractMetric} = MetricManifold(M.manifold, T()) @doc raw""" RiemannianMetric <: AbstractMetric @@ -51,6 +56,145 @@ inner product ``g(X, X) > 0`` whenever ``X`` is not the zero vector. """ abstract type RiemannianMetric <: AbstractMetric end +@doc raw""" + change_metric(M::AbstractcManifold, G2::AbstractMetric, p, X) + +On the [`AbstractManifold`](@ref) `M` with implicitly given metric ``g_1`` +and a second [`AbstractMetric`](@ref) ``g_2`` this function performs a change of metric in the +sense that it returns the tangent vector ``Z=BX`` such that the linear map ``B`` fulfills + +````math +g_2(Y_1,Y_2) = g_1(BY_1,BY_2) \quad \text{for all } Y_1, Y_2 ∈ T_p\mathcal M. +```` + +If both metrics are given in their [`local_metric`](@ref) (symmetric positive defintie) matrix +representations ``G_1 = C_1C_1^{\mathrm{H}}`` and ``G_2 = C_2C_2^{\mathrm{H}}``, where ``C_1,C_2`` denote their respective +Cholesky factors, then solving ``C_2C_2^{\mathrm{H}} = G_2 = B^{\mathrm{H}}G_1B = B^{\mathrm{H}}C_1C_1^{\mathrm{H}}B`` yields ``B = (C_1 \backslash C_2)^{\mathrm{H}}``, +where ``\cdot^{\mathrm{H}}`` denotes the conjugate transpose. + +This function returns `Z = BX`. + +# Examples + + change_metric(Sphere(2), EuclideanMetric(), p, X) + +Since the metric in ``T_p\mathbb S^2`` is the Euclidean metric from the embedding restricted to ``T_p\mathbb S^2``, this just returns `X` + + change_metric(SymmetricPOsitiveDefinite(3), EuclideanMetric, p, X) + +Here, the default metric in ``\mathcal P(3)`` is the [`LinearAffineMetric`](@ref) and the transformation can be computed as ``B=p`` +""" +change_metric(::AbstractManifold, ::AbstractMetric, ::Any, ::Any) + +function change_metric(M::AbstractManifold, G::AbstractMetric, p, X) + Y = allocate_result(M, change_metric, X, p) # this way we allocate a tangent + return change_metric!(M, Y, G, p, X) +end +function change_metric!(M::AbstractManifold, Y, G::AbstractMetric, p, X) + is_default_metric(M, G) && return copyto!(M, Y, p, X) + M.metric === G && return copyto!(M, Y, p, X) # no metric change + # TODO: For local metric, inverse_local metric, det_local_metric: Introduce a default basis? + B = DefaultOrthogonalBasis() + G1 = local_metric(M, p, B) + G2 = local_metric(G(M), p, B) + x = get_coordinates(M, p, X, B) + C1 = cholesky(G1).L + C2 = cholesky(G2).L + z = (C1 \ C2)'x + return get_vector!(M, Y, p, z, B) +end + +@decorator_transparent_signature change_metric( + M::AbstractDecoratorManifold, + G::AbstractMetric, + X, + p, +) +@decorator_transparent_signature change_metric!( + M::AbstractDecoratorManifold, + Y, + G::AbstractMetric, + X, + p, +) + +@doc raw""" + change_representer(M::AbstractManifold, G2::AbstractMetric, p, X) + +Convert the representer `X` of a linear function (in other words a cotangent vector at `p`) +in the tangent space at `p` on the [`AbstractManifold`](@ref) `M` given with respect to the +[`AbstractMetric`](@ref) `G2` into the representer with respect to the (implicit) metric of `M`. + +In order to convert `X` into the representer with respect to the (implicitly given) metric ``g_1`` of `M`, +we have to find the conversion function ``c: T_p\mathcal M \to T_p\mathcal M`` such that + +```math + g_2(X,Y) = g_1(c(X),Y) +``` + +If both metrics are given in their [`local_metric`](@ref) (symmetric positive defintie) matrix +representations ``G_1`` and ``G_2`` and ``x,y`` are the local coordinates with respect to +the same basis of the tangent space, the equation reads + +```math + x^{\mathrm{H}}G_2y = c(x)^{\mathrm{H}}G_1 y \quad \text{for all } y \in ℝ^d, +``` +where ``\cdot^{\mathrm{H}}`` denotes the conjugate transpose. +We obtain ``c(X) = (G_1\backslash G_2)^{\mathrm{H}}X`` + +For example `X` could be the gradient ``\operatorname{grad}f`` of a real-valued function +``f: \mathcal M \to ℝ``, i.e. + +```math + g_2(X,Y) = Df(p)[Y] \quad \text{for all } Y ∈ T_p\mathcal M. +``` + +and we would change the Riesz representer `X` to the representer with respect to the metric ``g_1``. + +# Examples + + change_representer(Sphere(2), EuclideanMetric(), p, X) + +Since the metric in ``T_p\mathbb S^2`` is the Euclidean metric from the embedding restricted to ``T_p\mathbb S^2``, this just returns `X` + + change_representer(SymmetricPositiveDefinite(3), EuclideanMetric(), p, X) + +Here, the default metric in ``\mathcal P(3)`` is the [`LinearAffineMetric`](@ref) and the transformation can be computed as ``pXp`` +""" +change_representer(::AbstractManifold, ::AbstractMetric, ::Any, ::Any) + +function change_representer(M::AbstractManifold, G::AbstractMetric, p, X) + Y = allocate_result(M, change_representer, X, p) # this way we allocate a tangent + return change_representer!(M, Y, G, p, X) +end + +@decorator_transparent_signature change_representer( + M::AbstractDecoratorManifold, + G::AbstractMetric, + X, + p, +) +@decorator_transparent_signature change_representer!( + M::AbstractDecoratorManifold, + Y, + G::AbstractMetric, + X, + p, +) + +# Default fallback I: compute in local metric representations +function change_representer!(M::AbstractManifold, Y, G::AbstractMetric, p, X) + is_default_metric(M, G) && return copyto!(M, Y, p, X) + M.metric === G && return copyto!(M, Y, p, X) # no metric change + # TODO: For local metric, inverse_local metric, det_local_metric: Introduce a default basis? + B = DefaultOrthogonalBasis() + G1 = local_metric(M, p, B) + G2 = local_metric(G(M), p, B) + x = get_coordinates(M, p, X, B) + z = (G1 \ G2)'x + return get_vector!(M, Y, p, z, B) +end + @doc raw""" christoffel_symbols_first( M::MetricManifold, @@ -62,9 +206,9 @@ abstract type RiemannianMetric <: AbstractMetric end Compute the Christoffel symbols of the first kind in local coordinates of basis `B`. The Christoffel symbols are (in Einstein summation convention) -```math +````math Γ_{ijk} = \frac{1}{2} \Bigl[g_{kj,i} + g_{ik,j} - g_{ij,k}\Bigr], -``` +```` where ``g_{ij,k}=\frac{∂}{∂ p^k} g_{ij}`` is the coordinate derivative of the local representation of the metric tensor. The dimensions of @@ -179,12 +323,16 @@ flat(::MetricManifold, ::Any...) end @doc raw""" - inverse_local_metric(M::AbstractcManifold, p, B::AbstractBasis) + inverse_local_metric(M::AbstractcManifold{𝔽}, p, B::AbstractBasis) -Return the local matrix representation of the inverse metric (cometric) tensor, usually -written ``g^{ij}``. +Return the local matrix representation of the inverse metric (cometric) tensor +of the tangent space at `p` on the [`AbstractManifold`](@ref) `M` with respect +to the [`AbstractBasis`](@ref) basis `B`. -See also [`local_metric`](@ref) +The metric tensor (see [`local_metric`](@ref)) is usually denoted by ``G = (g_{ij}) ∈ 𝔽^{d×d}``, +where ``d`` is the dimension of the manifold. + +Then the inverse local metric is denoted by ``G^{-1} = g^{ij}``. """ inverse_local_metric(::AbstractManifold, ::Any, ::AbstractBasis) function inverse_local_metric(M::AbstractManifold, p, B::AbstractBasis) @@ -280,13 +428,16 @@ inner(::MetricManifold, ::Any, ::Any, ::Any) end @doc raw""" - local_metric(M::AbstractManifold, p, B::AbstractBasis) + local_metric(M::AbstractManifold{𝔽}, p, B::AbstractBasis) Return the local matrix representation at the point `p` of the metric tensor ``g`` with -respect to the [`AbstractBasis`](@ref) `B` on the [`AbstractManifold`](@ref) `M`, usually written ``g_{ij}``. -The matrix has the property that ``g(X, Y)=X^\mathrm{T} [g_{ij}] Y = g_{ij} X^i Y^j``, -where the latter expression uses Einstein summation convention. -The metric tensor is such that the formula works for the given [`AbstractBasis`](@ref) `B`. +respect to the [`AbstractBasis`](@ref) `B` on the [`AbstractManifold`](@ref) `M`. +Let ``d``denote the dimension of the manifold and $b_1,\ldots,b_d$ the basis vectors. +Then the local matrix representation is a matrix ``G\in 𝔽^{n\times n}`` whose entries are +given by ``g_{ij} = g_p(b_i,b_j), i,j\in\{1,…,d\}``. + +This yields the property for two tangent vectors (using Einstein summation convention) +``X = X^ib_i, Y=Y^ib_i \in T_p\mathcal M`` we get ``g_p(X, Y) = g_{ij} X^i Y^j``. """ local_metric(::AbstractManifold, ::Any, ::AbstractBasis) @decorator_transparent_signature local_metric( @@ -436,7 +587,7 @@ for f in [ quote function decorator_transparent_dispatch( ::typeof($f), - ::AbstractConnectionManifold, + M::AbstractConnectionManifold, args..., ) return Val(:parent) @@ -445,6 +596,19 @@ for f in [ ) end +for f in [change_metric, change_representer, change_metric!, change_representer!] + eval( + quote + function decorator_transparent_dispatch( + ::typeof($f), + ::AbstractManifold, + args..., + ) + return Val(:parent) + end + end, + ) +end function decorator_transparent_dispatch( ::typeof(christoffel_symbols_second), ::MetricManifold, diff --git a/src/manifolds/Multinomial.jl b/src/manifolds/Multinomial.jl index 36f86eef64..8db4e6e63d 100644 --- a/src/manifolds/Multinomial.jl +++ b/src/manifolds/Multinomial.jl @@ -53,6 +53,7 @@ function check_point(M::MultinomialMatrices{n,m}, p; kwargs...) where {n,m} end return check_point(PowerManifold(M.manifold, m), p; kwargs...) end + @doc raw""" check_vector(M::MultinomialMatrices p, X; kwargs...) diff --git a/src/manifolds/PositiveNumbers.jl b/src/manifolds/PositiveNumbers.jl index 05f16c3c05..c3b22362fa 100644 --- a/src/manifolds/PositiveNumbers.jl +++ b/src/manifolds/PositiveNumbers.jl @@ -38,6 +38,55 @@ This manifold is modeled as a [`PowerManifold`](@ref) of [`PositiveNumbers`](@re """ PositiveArrays(n::Vararg{Int,I}) where {I} = PositiveNumbers()^(n) +@doc raw""" + change_representer(M::PositiveNumbers, E::EuclideanMetric, p, X) + +Given a tangent vector ``X ∈ T_p\mathcal M`` representing a linear function with respect +to the [`EuclideanMetric`](@ref) `g_E`, this function changes the representer into the one +with respect to the positivity metric representation of +[`PositiveNumbers`](@ref) `M`. + +For all tangent vectors ``Y`` the result ``Z`` has to fulfill + +```math + ⟨X,Y⟩ = XY = \frac{ZY}{p^2} = g_p(YZ) +``` + +and hence ``Z = p^2X`` + +""" +change_representer(::PositiveNumbers, ::EuclideanMetric, ::Any, ::Any) +change_representer(::PositiveNumbers, ::EuclideanMetric, p::Real, X::Real) = p * X * p + +function change_representer!(::PositiveNumbers, Y, ::EuclideanMetric, p, X) + Y .= p .* X .* p + return Y +end + +@doc raw""" + change_metric(M::PositiveNumbers, E::EuclideanMetric, p, X) + +Given a tangent vector ``X ∈ T_p\mathcal M`` representing a linear function with respect to +the [`EuclideanMetric`](@ref) `g_E`, +this function changes the representer into the one with respect to the positivity metric +of [`PositiveNumbers`](@ref) `M`. + +For all ``Z,Y`` we are looking for the function ``c`` on the tangent space at ``p`` such that + +```math + ⟨Z,Y⟩ = XY = \frac{c(Z)c(Y)}{p^2} = g_p(c(Y),c(Z)) +``` + +and hence ``C(X) = pX``. +""" +change_metric(::PositiveNumbers, ::EuclideanMetric, ::Any, ::Any) +change_metric(::PositiveNumbers, ::EuclideanMetric, p::Real, X::Real) = p * X + +function change_metric!(::PositiveNumbers, Y, ::EuclideanMetric, p, X) + Y .= p .* X + return Y +end + @doc raw""" check_point(M::PositiveNumbers, p) diff --git a/src/manifolds/PowerManifold.jl b/src/manifolds/PowerManifold.jl index 833a952d04..e5cc8c6ee8 100644 --- a/src/manifolds/PowerManifold.jl +++ b/src/manifolds/PowerManifold.jl @@ -78,6 +78,48 @@ end default_metric_dispatch(::AbstractPowerManifold, ::PowerMetric) = Val(true) +""" + change_representer(M::AbstractPowerManifold, ::AbstractMetric, p, X) + +Since the metric on a power manifold decouples, the change of a representer can be done elementwise +""" +change_representer(::AbstractPowerManifold, ::AbstractMetric, ::Any, ::Any) + +function change_representer!(M::AbstractPowerManifold, Y, G::AbstractMetric, p, X) + rep_size = representation_size(M.manifold) + for i in get_iterator(M) + change_representer!( + M.manifold, + _write(M, rep_size, Y, i), + G, + _read(M, rep_size, p, i), + _read(M, rep_size, X, i), + ) + end + return Y +end + +""" + change_metric(M::AbstractPowerManifold, ::AbstractMetric, p, X) + +Since the metric on a power manifold decouples, the change of metric can be done elementwise. +""" +change_metric(M::AbstractPowerManifold, ::AbstractMetric, ::Any, ::Any) + +function change_metric!(M::AbstractPowerManifold, Y, G::AbstractMetric, p, X) + rep_size = representation_size(M.manifold) + for i in get_iterator(M) + change_metric!( + M.manifold, + _write(M, rep_size, Y, i), + G, + _read(M, rep_size, p, i), + _read(M, rep_size, X, i), + ) + end + return Y +end + @doc raw""" flat(M::AbstractPowerManifold, p, X) diff --git a/src/manifolds/ProbabilitySimplex.jl b/src/manifolds/ProbabilitySimplex.jl index 557c344c25..c66aebff99 100644 --- a/src/manifolds/ProbabilitySimplex.jl +++ b/src/manifolds/ProbabilitySimplex.jl @@ -58,6 +58,44 @@ See for example the [`ProbabilitySimplex`](@ref). """ struct FisherRaoMetric <: AbstractMetric end +@doc raw""" + change_representer(M::ProbabilitySimplex, ::EuclideanMetric, p, X) + +Given a tangent vector with respect to the metric from the embedding, the [`EuclideanMetric`](@ref), +the representer of a linear functional on the tangent space is adapted as ``Z = p .* X``, since +this “compensates” for the divsion by ``p`` in the Riemannian metric on the [`ProbabilitySimplex`](@ref). + +To be precise for any ``Y ∈ T_pΔ^n`` we are looking for ``Z ∈ T_pΔ^n`` such that + +```math + ⟨X,Y⟩ = X^\mathrm{T}Y = \sum_{i=1}^{n+1}\frac{Z_iY_i}{p_i} = g_p(Z,Y) +``` + +and hence ``Z_i = X_ip_i, i=1,…,n+1``. +""" +change_representer(::ProbabilitySimplex, ::EuclideanMetric, ::Any, ::Any) + +function change_representer!(::ProbabilitySimplex, Y, ::EuclideanMetric, p, X) + return Y .= p .* X +end + +@doc raw""" + change_metric(M::ProbabilitySimplex, ::EuclideanMetric, p, X) + +To change the metric, we are looking for a function ``c\colon T_pΔ^n \to T_pΔ^n`` such that for all ``X,Y ∈ T_pΔ^n`` + +```math + ⟨X,Y⟩ = X^\mathrm{T}Y = \sum_{i=1}^{n+1}\frac{c(X)_ic(Y)_i}{p_i} = g_p(X,Y) +``` + +and hence ``C(X)_i = X_i\sqrt{p_i}, i=1,…,n+1``. +""" +change_metric(::ProbabilitySimplex, ::EuclideanMetric, ::Any, ::Any) + +function change_metric!(::ProbabilitySimplex, Y, ::EuclideanMetric, p, X) + return Y .= sqrt.(p) .* X +end + """ check_point(M::ProbabilitySimplex, p; kwargs...) diff --git a/src/manifolds/ProductManifold.jl b/src/manifolds/ProductManifold.jl index 4aa0e60833..aa26540412 100644 --- a/src/manifolds/ProductManifold.jl +++ b/src/manifolds/ProductManifold.jl @@ -135,6 +135,42 @@ function ProductVectorTransport(methods::AbstractVectorTransportMethod...) return ProductVectorTransport{typeof(methods)}(methods) end +""" + change_representer(M::ProductManifold, ::AbstractMetric, p, X) + +Since the metric on a product manifold decouples, the change of a representer can be done elementwise +""" +change_representer(::ProductManifold, ::AbstractMetric, ::Any, ::Any) + +function change_representer!(M::ProductManifold, Y, G::AbstractMetric, p, X) + map( + (m, y, P, x) -> change_representer!(m, y, G, P, x), + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + ) + return Y +end + +""" + change_metric(M::ProductManifold, ::AbstractMetric, p, X) + +Since the metric on a product manifold decouples, the change of metric can be done elementwise. +""" +change_metric(::ProductManifold, ::AbstractMetric, ::Any, ::Any) + +function change_metric!(M::ProductManifold, Y, G::AbstractMetric, p, X) + map( + (m, y, P, x) -> change_metric!(m, y, G, P, x), + M.manifolds, + submanifold_components(M, Y), + submanifold_components(M, p), + submanifold_components(M, X), + ) + return Y +end + """ check_point(M::ProductManifold, p; kwargs...) diff --git a/src/manifolds/Spectrahedron.jl b/src/manifolds/Spectrahedron.jl index 53289fccc8..efd14f7f12 100644 --- a/src/manifolds/Spectrahedron.jl +++ b/src/manifolds/Spectrahedron.jl @@ -138,7 +138,7 @@ project!(::Spectrahedron, r, q) = copyto!(r, q ./ norm(q)) """ project(M::Spectrahedron, q, Y) -Project `Y` onto the tangent space at `q`, i.e. row-wise onto the oblique manifold. +Project `Y` onto the tangent space at `q`, i.e. row-wise onto the Spectrahedron manifold. """ project(::Spectrahedron, ::Any...) diff --git a/src/manifolds/Sphere.jl b/src/manifolds/Sphere.jl index 594e05bc6c..92f07d1a25 100644 --- a/src/manifolds/Sphere.jl +++ b/src/manifolds/Sphere.jl @@ -144,7 +144,9 @@ end function decorated_manifold(M::AbstractSphere{𝔽}) where {𝔽} return Euclidean(representation_size(M)...; field=𝔽) end -get_embedding(M::AbstractSphere{𝔽}) where {𝔽} = decorated_manifold(M) + +# Since on every tangent space the Euclidean matric (restricted to this space) is used, this should be fine +default_metric_dispatch(::AbstractSphere, ::EuclideanMetric) = Val(true) @doc raw""" distance(M::AbstractSphere, p, q) @@ -227,6 +229,8 @@ function get_coordinates!( return Y end +get_embedding(M::AbstractSphere{𝔽}) where {𝔽} = decorated_manifold(M) + @doc raw""" get_vector(M::AbstractSphere{ℝ}, p, X, B::DefaultOrthonormalBasis) diff --git a/src/manifolds/SymmetricPositiveDefiniteLinearAffine.jl b/src/manifolds/SymmetricPositiveDefiniteLinearAffine.jl index 237dce32ce..e072a86676 100644 --- a/src/manifolds/SymmetricPositiveDefiniteLinearAffine.jl +++ b/src/manifolds/SymmetricPositiveDefiniteLinearAffine.jl @@ -6,6 +6,53 @@ matrix logarithms and exponentials, which yields a linear and affine metric. """ struct LinearAffineMetric <: RiemannianMetric end +@doc raw""" + change_representer(M::SymmetricPositiveDefinite, 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`](@ref) `g_E`, +this is turned into the representer with respect to the (default) metric, +the [`LinearAffineMetric`](@ref) on the [`SymmetricPositiveDefinite`](@ref) `M`. + +To be precise we are looking for ``Z∈T_p\mathcal P(n)`` such that for all ``Y∈T_p\mathcal P(n)``` +it holds + +```math +⟨X,Y⟩ = \operatorname{tr}(XY) = \operatorname{tr}(p^{-1}Zp^{-1}Y) = g_p(Z,Y) +``` + +and hence ``Z = pXp``. +""" +change_representer(::SymmetricPositiveDefinite, ::EuclideanMetric, ::Any, ::Any) + +function change_representer!(::SymmetricPositiveDefinite, Y, ::EuclideanMetric, p, X) + Y .= p * X * p + return Y +end + +@doc raw""" + change_metric(M::SymmetricPositiveDefinite{n}, E::EuclideanMetric, p, X) + +Given a tangent vector ``X ∈ T_p\mathcal P(n)`` with respect to the [`EuclideanMetric`](@ref) `g_E`, +this function changes into the [`LinearAffineMetric`](@ref) (default) metric on the +[`SymmetricPositiveDefinite`](@ref) `M`. + +To be precise we are looking for ``c\colon T_p\mathcal P(n) \to T_p\mathcal P(n) `` +such that for all ``Y,Z ∈ T_p\mathcal P(n)``` it holds + +```math +⟨Y,Z⟩ = \operatorname{tr}(YZ) = \operatorname{tr}(p^{-1}c(Y)p^{-1}c(Z)) = g_p(c(Z),c(Y)) +``` + +and hence ``c(X) = pX`` is computed. +""" +change_metric(::SymmetricPositiveDefinite, ::EuclideanMetric, ::Any, ::Any) + +function change_metric!(::SymmetricPositiveDefinite, Y, ::EuclideanMetric, p, X) + Y .= p * X + return Y +end + default_metric_dispatch(::SymmetricPositiveDefinite, ::LinearAffineMetric) = Val(true) @doc raw""" diff --git a/src/riemannian_diff.jl b/src/riemannian_diff.jl index de8654bf98..d3cc0d1bea 100644 --- a/src/riemannian_diff.jl +++ b/src/riemannian_diff.jl @@ -232,7 +232,7 @@ end function gradient(M::AbstractManifold, f, p, backend::RiemannianProjectionGradientBackend) amb_grad = _gradient(f, p, backend.diff_backend) - return project(M, p, amb_grad) + return change_representer(M, EuclideanMetric(), p, project(M, p, amb_grad)) end function gradient!( @@ -244,5 +244,6 @@ function gradient!( ) amb_grad = embed(M, p, X) _gradient!(f, amb_grad, p, backend.diff_backend) - return project!(M, X, p, amb_grad) + project!(M, X, p, amb_grad) + return change_representer!(M, X, EuclideanMetric(), p, X) end diff --git a/test/manifolds/generalized_grassmann.jl b/test/manifolds/generalized_grassmann.jl index 1e81a7ca41..b8c6905321 100644 --- a/test/manifolds/generalized_grassmann.jl +++ b/test/manifolds/generalized_grassmann.jl @@ -4,7 +4,7 @@ include("../utils.jl") @testset "Real" begin B = [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0] M = GeneralizedGrassmann(3, 2, B) - x = [1.0 0.0; 0.0 0.5; 0.0 0.0] + p = [1.0 0.0; 0.0 0.5; 0.0 0.0] @testset "Basics" begin @test repr(M) == "GeneralizedGrassmann(3, 2, [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0], ℝ)" @@ -13,62 +13,69 @@ include("../utils.jl") @test base_manifold(M) === M @test_throws DomainError is_point(M, [1.0, 0.0, 0.0, 0.0], true) @test_throws DomainError is_point(M, 1im * [1.0 0.0; 0.0 1.0; 0.0 0.0], true) - @test !is_vector(M, x, [0.0, 0.0, 1.0, 0.0]) - @test_throws DomainError is_vector(M, x, [0.0, 0.0, 1.0, 0.0], true) - @test_throws DomainError is_vector(M, x, 1 * im * zero_vector(M, x), true) + @test !is_vector(M, p, [0.0, 0.0, 1.0, 0.0]) + @test_throws DomainError is_vector(M, p, [0.0, 0.0, 1.0, 0.0], true) + @test_throws DomainError is_vector(M, p, 1 * im * zero_vector(M, p), true) @test injectivity_radius(M) == π / 2 @test injectivity_radius(M, ExponentialRetraction()) == π / 2 - @test injectivity_radius(M, x) == π / 2 - @test injectivity_radius(M, x, ExponentialRetraction()) == π / 2 - @test mean(M, [x, x, x]) == x + @test injectivity_radius(M, p) == π / 2 + @test injectivity_radius(M, p, ExponentialRetraction()) == π / 2 + @test mean(M, [p, p, p]) == p end @testset "Embedding and Projection" begin - y = similar(x) - z = embed(M, x) - @test z == x - embed!(M, y, x) - @test y == z + q = similar(p) + p2 = embed(M, p) + @test p2 == p + embed!(M, q, p) + @test q == p2 a = [1.0 0.0; 0.0 2.0; 0.0 0.0] @test !is_point(M, a) b = similar(a) c = project(M, a) - @test c == x + @test c == p project!(M, b, a) - @test b == x + @test b == p X = [0.0 0.0; 0.0 0.0; -1.0 1.0] Y = similar(X) - Z = embed(M, x, X) - embed!(M, Y, x, X) + Z = embed(M, p, X) + embed!(M, Y, p, X) @test Y == X @test Z == X end - + @testset "gradient and metric conversion" begin + L = cholesky(B).L + X = [0.0 0.0; 0.0 0.0; 1.0 -1.0] + Y = change_metric(M, EuclideanMetric(), p, X) + @test Y == L \ X + Z = change_representer(M, EuclideanMetric(), p, X) + @test Z == B \ X + end types = [Matrix{Float64}] TEST_STATIC_SIZED && push!(types, MMatrix{3,2,Float64,6}) X = [0.0 0.0; 1.0 0.0; 0.0 2.0] Y = [0.0 0.0; -1.0 0.0; 0.0 2.0] - @test inner(M, x, X, Y) == 0 - y = retract(M, x, X) - z = retract(M, x, Y) - @test is_point(M, y) - @test is_point(M, z) - @test retract(M, x, X) == exp(M, x, X) + @test inner(M, p, X, Y) == 0 + q = retract(M, p, X) + r = retract(M, p, Y) + @test is_point(M, q) + @test is_point(M, r) + @test retract(M, p, X) == exp(M, p, X) - a = project(M, x + X) - c = retract(M, x, X, ProjectionRetraction()) - d = retract(M, x, X, PolarRetraction()) + a = project(M, p + X) + c = retract(M, p, X, ProjectionRetraction()) + d = retract(M, p, X, PolarRetraction()) @test a == c @test c == d e = similar(a) - retract!(M, e, x, X) - @test e == exp(M, x, X) - @test vector_transport_to(M, x, X, y, ProjectionTransport()) == project(M, y, X) + retract!(M, e, p, X) + @test e == exp(M, p, X) + @test vector_transport_to(M, p, X, q, ProjectionTransport()) == project(M, q, X) @testset "Type $T" for T in types - pts = convert.(T, [x, y, z]) - @test !is_point(M, 2 * x) - @test_throws DomainError !is_point(M, 2 * x, true) - @test !is_vector(M, x, y) - @test_throws DomainError is_vector(M, x, y, true) + pts = convert.(T, [p, q, r]) + @test !is_point(M, 2 * p) + @test_throws DomainError !is_point(M, 2 * r, true) + @test !is_vector(M, p, q) + @test_throws DomainError is_vector(M, p, q, true) test_manifold( M, pts, diff --git a/test/manifolds/hyperbolic.jl b/test/manifolds/hyperbolic.jl index af1272510e..c0caa4782c 100644 --- a/test/manifolds/hyperbolic.jl +++ b/test/manifolds/hyperbolic.jl @@ -44,8 +44,9 @@ include("../utils.jl") copyto!(M, pC, p) @test pC.value == p.value XC = allocate(X) - copyto!(M, XC, p, X) - @test XC.value == X.value + @test copyto!(M, XC, p, X) == X # does copyto return the right value? + @test XC == X # does copyto store the right value? + @test XC.value == X.value # another check end end @testset "Hyperbolic Representation Conversion I" begin @@ -261,4 +262,28 @@ include("../utils.jl") @test isa(Z2, PoincareHalfSpaceTVector) @test Z2 == X2 end + @testset "Metric conversion on Hyperboloid" begin + M = Hyperbolic(2) + p = [1.0, 1.0, sqrt(3)] + X = [1.0, 2.0, sqrt(3)] + Y = change_representer(M, EuclideanMetric(), p, X) + @test inner(M, p, X, Y) == inner(Euclidean(3), p, X, X) + # change metric not possible from Euclidean, since the embedding is Lorenzian + @test_throws ErrorException change_metric(M, EuclideanMetric(), p, X) + # but if we come from the same metric, we have the identity + @test change_metric(M, MinkowskiMetric(), p, X) == X + end + @testset "Metric conversion on Poincare Ball" begin + M = Hyperbolic(2) + p = convert(PoincareBallPoint, [1.0, 1.0, sqrt(3)]) + X = convert(PoincareBallTVector, [1.0, 1.0, sqrt(3)], [1.0, 2.0, sqrt(3)]) + Y = change_representer(M, EuclideanMetric(), p, X) + @test inner(M, p, X, Y) == inner(Euclidean(3), p, X.value, X.value) + α = 2 / (1 - norm(p.value)^2) + @test Y.value == X.value ./ α^2 + Z = change_metric(M, EuclideanMetric(), p, X) + @test Z.value == X.value ./ α + A = change_metric(M, MinkowskiMetric(), p, X) + @test A == X + end end diff --git a/test/manifolds/positive_numbers.jl b/test/manifolds/positive_numbers.jl index 6a37883a7e..3e4b00ee1a 100644 --- a/test/manifolds/positive_numbers.jl +++ b/test/manifolds/positive_numbers.jl @@ -23,6 +23,12 @@ include("../utils.jl") X = similar([1.0]) zero_vector!(M, X, 1.0) @test X == [0.0] + + @test change_metric(M, EuclideanMetric(), 2, 3) == 3 * 2 + @test change_representer(M, EuclideanMetric(), 2, 3) == 3 * 2^2 + N = PositiveVectors(2) + @test change_metric(M, EuclideanMetric(), [1, 2], [3, 4]) == [3, 4 * 2] + @test change_representer(M, EuclideanMetric(), [1, 2], [3, 4]) == [3, 4 * 2^2] end types = [Float64] TEST_FLOAT32 && push!(types, Float32) diff --git a/test/manifolds/power_manifold.jl b/test/manifolds/power_manifold.jl index bc7d03b95b..61ad2a633c 100644 --- a/test/manifolds/power_manifold.jl +++ b/test/manifolds/power_manifold.jl @@ -413,4 +413,23 @@ end @test isapprox(a, a2) @test_throws ErrorException get_point(M, A2, p, a2) end + + @testset "metric conversion" begin + M = SymmetricPositiveDefinite(3) + N = PowerManifold(M, NestedPowerRepresentation(), 2) + e = EuclideanMetric() + p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1] + q = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1] + P = [p, q] + X = [log(M, p, q), log(M, q, p)] + Y = change_metric(N, e, P, X) + Yc = [change_metric(M, e, p, log(M, p, q)), change_metric(M, e, q, log(M, q, p))] + @test norm(N, P, Y .- Yc) ≈ 0 + Z = change_representer(N, e, P, X) + Zc = [ + change_representer(M, e, p, log(M, p, q)), + change_representer(M, e, q, log(M, q, p)), + ] + @test norm(N, P, Z .- Zc) ≈ 0 + end end diff --git a/test/manifolds/probability_simplex.jl b/test/manifolds/probability_simplex.jl index 071903636b..5d2e09dc70 100644 --- a/test/manifolds/probability_simplex.jl +++ b/test/manifolds/probability_simplex.jl @@ -72,6 +72,12 @@ include("../utils.jl") Y = project(M, p, X2) @test isapprox(M, p, X, Y) + # Check adaption of metric and representer + Y1 = change_metric(M, EuclideanMetric(), p, X) + @test Y1 == X .* sqrt.(p) + Y2 = change_representer(M, EuclideanMetric(), p, X) + @test Y2 == X .* p + X = log(M, q, p) X2 = X + [1, 2, 3] Y = project(M, q, X2) diff --git a/test/manifolds/product_manifold.jl b/test/manifolds/product_manifold.jl index 5d9a742174..07dc94e807 100644 --- a/test/manifolds/product_manifold.jl +++ b/test/manifolds/product_manifold.jl @@ -693,4 +693,26 @@ end p2 = get_point(M, A, p, a) @test all(p2.parts .== p.parts) end + + @testset "metric conversion" begin + M = SymmetricPositiveDefinite(3) + N = ProductManifold(M, M) + e = EuclideanMetric() + p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1] + q = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1] + P = ProductRepr(p, q) + X = ProductRepr(log(M, p, q), log(M, q, p)) + Y = change_metric(N, e, P, X) + Yc = ProductRepr( + change_metric(M, e, p, log(M, p, q)), + change_metric(M, e, q, log(M, q, p)), + ) + @test norm(N, P, Y - Yc) ≈ 0 + Z = change_representer(N, e, P, X) + Zc = ProductRepr( + change_representer(M, e, p, log(M, p, q)), + change_representer(M, e, q, log(M, q, p)), + ) + @test norm(N, P, Z - Zc) ≈ 0 + end end diff --git a/test/manifolds/sphere.jl b/test/manifolds/sphere.jl index ff533f8bf1..f5889e9c26 100644 --- a/test/manifolds/sphere.jl +++ b/test/manifolds/sphere.jl @@ -214,4 +214,13 @@ using ManifoldsBase: TFVector end end end + + @testset "Metric conversion is the identity" begin + p = [1.0, 0.0, 0.0] + X = [0.0, 1.0, 1.0] + Y = change_representer(M, EuclideanMetric(), p, X) + @test Y == X + Z = change_metric(M, EuclideanMetric(), p, X) + @test Z == X + end end diff --git a/test/manifolds/symmetric_positive_definite.jl b/test/manifolds/symmetric_positive_definite.jl index eeeb447cb1..0b134d378e 100644 --- a/test/manifolds/symmetric_positive_definite.jl +++ b/test/manifolds/symmetric_positive_definite.jl @@ -87,35 +87,35 @@ using Manifolds: default_metric_dispatch end end end - x = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1] - y = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1] + p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1] + q = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1] @testset "Convert SPD to Cholesky" begin - v = log(M1, x, y) - (l, w) = Manifolds.spd_to_cholesky(x, v) + v = log(M1, p, q) + (l, w) = Manifolds.spd_to_cholesky(p, v) (xs, vs) = Manifolds.cholesky_to_spd(l, w) - @test isapprox(xs, x) + @test isapprox(xs, p) @test isapprox(vs, v) end @testset "Preliminary tests for LogEuclidean" begin @test representation_size(M4) == (3, 3) - @test isapprox(distance(M4, x, y), sqrt(2) * log(2)) + @test isapprox(distance(M4, p, q), sqrt(2) * log(2)) @test manifold_dimension(M4) == manifold_dimension(M1) end @testset "Test for tangent ONB on LinearAffineMetric" begin - v = log(M2, x, y) - donb = get_basis(base_manifold(M2), x, DiagonalizingOrthonormalBasis(v)) - X = get_vectors(base_manifold(M2), x, donb) + v = log(M2, p, q) + donb = get_basis(base_manifold(M2), p, DiagonalizingOrthonormalBasis(v)) + X = get_vectors(base_manifold(M2), p, donb) k = donb.data.eigenvalues @test isapprox(0.0, first(k)) for i in 1:length(X) - @test isapprox(1.0, norm(M2, x, X[i])) + @test isapprox(1.0, norm(M2, p, X[i])) for j in (i + 1):length(X) - @test isapprox(0.0, inner(M2, x, X[i], X[j])) + @test isapprox(0.0, inner(M2, p, X[i], X[j])) end end - d2onb = get_basis(M2, x, DiagonalizingOrthonormalBasis(v)) + d2onb = get_basis(M2, p, DiagonalizingOrthonormalBasis(v)) @test donb.data.eigenvalues == d2onb.data.eigenvalues - @test get_vectors(base_manifold(M2), x, donb) == get_vectors(M2, x, d2onb) + @test get_vectors(base_manifold(M2), p, donb) == get_vectors(M2, p, d2onb) end @testset "Vector transport and transport along with Schild and Pole ladder" begin A(α) = [1.0 0.0 0.0; 0.0 cos(α) sin(α); 0.0 -sin(α) cos(α)] @@ -141,4 +141,11 @@ using Manifolds: default_metric_dispatch @test inner(M, p1, X1, X2) ≈ inner(M, p2, Y1, Y4) # parallel transport isometric @test inner(M, p1, X1, X2) ≈ inner(M, p2, Y2, Y5) # pole ladder transport isometric end + @testset "Metric change for Linear Affine Metric" begin + X = log(M1, p, q) + Y = change_metric(M1, EuclideanMetric(), p, X) + @test Y == p * X + Z = change_representer(M1, EuclideanMetric(), p, X) + @test Z == p * X * p + end end diff --git a/test/metric.jl b/test/metric.jl index ae7d7bb6d6..a1562bfe2c 100644 --- a/test/metric.jl +++ b/test/metric.jl @@ -7,6 +7,7 @@ include("utils.jl") struct TestEuclidean{N} <: AbstractManifold{ℝ} end struct TestEuclideanMetric <: AbstractMetric end +struct TestScaledEuclideanMetric <: AbstractMetric end Manifolds.manifold_dimension(::TestEuclidean{N}) where {N} = N function Manifolds.local_metric( @@ -19,11 +20,57 @@ end function Manifolds.local_metric( M::MetricManifold{ℝ,<:TestEuclidean,<:TestEuclideanMetric}, ::Any, - ::InducedBasis, -) + ::T, +) where {T<:ManifoldsBase.AbstractOrthogonalBasis} return Diagonal(1.0:manifold_dimension(M)) end - +function Manifolds.local_metric( + M::MetricManifold{ℝ,<:TestEuclidean,<:TestScaledEuclideanMetric}, + ::Any, + ::T, +) where {T<:ManifoldsBase.AbstractOrthogonalBasis} + return 2 .* Diagonal(1.0:manifold_dimension(M)) +end +function Manifolds.get_coordinates!( + M::MetricManifold{ℝ,<:TestEuclidean,<:TestEuclideanMetric}, + c, + ::Any, + X, + ::DefaultOrthogonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, +) + c .= 1 ./ [1.0:manifold_dimension(M)...] .* X + return c +end +function Manifolds.get_vector!( + M::MetricManifold{ℝ,<:TestEuclidean,<:TestEuclideanMetric}, + X, + ::Any, + c, + ::DefaultOrthogonalBasis{ℝ,<:ManifoldsBase.TangentSpaceType}, +) + X .= [1.0:manifold_dimension(M)...] .* c + return X +end +function Manifolds.get_coordinates!( + M::MetricManifold{ℝ,<:TestEuclidean,<:TestScaledEuclideanMetric}, + c, + ::Any, + X, + ::DefaultOrthogonalBasis, +) + c .= 1 ./ (2 .* [1.0:manifold_dimension(M)...]) .* X + return c +end +function Manifolds.get_vector!( + M::MetricManifold{ℝ,<:TestEuclidean,<:TestScaledEuclideanMetric}, + X, + ::Any, + c, + ::DefaultOrthogonalBasis, +) + X .= 2 .* [1.0:manifold_dimension(M)...] .* c + return X +end struct TestSphere{N,T} <: AbstractManifold{ℝ} r::T end @@ -603,4 +650,20 @@ end ExponentialRetraction(), ) === Val{:parent}() end + + @testset "change metric and representer" begin + M = MetricManifold(TestEuclidean{2}(), TestEuclideanMetric()) + G = TestScaledEuclideanMetric() + M2 = TestScaledEuclideanMetric(M) + @test M2.manifold === M.manifold + @test M2.metric == G + p = ones(2) + X = 2 * ones(2) + @test change_metric(M, TestEuclideanMetric(), p, X) == X + Y = change_metric(M, G, p, X) + @test Y ≈ sqrt(2) .* X #scaled metric has a factor 2, removing introduces this factor + @test change_representer(M, TestEuclideanMetric(), p, X) == X + Y2 = change_representer(M, G, p, X) + @test Y2 ≈ 2 .* X #scaled metric has a factor 2, removing introduces this factor + end end