[WIP] HermitianPositiveDefinite

NEWS.md

@@ -5,6 +5,14 @@ All notable changes to ´Manifolds.jl´ 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.11.24] unreleased
+
+### Added
+
+* `HermitianPositiveDefinite` manifold of Hermitian positive definite matrices.
+* `MatrixSqrtManifoldPoint` as a generalization/replacement for the `SPDPoint`
+* `SymmetricPositiveDefinite` is now a `const` of `HermitianPositiveDefinite` over the reals
+
 ## [0.11.21] 2025-04-03
 
 ### Changed

src/Manifolds.jl

@@ -469,6 +469,7 @@ include("manifolds/FlagStiefel.jl")
 include("manifolds/GeneralizedGrassmann.jl")
 include("manifolds/GeneralizedStiefel.jl")
 include("manifolds/HeisenbergMatrices.jl")
+include("manifolds/HermitianPositiveDefinite.jl")
 include("manifolds/Hyperbolic.jl")
 include("manifolds/Hyperrectangle.jl")
 include("manifolds/InvertibleMatrices.jl")
@@ -626,6 +627,7 @@ export Euclidean,
     Grassmann,
     HamiltonianMatrices,
     HeisenbergMatrices,
+    HermitianPositiveDefinite,
     Hyperbolic,
     Hyperrectangle,
     InvertibleMatrices,

src/manifolds/HermitianPositiveDefinite.jl

