Bases and rand for InvertibleMatrices and HeisenbergMatrices (#777)

* Bases and rand for InvertibleMatrices and HeisenbergMatrices

* improve coverage

* improve coverage

NEWS.md

@@ -5,6 +5,12 @@ 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.10.11] - 2025-01-02
+
+### Added
+
+* Bases and rand for `HeisenbergMatrices` and `InvertibleMatrices`.
+
 ## [0.10.10] - 2024-12-20
 
 ### Added

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.10.10"
+version = "0.10.11"
 
 [deps]
 Einsum = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0"

src/manifolds/HeisenbergMatrices.jl

@@ -31,6 +31,15 @@ function HeisenbergMatrices(n::Int; parameter::Symbol=:type)
     return HeisenbergMatrices{typeof(size)}(size)
 end
 
+function _heisenberg_a_view(M::HeisenbergMatrices, p)
+    n = get_parameter(M.size)[1]
+    return view(p, 1, 2:(n + 1))
+end
+function _heisenberg_b_view(M::HeisenbergMatrices, p)
+    n = get_parameter(M.size)[1]
+    return view(p, 2:(n + 1), n + 2)
+end
+
 function active_traits(f, ::HeisenbergMatrices, args...)
     return merge_traits(IsEmbeddedSubmanifold())
 end
@@ -97,6 +106,33 @@ end
 embed(::HeisenbergMatrices, p) = p
 embed(::HeisenbergMatrices, p, X) = X
 
+@doc raw"""
+    get_coordinates(M::HeisenbergMatrices, p, X, ::DefaultOrthonormalBasis{ℝ,TangentSpaceType})
+
+Get coordinates of tangent vector `X` at point `p` from the [`HeisenbergMatrices`](@ref) `M`.
+Given a matrix
+```math
+\begin{bmatrix} 1 & \mathbf{a} & c \\
+\mathbf{0} & I_n & \mathbf{b} \\
+0 & \mathbf{0} & 1 \end{bmatrix}
+```
+the coordinates are concatenated vectors ``\mathbf{a}``, ``\mathbf{b}``, and number ``c``.
+"""
+get_coordinates(::HeisenbergMatrices, p, X, ::DefaultOrthonormalBasis{ℝ,TangentSpaceType})
+
+function get_coordinates_orthonormal(M::HeisenbergMatrices, p, X, ::RealNumbers)
+    n = get_parameter(M.size)[1]
+    return vcat(_heisenberg_a_view(M, X), _heisenberg_b_view(M, X), X[1, n + 2])
+end
+
+function get_coordinates_orthonormal!(M::HeisenbergMatrices, Xⁱ, p, X, ::RealNumbers)
+    n = get_parameter(M.size)[1]
+    Xⁱ[1:n] .= _heisenberg_a_view(M, X)
+    Xⁱ[(n + 1):(2 * n)] .= _heisenberg_b_view(M, X)
+    Xⁱ[2 * n + 1] = X[1, n + 2]
+    return Xⁱ
+end
+
 function get_embedding(::HeisenbergMatrices{TypeParameter{Tuple{n}}}) where {n}
     return Euclidean(n + 2, n + 2)
 end
@@ -105,6 +141,37 @@ function get_embedding(M::HeisenbergMatrices{Tuple{Int}})
     return Euclidean(n + 2, n + 2; parameter=:field)
 end
 
