Introduce r-norms on manifolds with components

NEWS.md

@@ -5,6 +5,15 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.15.18] 18/10/2024
+
+### Added
+
+* `distance(M, p, q, r)` to compute `r`-norms on manifolds that have components.
+* `distance(M, p, q, m, r)` to compute (approximate) `r`-norms on manifolds that have components
+  using an `AbstractInverseRetractionMethod m` within every (inner) distance call.
+* `norm(M, p, X, r)` to compute `r`-norms on manifolds that have components.
+
 ## [0.15.17] 04/10/2024
 
 ### Changed

Project.toml

@@ -1,7 +1,7 @@
 name = "ManifoldsBase"
 uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
 authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>", "Antoine Levitt <[email protected]>"]
-version = "0.15.17"
+version = "0.15.18"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/DefaultManifold.jl

@@ -57,7 +57,7 @@ function check_approx(M::DefaultManifold, p, X, Y; kwargs...)
     return ApproximatelyError(v, s)
 end
 
-distance(::DefaultManifold, p, q) = norm(p - q)
+distance(::DefaultManifold, p, q, r::Real = 2.0) = norm(p - q, r)
 
 embed!(::DefaultManifold, q, p) = copyto!(q, p)
 
@@ -123,6 +123,8 @@ function get_vector_orthonormal!(M::DefaultManifold{ℂ}, Y, p, c, ::RealNumbers
     return copyto!(Y, reshape(c[1:n] + c[(n + 1):(2n)] * 1im, representation_size(M)))
 end
 
+has_components(::DefaultManifold) = true
+
 injectivity_radius(::DefaultManifold) = Inf
 
 @inline inner(::DefaultManifold, p, X, Y) = dot(X, Y)
@@ -138,7 +140,7 @@ end
 
 number_system(::DefaultManifold{𝔽}) where {𝔽} = 𝔽
 
-norm(::DefaultManifold, p, X) = norm(X)
+norm(::DefaultManifold, p, X, r = 2.0) = norm(X, r)
 
 project!(::DefaultManifold, q, p) = copyto!(q, p)
 project!(::DefaultManifold, Y, p, X) = copyto!(Y, X)

src/ManifoldsBase.jl

@@ -594,6 +594,15 @@ function embed_project!(M::AbstractManifold, Y, p, X)
     return project!(M, Y, p, embed(M, p, X))
 end
 
+"""
+    has_components(M::AbstractManifold)
+
+Return whether the [`AbstractManifold`](@ref)`(M)` consists of components,
+like the [`PowerManifold`](@ref) or the [`ProductManifold`](@ref), that one can iterate over.
+By default, this function returns `false`.
+"""
+has_components(M::AbstractManifold) = false
+
 @doc raw"""
     injectivity_radius(M::AbstractManifold)
 
@@ -1354,6 +1363,7 @@ export ×,
     get_vector,
     get_vector!,
     get_vectors,
+    has_components,
     hat,
     hat!,
     injectivity_radius,

src/PowerManifold.jl

@@ -530,20 +530,103 @@
     return default_vector_transport_method(M.manifold, eltype(t))
 end
 
-@doc raw"""
-    distance(M::AbstractPowerManifold, p, q)
+_doc_distance_pow = """
+    distance(M::AbstractPowerManifold, p, q, r::Real=2.0)
+    distance(M::AbstractPowerManifold, p, q, m::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(), r::Real=2.0)
+
+Compute the distance between `q` and `p` on an [`AbstractPowerManifold`](@ref).
 
-Compute the distance between `q` and `p` on an [`AbstractPowerManifold`](@ref),
-i.e. from the element wise distances the Forbenius norm is computed.
+First, the componentwise distances are computed. These can be approximated using the
+`norm` of an [`AbstractInverseRetractionMethod`](@ref) `m`.
+Then, the `r`-norm of these elements is computed.
 """
+
 function distance(M::AbstractPowerManifold, p, q)
-    sum_squares = zero(number_eltype(p))
+    return _distance_r(M, p, q, 2.0)
+end
+
+@doc "$(_doc_distance_pow)"
+function distance(M::AbstractPowerManifold, p, q, r::Real)
+    isinf(r) && return _distance_Inf(M, p, q)
+    (r == 1) && return _distance_1(M, p, q)
+    return _distance_r(M, p, q, r)
+end
+function distance(
+    M::AbstractPowerManifold,
+    p,
+    q,
+    ::LogarithmicInverseRetraction,
+    r::Real = 2.0,
+)
+    return distance(M, p, q, r)
+end
+
+@doc "$(_doc_distance_pow)"
+function distance(
+    M::AbstractPowerManifold,
+    p,
+    q,
+    m::AbstractInverseRetractionMethod,
+    r::Real = 2.0,
+)
+    isinf(r) && return _distance_Inf(M, p, q, m)
+    (r == 1) && return _distance_1(M, p, q, m)
+    return _distance_Inf(M, p, q, r)
+end
+#
+#
+# The three special cases
+function _distance_r(M::AbstractPowerManifold, p, q, m::AbstractInverseRetractionMethod, r)
+    rep_size = representation_size(M.manifold)
+    values = [
+        distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i), m) for
+        i in get_iterator(M)
+    ]
+    return norm(values, r)
+end
+function _distance_r(M::AbstractPowerManifold, p, q, r)
+    rep_size = representation_size(M.manifold)
+    values = [
+        distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i)) for
+        i in get_iterator(M)
+    ]
+    return norm(values, r)
+end
+function _distance_1(M::AbstractPowerManifold, p, q, m::AbstractInverseRetractionMethod)
+    s = zero(number_eltype(p))
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        s += distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i), m)
+    end
+    return s
+end
+function _distance_1(M::AbstractPowerManifold, p, q)
+    s = zero(number_eltype(p))
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        s += distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i))
+    end
+    return s
+end
+function _distance_Inf(M::AbstractPowerManifold, p, q, m::AbstractInverseRetractionMethod)
+    d = 0.0
+    v = 0.0
     rep_size = representation_size(M.manifold)
     for i in get_iterator(M)
