Merge pull request #189 from r-dornig/SphereSymmetricMatrices

Introduces the manifold of unit-norm symmetric matrices

Project.toml

@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>", "Antoine Levitt <[email protected]>"]
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

docs/make.jl

@@ -41,6 +41,8 @@ makedocs(
                 "Symmetric positive definite" =>
                     "manifolds/symmetricpositivedefinite.md",
                 "Torus" => "manifolds/torus.md",
+                "Unit-norm symmetric matrices" =>
+                    "manifolds/spheresymmetricmatrices.md",
             ],
             "Combined manifolds" => [
                 "Graph manifold" => "manifolds/graph.md",

docs/src/manifolds/spheresymmetricmatrices.md

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

src/Manifolds.jl

@@ -129,6 +129,7 @@ include("manifolds/Rotations.jl")
 include("manifolds/SkewSymmetric.jl")
 include("manifolds/Stiefel.jl")
 include("manifolds/Sphere.jl")
+include("manifolds/SphereSymmetricMatrices.jl")
 include("manifolds/Symmetric.jl")
 include("manifolds/SymmetricPositiveDefinite.jl")
 include("manifolds/SymmetricPositiveDefiniteLinearAffine.jl")
@@ -187,6 +188,7 @@ export Euclidean,
     Rotations,
     SkewSymmetricMatrices,
     Sphere,
+    SphereSymmetricMatrices,
     Stiefel,
     SymmetricMatrices,
     SymmetricPositiveDefinite,

src/manifolds/Sphere.jl

@@ -74,7 +74,7 @@ several functions like the [`inner`](@ref inner(::Euclidean, ::Any...)) product
 
 # Constructor
 
-    ArraySphere(n₁,n₂,...,nᵢ; field=ℝ))
+    ArraySphere(n₁,n₂,...,nᵢ; field=ℝ)
 
 Generate sphere in $𝔽^{n_1, n_2, …, n_i}$, where 𝔽 defaults to the real-valued case ℝ.
 """

src/manifolds/SphereSymmetricMatrices.jl

@@ -0,0 +1,147 @@
+@doc raw"""
+    SphereSymmetricMatrices{n,𝔽} <: AbstractEmbeddedManifold{ℝ,TransparentIsometricEmbedding}
+
+The [`Manifold`](@ref) consisting of the $n × n$ symmetric matrices 
+of unit Frobenius norm, i.e. 
+````math
+\mathcal{S}_{\text{sym}} :=\bigl\{p  ∈ 𝔽^{n × n}\ \big|\ p^{\mathrm{H}} = p, \lVert p \rVert = 1 \bigr\},
+````
+where $\cdot^{\mathrm{H}}$ denotes the Hermitian, i.e. complex conjugate transpose,
+and the field $𝔽 ∈ \{ ℝ, ℂ\}$.
+
+# Constructor
+    SphereSymmetricMatrices(n[, field=ℝ])
+
+Generate the manifold of `n`-by-`n` symmetric matrices of unit Frobenius norm.
+"""
+struct SphereSymmetricMatrices{N,𝔽} <:
+       AbstractEmbeddedManifold{𝔽,TransparentIsometricEmbedding} end
+
+function SphereSymmetricMatrices(n::Int, field::AbstractNumbers = ℝ)
+    return SphereSymmetricMatrices{n,field}()
+end
+
+@doc raw"""
+    check_manifold_point(M::SphereSymmetricMatrices{n,𝔽}, p; kwargs...) 
+
+Check whether the matrix is a valid point on the [`SphereSymmetricMatrices`](@ref) `M`, 
+i.e. is an `n`-by-`n` symmetric matrix of unit Frobenius norm.
+
+The tolerance for the symmetry of `p` can be set using `kwargs...`.
+"""
+function check_manifold_point(M::SphereSymmetricMatrices{n,𝔽}, p; kwargs...) where {n,𝔽}
+    mpv =
+        invoke(check_manifold_point, Tuple{supertype(typeof(M)),typeof(p)}, M, p; kwargs...)
+    mpv === nothing || return mpv
+    if !isapprox(norm(p - p'), 0.0; kwargs...)
+        return DomainError(
+            norm(p - p'),
+            "The point $(p) does not lie on $M, since it is not symmetric.",
+        )
+    end
+    return nothing
+end
+
+
+"""
+    check_tangent_vector(M::SphereSymmetricMatrices{n,𝔽}, p, X; check_base_point = true, kwargs... )
+
+Check whether `X` is a tangent vector to manifold point `p` on the
+[`SphereSymmetricMatrices`](@ref) `M`, i.e. `X` has to be a symmetric matrix of size `(n,n)`
+of unit Frobenius norm.
+The optional parameter `check_base_point` indicates, whether to call
+ [`check_manifold_point`](@ref)  for `p`.
+
+The tolerance for the symmetry of `p` and `X` can be set using `kwargs...`.
+"""
+function check_tangent_vector(
+    M::SphereSymmetricMatrices{n,𝔽},
+    p,
+    X;
+    check_base_point = true,
+    kwargs...,
+) where {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, # already checked above
+        kwargs...,
+    )
+    mpv === nothing || return mpv
+    if !isapprox(norm(X - X'), 0.0; kwargs...)
+        return DomainError(
+            norm(X - X'),
+            "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not symmetric.",
+        )
+    end
+    return nothing
+end
+
+function decorated_manifold(M::SphereSymmetricMatrices{n,𝔽}) where {n,𝔽}
+    return ArraySphere(n, n; field = 𝔽)
+end
+
+embed!(M::SphereSymmetricMatrices, q, p) = copyto!(q, p)
+embed!(M::SphereSymmetricMatrices, Y, p, X) = copyto!(Y, X)
+
+@doc raw"""
+    manifold_dimension(M::SphereSymmetricMatrices{n,𝔽})
+
+Return the manifold dimension of the [`SphereSymmetricMatrices`](@ref) `n`-by-`n` symmetric matrix `M` of unit
+Frobenius norm over the number system `𝔽`, i.e.
+
+````math
+\begin{aligned}
+\dim(\mathcal{S}_{\text{sym}})(n,ℝ) &= \frac{n(n+1)}{2} - 1,\\
+\dim(\mathcal{S}_{\text{sym}})(n,ℂ) &= 2\frac{n(n+1)}{2} - n -1.
+\end{aligned}
+````
+"""
+function manifold_dimension(::SphereSymmetricMatrices{n,𝔽}) where {n,𝔽}
+    return div(n * (n + 1), 2) * real_dimension(𝔽) - (𝔽 === ℂ ? n : 0) - 1
+end
+
+@doc raw"""
+    project(M::SphereSymmetricMatrices, p)
+
+Projects `p` from the embedding onto the [`SphereSymmetricMatrices`](@ref) `M`, i.e.
+
+````math
+\operatorname{proj}_{\mathcal{S}_{\text{sym}}}(p) = \frac{1}{2} \bigl( p + p^{\mathrm{H}} \bigr),
+````
+where $\cdot^{\mathrm{H}}$ denotes the Hermitian, i.e. complex conjugate transposed.
+"""
+project(::SphereSymmetricMatrices, ::Any)
+
+function project!(M::SphereSymmetricMatrices, q, p)
+    return project!(get_embedding(M), q, (p + p') ./ 2)
+end
+
+@doc raw"""
+    project(M::SphereSymmetricMatrices, p, X)
+
+Project the matrix `X` onto the tangent space at `p` on the [`SphereSymmetricMatrices`](@ref) `M`, i.e.
+
+````math
+\operatorname{proj}_p(X) = \frac{X + X^{\mathrm{H}}}{2} - ⟨p, \frac{X + X^{\mathrm{H}}}{2}⟩p,
+````
+where $\cdot^{\mathrm{H}}$ denotes the Hermitian, i.e. complex conjugate transposed.
+"""
+project(::SphereSymmetricMatrices, ::Any, ::Any)
+
+function project!(M::SphereSymmetricMatrices, Y, p, X)
+    return project!(get_embedding(M), Y, p, (X .+ X') ./ 2)
+end
+
+@generated representation_size(::SphereSymmetricMatrices{n,𝔽}) where {n,𝔽} = (n, n)
+
+function Base.show(io::IO, ::SphereSymmetricMatrices{n,𝔽}) where {n,𝔽}
+    return print(io, "SphereSymmetricMatrices($(n), $(𝔽))")
+end

test/grassmann.jl

@@ -173,7 +173,7 @@ include("utils.jl")
                     QRInverseRetraction(),
                 ],
                 exp_log_atol_multiplier = 10.0,
-                is_tangent_atol_multiplier = 10.0,
+                is_tangent_atol_multiplier = 20.0,
             )
 
             @testset "inner/norm" begin

test/runtests.jl

@@ -42,6 +42,7 @@ include("probability_simplex.jl")
 include("rotations.jl")
 include("skewsymmetric.jl")
 include("sphere.jl")
+include("sphere_symmetric_matrices.jl")
 include("stiefel.jl")
 include("symmetric.jl")
 include("symmetric_positive_definite.jl")

test/sphere_symmetric_matrices.jl

@@ -0,0 +1,78 @@
+include("utils.jl")
+
+@testset "SphereSymmetricMatrices" begin
+    M = SphereSymmetricMatrices(3)
+    M_complex = SphereSymmetricMatrices(3, ℂ)
+    A = [1 2 3; 2 4 -5; 3 -5 6] / norm([1 2 3; 2 4 -5; 3 -5 6])
+    B = [1 2; 2 4] / norm([1 2; 2 4])                                     #wrong dimensions
+    C = [-3 -im 4; im 5 im; 4 -im 0] / norm([-3 -im 4; im 5 im; 4 -im 0]) #complex
+    D = [1 2 3; 4 5 6; 7 8 9] / norm([1 2 3; 4 5 6; 7 8 9])               #not symmetric
+    E = [1 2 3; 2 4 -5; 3 -5 6]                                           #not of unit norm
+    J = [1.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
+    K = [0.0 1.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
+    @testset "Real Sphere Symmetric Matrices Basics" begin
+        @test repr(M) == "SphereSymmetricMatrices(3, ℝ)"
+        @test representation_size(M) == (3, 3)
+        @test base_manifold(M) === M
+        @test typeof(get_embedding(M)) === ArraySphere{Tuple{3,3},ℝ}
+        @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_throws DomainError is_manifold_point(M, E, true)
+        @test check_tangent_vector(M, A, zeros(3, 3)) === nothing
+        @test_throws DomainError is_tangent_vector(M, A, B, true)
+        @test_throws DomainError is_tangent_vector(M, A, C, true)
+        @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, E, true)
+        @test_throws DomainError is_tangent_vector(M, J, K, true)
+        @test manifold_dimension(M) == 5
+        A2 = similar(A)
+        @test A == project!(M, A2, A)
+        @test A == project(M, A)
+        embed!(M, A2, A)
+        A3 = embed(M, A)
+        @test A2 == A
+        @test A3 == A
+        F = [2 4 -1; 4 9 7; -1 7 5] / norm([2 4 -1; 4 9 7; -1 7 5])
+        G = [-10 9 -8; 9 7 6; -8 6 5] / norm([-10 9 -8; 9 7 6; -8 6 5])
+        test_manifold(
+            M,
+            [A, F, G],
+            test_injectivity_radius = false,
+            test_forward_diff = false,
+            test_reverse_diff = false,
+            test_vector_spaces = true,
+            test_project_tangent = true,
+            test_musical_isomorphisms = true,
+            test_vector_transport = true,
+            is_tangent_atol_multiplier = 2,
+        )
+    end
+    @testset "Complex Sphere Symmetric Matrices Basics" begin
+        @test repr(M_complex) == "SphereSymmetricMatrices(3, ℂ)"
+        @test manifold_dimension(M_complex) == 8
+        H =
+            [7.0 -1.5im 0.0; 1.5im 3.0 -1.0im; 0.0 1.0im 4.5] /
+            norm([7.0 -1.5im 0.0; 1.5im 3.0 -1.0im; 0.0 1.0im 4.5])
+        I =
+            [2.0 4.0 5.0im; 4.0 -1.0 -9.0; -5.0im -9.0 2.5] /
+            norm([2.0 4.0 5.0im; 4.0 -1.0 -9.0; -5.0im -9.0 2.5])
+        test_manifold(
+            M_complex,
+            [C, H, I],
+            test_injectivity_radius = false,
+            test_forward_diff = false,
+            test_reverse_diff = false,
+            test_vector_spaces = true,
+            test_project_tangent = true,
+            test_musical_isomorphisms = true,
+            test_vector_transport = true,
+            is_tangent_atol_multiplier = 2,
+            is_point_atol_multiplier = 2,
+            projection_atol_multiplier = 2,
+            exp_log_atol_multiplier = 2,
+        )
+    end
+end