@@ -0,0 +1,206 @@
+@doc raw"""
+    HermitianPositiveDefinite{𝔽,T} <: AbstractDecoratorManifold{𝔽}
+
+The manifold of hermitian positive definite matrices
+
+````math
+\mathcal H(n) :=
+\bigl\{
+p ∈ 𝔽^{n×n}\ \big|\ a^\mathrm{T}pa > 0 \text{ for all } a ∈ 𝔽^{n}\backslash\{0\}
+\bigr\},
+````
+where usually ``𝔽=ℂ``.
+For the case ``𝔽=ℝ`` this manifold simplifies to the [`SymmetricPositiveDefinite`](@ref)
+
+The tangent space at ``p∈\mathcal H(n)`` reads
+
+```math
+    T_p\mathcal H(n) =
+    \bigl\{
+        X \in 𝔽^{n×n} \big|\ X=X^\mathrm{H}
+    \bigr\},
+```
+i.e. the set of hermitian matrices.
+
+# Constructor
+
+    HermitianPositiveDefinite(n, 𝔽=ℂ; parameter::Symbol=:type)
+
+generates the manifold of hermitian positive definite matrices ``\mathcal H(n) \subset 𝔽^{n×n}``.
+"""
+struct HermitianPositiveDefinite{𝔽, T} <: AbstractDecoratorManifold{𝔽}
+    size::T
+end
+function HermitianPositiveDefinite(n, 𝔽::AbstractNumbers = ℂ; parameter::Symbol = :type)
+    size = wrap_type_parameter(parameter, (n,))
+    return HermitianPositiveDefinite{𝔽, typeof(size)}(size)
+end
+
+@doc raw"""
+    MatrixSqrtManifoldPoint{A,P,Q,R,E} <: AbstractManifoldsPoint
+
+A point on a manifold, that is represented by a matrix ``p ∈ 𝔽^{n×n}`` over a field ``𝔽 ∈ \{ℂ,ℝ\}``,
+which can cache the computation of its Eigen values and Eigen vectors as well as
+its matrix square root and its inverse square root.
+
+This is for example the case for the [`HermitianPositiveDefinite`](@ref) manifold.
+
+# Fields
+
+* `p::P`` a matrix representation of the point
+* `eigen::E` the eigenvalues of `p`
+* `sqrt::Q` the matrix square root of `p`
+* `sqrt_inv::R` the inverse of the matrix square root of `p`
+
+Any of the fields `p`, `sqrt`, `sqrt_inv` can be set to not be stored and is then `missing` (its field type hence `Missing`)
+to indicate that that field should not be stored/cached. If given, the field types `P`, `Q`, and `R` have to be the same as the type `A`.
+The result of `eigen(p)` will always be stored. This way the other three can always be provided from this field,
+but it might be beneficial to cache them as well.
+
+# Constructor
+
+    MatrixSqrtManifoldPoint(
+        p::AbstractMatrix; store_p=true, store_sqrt=true, store_sqrt_inv=true
+    )
+
+Create an `MatrixSqrtManifoldPoint` point using an matrix `p`, where you can optionally store `p`, `sqrt` and `sqrt_inv`
+"""
+struct MatrixSqrtManifoldPoint{
+        A <: AbstractMatrix, P <: Union{A, Missing}, Q <: Union{A, Missing}, R <: Union{A, Missing},
+        E <: Eigen,
+    } <: AbstractManifoldPoint
+    p::P
+    eigen::E
+    sqrt::Q
+    sqrt_inv::R
+end
+
+MatrixSqrtManifoldPoint(p::MatrixSqrtManifoldPoint) = p
+function MatrixSqrtManifoldPoint(
+        p::A;
+        store_p = true, store_sqrt = true, store_sqrt_inv = true,
+    ) where {A}
+    e = eigen(Symmetric(p))
+    U = e.vectors
+    S = max.(e.values, floatmin(eltype(e.values)))
+    if store_sqrt
+        s_sqrt = Diagonal(sqrt.(S))
+        p_sqrt = U * s_sqrt * transpose(U)
+    else
+        p_sqrt = missing
+    end
+    if store_sqrt_inv
+        s_sqrt_inv = Diagonal(1 ./ sqrt.(S))
+        p_sqrt_inv = U * s_sqrt_inv * transpose(U)
+    else
+        p_sqrt_inv = missing
+    end
+    if store_p
+        q = p
+    else
+        q = missing
+    end
+    return MatrixSqrtManifoldPoint{A, typeof(q), typeof(p_sqrt), typeof(p_sqrt_inv), typeof(e)}(
+        q, e, p_sqrt, p_sqrt_inv,
+    )
+end
+convert(::Type{MatrixSqrtManifoldPoint}, p::AbstractMatrix) = MatrixSqrtManifoldPoint(p)
+
+function Base.:(==)(p::MatrixSqrtManifoldPoint, q::MatrixSqrtManifoldPoint)
+    return p.eigen == q.eigen
+end
+
+function allocate(p::MatrixSqrtManifoldPoint)
+    return MatrixSqrtManifoldPoint(
+        ismissing(p.p) ? missing : allocate(p.p),
+        Eigen(allocate(p.eigen.values), allocate(p.eigen.vectors)),
+        ismissing(p.sqrt) ? missing : allocate(p.sqrt),
+        ismissing(p.sqrt_inv) ? missing : allocate(p.sqrt_inv),
+    )
+end
+function allocate(p::MatrixSqrtManifoldPoint, ::Type{T}) where {T}
+    return MatrixSqrtManifoldPoint(
+        ismissing(p.p) ? missing : allocate(p.p, T),
+        Eigen(allocate(p.eigen.values, T), allocate(p.eigen.vectors, T)),
+        ismissing(p.sqrt) ? missing : allocate(p.sqrt, T),
+        ismissing(p.sqrt_inv) ? missing : allocate(p.sqrt_inv, T),
+    )
+end
+
+function allocate_result(
+        M::HermitianPositiveDefinite, ::typeof(zero_vector), p::MatrixSqrtManifoldPoint,
+    )
+    return allocate_result(M, zero_vector, convert(AbstractMatrix, p))
+end
+function allocate_coordinates(
+        M::HermitianPositiveDefinite, p::MatrixSqrtManifoldPoint, T, n::Int,
+    )
+    return allocate_coordinates(M, convert(AbstractMatrix, p), T, n)
+end
+function allocate_result(
+        M::HermitianPositiveDefinite, ::typeof(get_vector), p::MatrixSqrtManifoldPoint, c,
+    )
+    return allocate_result(M, get_vector, convert(AbstractMatrix, p), c)
+end
+
+@doc raw"""
+    check_point(M::HermitianPositiveDefinite, p; kwargs...)
+
+checks, whether `p` is a valid point on the [`HermitianPositiveDefinite`](@ref) `M`,
+i.e. is a matrix ``p ∈ 𝔽^{n×n}`` over a field ``𝔽 ∈ \{ℂ,ℝ\}`` and is hermitian (``p^{\mathrm{H}} = p``)
+and positive definite, that is a^\mathrm{T}pa > 0$ for all $a ∈ 𝔽^{n}\backslash\{0\}$.
+The tolerance for the second to last test can be set using the `kwargs...`.
+"""
+function check_point(M::HermitianPositiveDefinite, p; kwargs...)
+    if !isapprox(p, p'; kwargs...)
+        return DomainError(
+            norm(p - p'), "The point $(p) does not lie on $(M) since its not a hermitian matrix.",
+        )
+    end
+    if !isposdef(p)
+        return DomainError(
+            eigvals(p), "The point $p does not lie on $(M) since its not a positive definite matrix.",
+        )
+    end
+    return nothing
+end
+
+"""
+    check_vector(M::HermitianPositiveDefinite, p, X; kwargs... )
+
+Check whether `X` is a tangent vector to `p` on the [`HermitianPositiveDefinite`](@ref) `M`,
+i.e. a symmetric matrix.
+The tolerance for the last test can be set using the `kwargs...`.
+"""
+function check_vector(M::HermitianPositiveDefinite, p, X; kwargs...)
+    if !isapprox(X, X'; kwargs...)
+        return DomainError(
+            X,
+            "The vector $(X) is not a tangent to a point on $(M) (represented as an element of the Lie algebra) since its not hermitian.",
+        )
+    end
+    return nothing
+end
+
+# Internal function for nicer printing.
+get_parameter_type(::HermitianPositiveDefinite{𝔽, <:TypeParameter}) where {𝔽} = :type
+get_parameter_type(::HermitianPositiveDefinite{𝔽, Tuple{Int}}) where {𝔽} = :field
+
+@doc raw"""
+    representation_size(M::HermitianPositiveDefinite)
+
+Return the size of an array representing an element on the
+[`HermitianPositiveDefinite`](@ref) manifold `M`, i.e. ``n×n``, the size of such a
+hermitian positive definite matrix on ``\mathcal M = \mathcal H(n)``.
+"""
+function representation_size(M::HermitianPositiveDefinite)
+    n = get_parameter(M.size)[1]
+    return (n, n)
+end
+
+function Base.show(io::IO, M::HermitianPositiveDefinite{𝔽}) where {𝔽}
+    n = get_parameter(M.size)[1]
+    p = 𝔽 === ℂ ? "" : ", $𝔽"
+    kw = get_parameter_type(M) === :type ? "" : "; parameter=:$(get_parameter_type(M))"
+    return print(io, "HermitianPositiveDefinite($n$p$kw)")
+end

src/manifolds/SymmetricPositiveDefinite.jl

@@ -4,19 +4,16 @@
 The manifold of symmetric positive definite matrices, i.e.
 
 ````math
-\mathcal P(n) =
-\bigl\{
-p ∈ ℝ^{n×n}\ \big|\ a^\mathrm{T}pa > 0 \text{ for all } a ∈ ℝ^{n}\backslash\{0\}
+\mathcal P(n) = \bigl\{
+  p ∈ ℝ^{n×n}\ \big|\ a^\mathrm{T}pa > 0 \text{ for all } a ∈ ℝ^{n}\backslash\{0\}
 \bigr\}
 ````
 
 The tangent space at ``T_p\mathcal P(n)`` reads
 
 ```math
     T_p\mathcal P(n) =
-    \bigl\{
-        X \in \mathbb R^{n×n} \big|\ X=X^\mathrm{T}
-    \bigr\},
+    \bigl\{ X \in \mathbb R^{n×n} \big|\ X=X^\mathrm{T} \bigr\},
 ```
 i.e. the set of symmetric matrices,
 
@@ -26,9 +23,7 @@ i.e. the set of symmetric matrices,
 
 generates the manifold ``\mathcal P(n) \subset ℝ^{n×n}``
 """
-struct SymmetricPositiveDefinite{T} <: AbstractDecoratorManifold{ℝ}
-    size::T
-end
+const SymmetricPositiveDefinite{T} = HermitianPositiveDefinite{ℝ, T}
 
 function SymmetricPositiveDefinite(n::Int; parameter::Symbol = :type)
     size = wrap_type_parameter(parameter, (n,))
@@ -135,7 +130,7 @@ function check_point(M::SymmetricPositiveDefinite, p; kwargs...)
     if !isapprox(norm(p - transpose(p)), 0.0; kwargs...)
         return DomainError(
             norm(p - transpose(p)),
-            "The point $(p) does not lie on $(M) since its not a symmetric matrix:",
+            "The point $(p) does not lie on $(M) since its not a symmetric matrix.",
         )
     end
     if !isposdef(p)
@@ -481,12 +476,10 @@ function representation_size(M::SymmetricPositiveDefinite)
     return (N, N)
 end
 
-function Base.show(io::IO, ::SymmetricPositiveDefinite{TypeParameter{Tuple{n}}}) where {n}
-    return print(io, "SymmetricPositiveDefinite($(n))")
-end
-function Base.show(io::IO, M::SymmetricPositiveDefinite{Tuple{Int}})
+function Base.show(io::IO, M::SymmetricPositiveDefinite)
     n = get_parameter(M.size)[1]
-    return print(io, "SymmetricPositiveDefinite($(n); parameter=:field)")
+    kw = get_parameter_type(M) === :type ? "" : "; parameter=:$(get_parameter_type(M))"
+    return print(io, "SymmetricPositiveDefinite($n$kw)")
 end
 
 function Base.show(io::IO, ::MIME"text/plain", p::SPDPoint)

test/manifolds/test_hermitian_symmetric_positive_definite.jl

@@ -0,0 +1,10 @@
+using Manifolds, Test
+
+@testset "HermitianPositiveDefinite" begin
+    M = HermitianPositiveDefinite(2)
+    M2 = HermitianPositiveDefinite(2, ℂ)
+    M3 = HermitianPositiveDefinite(2; parameter = :field)
+
+    A = [4.0 -2im; 2.0im 4.0]
+    @test is_point(M, A, true)
+end

test/runtests.jl

@@ -41,6 +41,8 @@ end
         include_test("manifolds/test_heisenberg_matrices.jl")
         include_test("manifolds/test_stiefel.jl")
         include_test("manifolds/test_unitary_matrices.jl")
+        include_test("manifolds/test_hermitian_symmetric_positive_definite.jl")
+
     end
     (TEST_SET ∈ ["all", "manifolds"]) && Test.@testset "Manifolds.jl Old Tests" begin
         # starting with tests of simple manifolds