Introduces the manifold of centered matrices

docs/make.jl

@@ -19,6 +19,7 @@ makedocs(
         "ManifoldsBase.jl" => "interface.md",
         "Manifolds" => [
             "Basic manifolds" => [
+                "Centered matrices" => "manifolds/centeredmatrices.md",
                 "Cholesky space" => "manifolds/choleskyspace.md",
                 "Circle" => "manifolds/circle.md",
                 "Euclidean" => "manifolds/euclidean.md",

docs/src/manifolds/centeredmatrices.md

@@ -0,0 +1,7 @@
+# Centered matrices
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["manifolds/CenteredMatrices.jl"]
+Order = [:type, :function]
+```

src/Manifolds.jl

@@ -114,6 +114,7 @@ include("manifolds/GraphManifold.jl")
 include("manifolds/Euclidean.jl")
 include("manifolds/Lorentz.jl")
 
+include("manifolds/CenteredMatrices.jl")
 include("manifolds/CholeskySpace.jl")
 include("manifolds/Circle.jl")
 include("manifolds/FixedRankMatrices.jl")
@@ -169,6 +170,7 @@ export CoTVector, Manifold, MPoint, TVector, Manifold
 export AbstractSphere
 export Euclidean,
     ArraySphere,
+    CenteredMatrices,
     CholeskySpace,
     Circle,
     FixedRankMatrices,

src/manifolds/CenteredMatrices.jl

@@ -0,0 +1,144 @@
+@doc raw"""
+    CenteredMatrices{m,n,𝔽} <: AbstractEmbeddedManifold{𝔽,TransparentIsometricEmbedding}
+
+The manifold of $m × n$ real-valued or complex-valued matrices whose columns sum to zero, i.e.
+````math
+\bigl\{ p ∈ 𝔽^{m × n}\ \big|\ [1 … 1] * p = [0 … 0] \bigr\},
+````
+where $𝔽 ∈ \{ℝ,ℂ\}$.
+
+# Constructor
+    CenteredMatrices(m, n[, field=ℝ])
+
+Generate the manifold of `m`-by-`n` (`field`-valued) matrices whose columns sum to zero.
+"""
+struct CenteredMatrices{M,N,𝔽} <: AbstractEmbeddedManifold{𝔽,TransparentIsometricEmbedding} end
+
+function CenteredMatrices(m::Int, n::Int, field::AbstractNumbers = ℝ)
+    return CenteredMatrices{m,n,field}()
+end
+
+@doc raw"""
+    check_manifold_point(M::CenteredMatrices{m,n,𝔽}, p; kwargs...) 
+
+Check whether the matrix is a valid point on the
+[`CenteredMatrices`](@ref) `M`, i.e. is an `m`-by-`n` matrix whose columns sum to 
+zero.
+
+The tolerance for the column sums of `p` can be set using `kwargs...`.
+"""
+function check_manifold_point(M::CenteredMatrices{m,n,𝔽}, p; kwargs...) where {m,n,𝔽}
+    mpv =
+        invoke(check_manifold_point, Tuple{supertype(typeof(M)),typeof(p)}, M, p; kwargs...)
+    mpv === nothing || return mpv
+    if !isapprox(sum(p, dims = 1), zeros(1, n); kwargs...)
+        return DomainError(
+            p,
+            string("The point $(p) does not lie on $(M), since its columns do not sum to zero."),
+        )
+    end
+    return nothing
+end
+
+
+"""
+    check_tangent_vector(M::CenteredMatrices{m,n,𝔽}, p, X; check_base_point = true, kwargs... )
+
+Check whether `X` is a tangent vector to manifold point `p` on the
+[`CenteredMatrices`](@ref) `M`, i.e. that `X` is a matrix of size `(m,n)` whose columns
+sum to zero and its values are from the correct [`AbstractNumbers`](@ref).
+The optional parameter `check_base_point` indicates, whether to call
+ [`check_manifold_point`](@ref)  for `p`.
+The tolerance for the column sums of `p` and `X` can be set using `kwargs...`.
+"""
+function check_tangent_vector(
+    M::CenteredMatrices{m,n,𝔽},
+    p,
+    X;
+    check_base_point = true,
+    kwargs...,
+) where {m,n,𝔽}
+    if check_base_point
+        mpe = check_manifold_point(M, p; kwargs...)
+        mpe === nothing || return mpe
+    end
+    mpv = invoke(
+        check_tangent_vector,
+        Tuple{supertype(typeof(M)),typeof(p),typeof(X)},
+        M,
+        p,
+        X;
+        check_base_point = false,
+        kwargs...,
+    )
+    mpv === nothing || return mpv
+    if !isapprox(sum(X, dims = 1), zeros(1, n); kwargs...)
+        return DomainError(
+            X,
+            "The vector $(X) is not a tangent vector to $(p) on $(M), since its columns do not sum to zero.",
+        )
+    end
+    return nothing
+end
+
+decorated_manifold(M::CenteredMatrices{m,n,𝔽}) where {m,n,𝔽} = Euclidean(m, n; field = 𝔽)
+
+embed!(M::CenteredMatrices, q, p) = copyto!(q, p)
+embed!(M::CenteredMatrices, Y, p, X) = copyto!(Y, X)
+
+@doc raw"""
+    manifold_dimension(M::CenteredMatrices{m,n,𝔽})
+
+Return the manifold dimension of the [`CenteredMatrices`](@ref) `m`-by-`n` matrix `M` over the number system
+`𝔽`, i.e.
+
+````math
+\dim(\mathcal M) = (m*n - n) \dim_ℝ 𝔽,
+````
+where $\dim_ℝ 𝔽$ is the [`real_dimension`](@ref) of `𝔽`.
+"""
+function manifold_dimension(::CenteredMatrices{m,n,𝔽}) where {m,n,𝔽}
+    return (m * n - n) * real_dimension(𝔽)
+end
+
+@doc raw"""
+    project(M::CenteredMatrices, p)
+
+Projects `p` from the embedding onto the [`CenteredMatrices`](@ref) `M`, i.e.
+
+````math
+\operatorname{proj}_{\mathcal M}(p) = p - \begin{bmatrix} 
+1\\
+⋮\\
+1
+\end{bmatrix} * [c_1 \dots c_n],
+````
+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))
+
+@doc raw"""
+    project(M::CenteredMatrices, p, X)
+
+Project the matrix `X` onto the tangent space at `p` on the [`CenteredMatrices`](@ref) `M`, i.e.
+
+````math
+\operatorname{proj}_p(X) = X - \begin{bmatrix} 
+1\\
+⋮\\
+1
+\end{bmatrix} * [c_1 \dots c_n],
+````
+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))
+
+@generated representation_size(::CenteredMatrices{m,n,𝔽}) where {m,n,𝔽} = (m, n)
+
+function Base.show(io::IO, ::CenteredMatrices{m,n,𝔽}) where {m,n,𝔽}
+    return print(io, "CenteredMatrices($(m), $(n), $(𝔽))")
+end

test/centered_matrices.jl

@@ -0,0 +1,46 @@
+include("utils.jl")
+
+@testset "CenteredMatrices" begin
+    M = CenteredMatrices(3, 2)
+    M_complex = CenteredMatrices(3, 2, ℂ)
+    @test repr(M_complex) == "CenteredMatrices(3, 2, ℂ)"
+    @test manifold_dimension(M_complex) == 8
+    A = [1 2; 4 5; -5 -7]
+    B = [1 2 3; 4 5 6; -5 -7 -9]    #wrong dimensions
+    C = [-3 -im; 2 im; 1 0]         #complex
+    D = [1 2; 3 4; 5 6]             #not centered 
+    @testset "Real Centered Matrices Basics" begin
+        @test repr(M) == "CenteredMatrices(3, 2, ℝ)"
+        @test representation_size(M) == (3, 2)
+        @test base_manifold(M) === M
+        @test typeof(get_embedding(M)) === Euclidean{Tuple{3,2},ℝ}
+        @test check_manifold_point(M, A) === nothing
+        @test_throws DomainError is_manifold_point(M, B, true)
+        @test_throws DomainError is_manifold_point(M, C, true)
+        @test_throws DomainError is_manifold_point(M, D, true)
+        @test check_tangent_vector(M, A, A) === nothing
+        @test_throws DomainError is_tangent_vector(M, A, D, true)
+        @test_throws DomainError is_tangent_vector(M, D, A, true)
+        @test_throws DomainError is_tangent_vector(M, A, B, true)
+        @test manifold_dimension(M) == 4
+        @test A == project!(M, A, A)
+        @test A == project(M, A, A)
+        A2 = similar(A)
+        embed!(M, A2, A)
+        A3 = embed(M, A)
+        @test A2 == A
+        @test A3 == A
+    end
+    types = [Matrix{Float64}]
+    E = [0 0; 1 -1; -1 1]
+    F = [0.5 1; -1 -0.7; 0.5 -0.3]
+    test_manifold(
+        M,
+        [A, E, F],
+        test_injectivity_radius = false,
+        test_reverse_diff = false,
+        test_project_tangent = true,
+        test_musical_isomorphisms = true,
+        test_vector_transport = true,
+    )
+end

test/runtests.jl

@@ -29,6 +29,7 @@ include("groups/group_utils.jl")
 include("sized_abstract_array.jl")
 
 # starting with tests of simple manifolds
+include("centered_matrices.jl")
 include("circle.jl")
 include("cholesky_space.jl")
 include("euclidean.jl")

test/utils.jl

@@ -92,6 +92,7 @@ function test_manifold(
     projection_atol_multiplier = 0,
     rand_tvector_atol_multiplier = 0,
     is_tangent_atol_multiplier = 1,
+    is_point_atol_multiplier = 1,
 )
 
     length(pts) ≥ 3 || error("Not enough points (at least three expected)")
@@ -152,9 +153,9 @@ function test_manifold(
 
     @testset "is_manifold_point" begin
         for pt in pts
-            @test is_manifold_point(M, pt)
-            @test check_manifold_point(M, pt) === nothing
-            @test check_manifold_point(M, pt) === nothing
+            atol = is_point_atol_multiplier * find_eps(pt)
+            @test is_manifold_point(M, pt; atol = atol)
+            @test check_manifold_point(M, pt; atol = atol) === nothing
         end
     end