+@doc raw"""
+    get_vector(M::HeisenbergMatrices, p, Xⁱ, ::DefaultOrthonormalBasis{ℝ,TangentSpaceType})
+
+Get tangent vector with coordinates `Xⁱ` at point `p` from the [`HeisenbergMatrices`](@ref) `M`.
+Given a vector of coordinates ``\begin{bmatrix}\mathbb{a} & \mathbb{b} & c\end{bmatrix}`` the tangent vector is equal to
+```math
+\begin{bmatrix} 1 & \mathbf{a} & c \\
+\mathbf{0} & I_n & \mathbf{b} \\
+0 & \mathbf{0} & 1 \end{bmatrix}
+```
+"""
+get_vector(M::HeisenbergMatrices, p, c, ::DefaultOrthonormalBasis{ℝ,TangentSpaceType})
+
+function get_vector_orthonormal(M::HeisenbergMatrices, p, Xⁱ, ::RealNumbers)
+    n = get_parameter(M.size)[1]
+    return [
+        0 Xⁱ[1:n] Xⁱ[2 * n + 1]
+        zeros(n, n + 1) Xⁱ[(n + 1):(2 * n)]'
+        zeros(1, n + 2)
+    ]
+end
+
+function get_vector_orthonormal!(M::HeisenbergMatrices, X, p, Xⁱ, ::RealNumbers)
+    n = get_parameter(M.size)[1]
+    fill!(X, 0)
+    X[1, 2:(n + 1)] .= Xⁱ[1:n]
+    X[2:(n + 1), n + 2] .= Xⁱ[(n + 1):(2 * n)]
+    X[1, n + 2] = Xⁱ[2 * n + 1]
+    return X
+end
+
 """
     is_flat(::HeisenbergMatrices)
 
@@ -119,6 +186,44 @@ Return the dimension of [`HeisenbergMatrices`](@ref)`(n)`, which is equal to ``2
 """
 manifold_dimension(M::HeisenbergMatrices) = 2 * get_parameter(M.size)[1] + 1
 
+@doc raw"""
+    Random.rand(M::HeisenbergMatrices; vector_at = nothing, σ::Real=1.0)
+
+If `vector_at` is `nothing`, return a random point on the [`HeisenbergMatrices`](@ref) `M`
+by sampling elements of the first row and the last column from the normal distribution with
+mean 0 and standard deviation `σ`.
+
+If `vector_at` is not `nothing`, return a random tangent vector from the tangent space of
+the point `vector_at` on the [`HeisenbergMatrices`](@ref) by using a normal distribution with
+mean 0 and standard deviation `σ`.
+"""
+rand(M::HeisenbergMatrices; vector_at=nothing, σ::Real=1.0)
+
+function Random.rand!(
+    rng::AbstractRNG,
+    M::HeisenbergMatrices,
+    pX;
+    σ::Real=one(eltype(pX)),
+    vector_at=nothing,
+)
+    n = ManifoldsBase.get_parameter(M.size)[1]
+    if vector_at === nothing
+        copyto!(pX, I)
+        va = view(pX, 1, 2:(n + 2))
+        randn!(rng, va)
+        va .*= σ
+        vb = view(pX, 2:(n + 1), n + 2)
+        randn!(rng, vb)
+        vb .*= σ
+    else
+        fill!(pX, 0)
+        randn!(rng, view(pX, 1, 2:(n + 2)))
+        randn!(rng, view(pX, 2:(n + 1), n + 2))
+        pX .*= σ
+    end
+    return pX
+end
+
 function Base.show(io::IO, ::HeisenbergMatrices{TypeParameter{Tuple{n}}}) where {n}
     return print(io, "HeisenbergMatrices($(n))")
 end

src/manifolds/InvertibleMatrices.jl

@@ -59,6 +59,25 @@ end
 embed(::InvertibleMatrices, p) = p
 embed(::InvertibleMatrices, p, X) = X
 
+function get_coordinates(
+    ::InvertibleMatrices{ℝ,<:Any},
+    p,
+    X,
+    ::DefaultOrthonormalBasis{ℝ,TangentSpaceType},
+)
+    return vec(X)
+end
+
+function get_coordinates!(
+    ::InvertibleMatrices{ℝ,<:Any},
+    Xⁱ,
+    p,
+    X,
+    ::DefaultOrthonormalBasis{ℝ,TangentSpaceType},
+)
+    return copyto!(Xⁱ, X)
+end
+
 function get_embedding(::InvertibleMatrices{𝔽,TypeParameter{Tuple{n}}}) where {n,𝔽}
     return Euclidean(n, n; field=𝔽)
 end