-        sum_squares +=
-            distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i))^2
+        v = distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i), m)
+        d = (v > d) ? v : d
     end
-    return sqrt(sum_squares)
+    return d
+end
+function _distance_Inf(M::AbstractPowerManifold, p, q)
+    d = 0.0
+    v = 0.0
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        v = distance(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, q, i))
+        d = (v > d) ? v : d
+    end
+    return d
 end
 
 @doc raw"""
@@ -881,6 +964,13 @@
     return TangentSpace(M.manifold, p[M, I...])
 end
 
+"""
+    has_components(::AbstractPowerManifold)
+
+Return `true, since points on an [`AbstractPowerManifold`](@ref) consist of components.
+"""
+has_components(::AbstractPowerManifold) = true
+
 @doc raw"""
     injectivity_radius(M::AbstractPowerManifold[, p])
 
@@ -1092,23 +1182,46 @@
 end
 
 @doc raw"""
-    norm(M::AbstractPowerManifold, p, X)
+    norm(M::AbstractPowerManifold, p, X, r::REal=2.0)
 
 Compute the norm of `X` from the tangent space of `p` on an
-[`AbstractPowerManifold`](@ref) `M`, i.e. from the element wise norms the
-Frobenius norm is computed.
+[`AbstractPowerManifold`](@ref) `M`, i.e. from the element wise norms `r`-norm is computed,
+where the default `r=2.0` yields the Frobenius norm is computed.
 """
-function LinearAlgebra.norm(M::AbstractPowerManifold, p, X)
-    sum_squares = zero(number_eltype(X))
+function LinearAlgebra.norm(M::AbstractPowerManifold, p, X, r::Real = 2.0)
+    isinf(r) && return _norm_Inf(M, p, X)
+    (r == 1) && return _norm_1(M, p, X)
+    return _norm_r(M, p, X, r)
+end
+function _norm_r(M::AbstractPowerManifold, p, X, r)
+    rep_size = representation_size(M.manifold)
+    values = [
+        norm(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, X, i)) for
+        i in get_iterator(M)
+    ]
+    return norm(values, r)
+end
+function _norm_1(M::AbstractPowerManifold, p, X)
+    s = zero(number_eltype(p))
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        s += norm(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, X, i))
+    end
+    return s
+end
+function _norm_Inf(M::AbstractPowerManifold, p, X)
+    d = 0.0
+    v = 0.0
     rep_size = representation_size(M.manifold)
     for i in get_iterator(M)
-        sum_squares +=
-            norm(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, X, i))^2
+        v = norm(M.manifold, _read(M, rep_size, p, i), _read(M, rep_size, X, i))
+        d = (v > d) ? v : d
     end
-    return sqrt(sum_squares)
+    return d
 end
 
 
+
 function parallel_transport_direction!(M::AbstractPowerManifold, Y, p, X, d)
     rep_size = representation_size(M.manifold)
     for i in get_iterator(M)

src/ProductManifold.jl

@@ -372,22 +372,49 @@
     return ProductVectorTransport(map(default_vector_transport_method, M.manifolds)...)
 end
 
-@doc raw"""
-    distance(M::ProductManifold, p, q)
-
-Compute the distance between two points `p` and `q` on the [`ProductManifold`](@ref) `M`, which is
-the 2-norm of the elementwise distances on the internal manifolds that build `M`.
-"""
-function distance(M::ProductManifold, p, q)
-    return sqrt(
-        sum(
-            map(
-                distance,
-                M.manifolds,
-                submanifold_components(M, p),
-                submanifold_components(M, q),
-            ) .^ 2,
+_doc_distance_prod = """
+    distance(M::ProductManifold, p, q, r::Real=2.0)
+    distance(M::ProductManifold, p, q, m::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(), r::Real=2.0)
+
+Compute the distance between `q` and `p` on an [`ProductManifold`](@ref).
+
+First, the componentwise distances are computed. These can be approximated using the
+`norm` of an [`AbstractInverseRetractionMethod`](@ref) `m`.
+Then, the `r`-norm of these elements is computed.
+"""
+
+@doc "$(_doc_distance_prod)"
+function distance(M::ProductManifold, p, q, r::Real = 2.0)
+    return norm(
+        map(
+            distance,
+            M.manifolds,
+            submanifold_components(M, p),
+            submanifold_components(M, q),
+        ),
+        r,
+    )
+end
+function distance(M::ProductManifold, p, q, ::LogarithmicInverseRetraction, r::Real = 2.0)
+    return distance(M, p, q, r)
+end
+
+@doc "$(_doc_distance_prod)"
+function distance(
+    M::ProductManifold,
+    p,
+    q,
+    m::AbstractInverseRetractionMethod,
+    r::Real = 2.0,
+)
+    return norm(
+        map(
+            (M, p, q) -> distance(M, p, q, m),
+            M.manifolds,
+            submanifold_components(M, p),
+            submanifold_components(M, q),
         ),
+        r,
     )
 end
 
@@ -557,6 +584,13 @@
     return X
 end
 
+"""
+    has_components(::ProductManifold)
+
+Return `true` since points on an [`ProductManifold`](@ref) consist of components.
+"""
+has_components(::ProductManifold) = true
+
 @doc raw"""
     injectivity_radius(M::ProductManifold)
     injectivity_radius(M::ProductManifold, x)
@@ -725,21 +759,15 @@
 end
 
 @doc raw"""
-    norm(M::ProductManifold, p, X)
+    norm(M::ProductManifold, p, X, r::Real=2.0)
 
-Compute the norm of `X` from the tangent space of `p` on the [`ProductManifold`](@ref),
+Compute the (`r`-)norm of `X` from the tangent space of `p` on the [`ProductManifold`](@ref),
 i.e. from the element wise norms the 2-norm is computed.
 """
-function LinearAlgebra.norm(M::ProductManifold, p, X)
-    norms_squared = (
-        map(
-            norm,
-            M.manifolds,
-            submanifold_components(M, p),
-            submanifold_components(M, X),
-        ) .^ 2
-    )
-    return sqrt(sum(norms_squared))
+function LinearAlgebra.norm(M::ProductManifold, p, X, r::Real = 2.0)
+    norms =
+        (map(norm, M.manifolds, submanifold_components(M, p), submanifold_components(M, X)))
+    return norm(norms, r)
 end
 
 """

test/default_manifold.jl

@@ -31,6 +31,7 @@ using ManifoldsBaseTestUtils
     @test isa(manifold_dimension(M), Integer)
     @test manifold_dimension(M) ≥ 0
     @test base_manifold(M) == M
+    @test has_components(M)
     @test number_system(M) == ManifoldsBase.ℝ
     @test ManifoldsBase.representation_size(M) == (3,)
 

test/empty_manifold.jl

@@ -11,6 +11,7 @@ using Test
     @test base_manifold(M) === M
     @test number_system(M) === ℝ
     @test representation_size(M) === nothing
+    @test !has_components(M)
 
     @test_throws MethodError manifold_dimension(M)