@@ -67,6 +86,26 @@ function get_embedding(M::InvertibleMatrices{𝔽,Tuple{Int}}) where {𝔽}
     return Euclidean(n, n; field=𝔽, parameter=:field)
 end
 
+function get_vector(
+    M::InvertibleMatrices{ℝ,<:Any},
+    p,
+    Xⁱ,
+    ::DefaultOrthonormalBasis{ℝ,TangentSpaceType},
+)
+    n = get_parameter(M.size)[1]
+    return reshape(Xⁱ, n, n)
+end
+
+function get_vector!(
+    ::InvertibleMatrices{ℝ,<:Any},
+    X,
+    p,
+    Xⁱ,
+    ::DefaultOrthonormalBasis{ℝ,TangentSpaceType},
+)
+    return copyto!(X, Xⁱ)
+end
+
 """
     is_flat(::InvertibleMatrices)
 
@@ -84,6 +123,27 @@ function manifold_dimension(M::InvertibleMatrices{<:Any,𝔽}) where {𝔽}
     return manifold_dimension(get_embedding(M))
 end
 
+@doc raw"""
+    Random.rand(M::InvertibleMatrices; vector_at=nothing, kwargs...)
+
+If `vector_at` is `nothing`, return a random point on the [`InvertibleMatrices`](@ref)
+manifold `M` by using `rand` in the embedding.
+
+If `vector_at` is not `nothing`, return a random tangent vector from the tangent space of
+the point `vector_at` on the [`InvertibleMatrices`](@ref) by using by using `rand` in the
+embedding.
+"""
+rand(M::InvertibleMatrices; kwargs...)
+
+function Random.rand!(M::InvertibleMatrices, pX; kwargs...)
+    rand!(get_embedding(M), pX; kwargs...)
+    return pX
+end
+function Random.rand!(rng::AbstractRNG, M::InvertibleMatrices, pX; kwargs...)
+    rand!(rng, get_embedding(M), pX; kwargs...)
+    return pX
+end
+
 function Base.show(io::IO, ::InvertibleMatrices{𝔽,TypeParameter{Tuple{n}}}) where {n,𝔽}
     return print(io, "InvertibleMatrices($(n), $(𝔽))")
 end

test/manifolds/heisenberg_matrices.jl

@@ -30,9 +30,12 @@ using LinearAlgebra, Manifolds, ManifoldsBase, Test
     test_manifold(
         M,
         pts;
+        basis_types_to_from=(DefaultOrthonormalBasis(),),
         parallel_transport=true,
         test_injectivity_radius=true,
         test_musical_isomorphisms=false,
+        test_rand_point=true,
+        test_rand_tvector=true,
     )
 
     @test all(

test/manifolds/invertible_matrices.jl

@@ -1,4 +1,4 @@
-using LinearAlgebra, Manifolds, ManifoldsBase, Test
+using LinearAlgebra, Manifolds, ManifoldsBase, Test, Random
 
 @testset "Invertible matrices" begin
     M = InvertibleMatrices(3, ℝ)
@@ -24,6 +24,20 @@ using LinearAlgebra, Manifolds, ManifoldsBase, Test
         @test embed(M, A, A) === A
         @test manifold_dimension(M) == 9
         @test Weingarten(M, A, A, A) == zero(A)
+
+        @test is_point(M, rand(M))
+        @test is_point(M, rand(Random.MersenneTwister(), M))
+        @test is_vector(M, A, rand(M; vector_at=A))
+
+        @test get_coordinates(M, A, A, DefaultOrthonormalBasis()) == vec(A)
+        c = similar(vec(A))
+        get_coordinates!(M, c, A, A, DefaultOrthonormalBasis())
+        @test isapprox(c, vec(A))
+
+        @test get_vector(M, A, vec(A), DefaultOrthonormalBasis()) == A
+        D = similar(A)
+        get_vector!(M, D, A, vec(A), DefaultOrthonormalBasis())
+        @test isapprox(D, A)
     end
     @testset "Complex invertible matrices" begin
         @test repr(Mc) == "InvertibleMatrices(3, ℂ)"