More random randomness on random (no actually all) Lie groups (#706)

* Sketch of more randomness.
* remove redundant doc
* Fix typo
* Bump version, fix date.

---------

Co-authored-by: Mateusz Baran <mateuszbaran89@gmail.com>

NEWS.md

@@ -5,6 +5,14 @@ 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.9.14] – 2024-01-31
+
+### Added
+
+* `rand` on `UnitaryMatrices`
+* `rand` on arbitrary `GroupManifold`s and manifolds with `IsGroupManifold` trait
+  generating points and elements from the Lie algebra, respectively
+
 ## [0.9.13] – 2024-01-24
 
 ### Added

Project.toml

@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <seth.axen@gmail.com>", "Mateusz Baran <mateuszbaran89@gmail.com>", "Ronny Bergmann <manopt@ronnybergmann.net>", "Antoine Levitt <antoine.levitt@gmail.com>"]
-version = "0.9.13"
+version = "0.9.14"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

src/groups/GroupManifold.jl

@@ -1,15 +1,17 @@
 """
     GroupManifold{𝔽,M<:AbstractManifold{𝔽},O<:AbstractGroupOperation} <: AbstractDecoratorManifold{𝔽}
 
-Decorator for a smooth manifold that equips the manifold with a group operation, thus making
-it a Lie group. See [`IsGroupManifold`](@ref) for more details.
+Concrete decorator for a smooth manifold that equips the manifold with a group operation,
+thus making it a Lie group. See [`IsGroupManifold`](@ref) for more details.
 
 Group manifolds by default forward metric-related operations to the wrapped manifold.
 
-
 # Constructor
 
     GroupManifold(manifold, op)
+
+Define the group operation `op` acting on the manifold `manifold`, hence if `op` acts smoothly,
+this forms a Lie group.
 """
 struct GroupManifold{𝔽,M<:AbstractManifold{𝔽},O<:AbstractGroupOperation} <:
        AbstractDecoratorManifold{𝔽}
@@ -31,6 +33,10 @@ end
         IsExplicitDecorator(),
     )
 end
+# This could maybe even moved to ManifoldsBase?
+@inline function active_traits(f, ::AbstractRNG, M::AbstractDecoratorManifold, args...)
+    return active_traits(f, M, args...)
+end
 
 decorated_manifold(G::GroupManifold) = G.manifold
 
@@ -105,6 +111,49 @@ function is_vector(
     return is_vector(G.manifold, identity_element(G), X, false; error=error, kwargs...)
 end
 
+@doc raw"""
+    rand(::GroupManifold; vector_at=nothing, σ=1.0)
+    rand!(::GroupManifold, pX; vector_at=nothing, kwargs...)
+    rand(::TraitList{IsGroupManifold}, M; vector_at=nothing, σ=1.0)
+    rand!(TraitList{IsGroupManifold}, M, pX; vector_at=nothing, kwargs...)
+
+Compute a random point or tangent vector on a Lie group.
+
+For points this just means to generate a random point on the
+underlying manifold itself.
+
+For tangent vectors, an element in the Lie Algebra is generated.
+"""
+Random.rand(::GroupManifold; kwargs...)
+
+function Random.rand!(
+    T::TraitList{<:IsGroupManifold},
+    G::AbstractDecoratorManifold,
+    pX;
+    kwargs...,
+)
+    return rand!(T, Random.default_rng(), G, pX; kwargs...)
+end
+
+function Random.rand!(
+    ::TraitList{<:IsGroupManifold},
+    rng::AbstractRNG,
+    G::AbstractDecoratorManifold,
+    pX;
+    vector_at=nothing,
+    kwargs...,
+)
+    M = base_manifold(G)
+    if vector_at === nothing
+        # points we produce the same as on manifolds
+        rand!(rng, M, pX, kwargs...)
+    else
+        # tangent vectors are represented in the Lie algebra
+        # => materialize the identity and produce a tangent vector there
+        rand!(rng, M, pX; vector_at=identity_element(G), kwargs...)
+    end
+end
+
 Base.show(io::IO, G::GroupManifold) = print(io, "GroupManifold($(G.manifold), $(G.op))")
 
 function Statistics.var(M::GroupManifold, x::AbstractVector; kwargs...)

src/groups/general_unitary_groups.jl

@@ -287,15 +287,6 @@ function manifold_volume(M::GeneralUnitaryMultiplicationGroup)
     return manifold_volume(M.manifold)
 end
 
-function Random.rand!(G::GeneralUnitaryMultiplicationGroup, pX; kwargs...)
-    rand!(G.manifold, pX; kwargs...)
-    return pX
-end
-function Random.rand!(rng::AbstractRNG, G::GeneralUnitaryMultiplicationGroup, pX; kwargs...)
-    rand!(rng, G.manifold, pX; kwargs...)
-    return pX
-end
-
 function translate_diff!(
     G::GeneralUnitaryMultiplicationGroup,
     Y,

src/groups/power_group.jl

@@ -48,6 +48,32 @@ function ManifoldsBase._access_nested(
     return Identity(M.manifold)
 end
 
+# lower level methods are added instead of top level ones to not have to deal
+# with `Identity` disambiguation
+
+_compose!(G::PowerGroup, x, p, q) = _compose!(G.manifold, x, p, q)
+function _compose!(M::AbstractPowerManifold, x, p, q)
+    N = M.manifold
+    rep_size = representation_size(N)
+    for i in get_iterator(M)
+        compose!(
+            N,
+            _write(M, rep_size, x, i),
+            _read(M, rep_size, p, i),
+            _read(M, rep_size, q, i),
+        )
+    end
+    return x
+end
+function _compose!(M::PowerManifoldNestedReplacing, x, p, q)
+    N = M.manifold
+    rep_size = representation_size(N)
+    for i in get_iterator(M)
+        x[i...] = compose(N, _read(M, rep_size, p, i), _read(M, rep_size, q, i))
+    end
+    return x
+end
+
 @inline function active_traits(f, M::PowerGroup, args...)
     if is_metric_function(f)
         #pass to manifold by default - but keep Group Decorator for the retraction
@@ -174,32 +200,6 @@ function inv_diff!(G::PowerGroupNestedReplacing, Y, p, X)
     return Y
 end
 
-# lower level methods are added instead of top level ones to not have to deal
-# with `Identity` disambiguation
-
-_compose!(G::PowerGroup, x, p, q) = _compose!(G.manifold, x, p, q)
-function _compose!(M::AbstractPowerManifold, x, p, q)
-    N = M.manifold
-    rep_size = representation_size(N)
-    for i in get_iterator(M)
-        compose!(
-            N,
-            _write(M, rep_size, x, i),
-            _read(M, rep_size, p, i),
-            _read(M, rep_size, q, i),
-        )
-    end
-    return x
-end
-function _compose!(M::PowerManifoldNestedReplacing, x, p, q)
-    N = M.manifold
-    rep_size = representation_size(N)
-    for i in get_iterator(M)
-        x[i...] = compose(N, _read(M, rep_size, p, i), _read(M, rep_size, q, i))
-    end
-    return x
-end
-
 function translate!(G::PowerGroup, x, p, q, conv::ActionDirectionAndSide)
     return translate!(G.manifold, x, p, q, conv)
 end

src/manifolds/ProbabilitySimplexEuclideanMetric.jl

@@ -23,20 +23,44 @@ end
 @doc raw"""
     rand(::MetricManifold{ℝ,<:ProbabilitySimplex,<:EuclideanMetric}; vector_at=nothing, σ::Real=1.0)
 
-
 When `vector_at` is `nothing`, return a random (uniform) point `x` on the [`ProbabilitySimplex`](@ref) with the Euclidean metric
 manifold `M` by normalizing independent exponential draws to unit sum, see [Devroye:1986](@cite), Theorems 2.1 and 2.2 on p. 207 and 208, respectively.
-
 When `vector_at` is not `nothing`, return a (Gaussian) random vector from the tangent space
 ``T_{p}\mathrm{\Delta}^n``by shifting a multivariate Gaussian with standard deviation `σ`
 to have a zero component sum.
-
 """
 rand(::MetricManifold{ℝ,<:ProbabilitySimplex,<:EuclideanMetric}; σ::Real=1.0)
+function Random.rand(
+    M::MetricManifold{ℝ,<:ProbabilitySimplex,<:EuclideanMetric};
+    vector_at=nothing,
+    kwargs...,
+)
+    if vector_at === nothing
+        pX = allocate_result(M, rand)
+    else
+        pX = allocate_result(M, rand, vector_at)
+    end
+    rand!(M, pX; vector_at=vector_at, kwargs...)
+    return pX
+end
+function Random.rand(
+    rng::AbstractRNG,
+    M::MetricManifold{ℝ,<:ProbabilitySimplex,<:EuclideanMetric};
+    vector_at=nothing,
+    kwargs...,
+)
+    if vector_at === nothing
+        pX = allocate_result(M, rand)
+    else
+        pX = allocate_result(M, rand, vector_at)
+    end
+    rand!(rng, M, pX; vector_at=vector_at, kwargs...)
+    return pX
+end
 
 function Random.rand!(
     rng::AbstractRNG,
-    M::MetricManifold{ℝ,<:ProbabilitySimplex,<:EuclideanMetric},
+    ::MetricManifold{ℝ,<:ProbabilitySimplex,<:EuclideanMetric},
     pX;
     vector_at=nothing,
     σ=one(eltype(pX)),
@@ -51,6 +75,14 @@ function Random.rand!(
     return pX
 end
 
+function Random.rand!(
+    M::MetricManifold{ℝ,<:ProbabilitySimplex,<:EuclideanMetric},
+    pX;
+    kwargs...,
+)
+    return rand!(Random.default_rng(), M, pX; kwargs...)
+end
+
 @doc raw"""
     volume_density(::MetricManifold{ℝ,<:ProbabilitySimplex,<:EuclideanMetric}, p, X)
 

src/manifolds/Unitary.jl

@@ -20,6 +20,7 @@ The tangent spaces are given by
 
 But note that tangent vectors are represented in the Lie algebra, i.e. just using ``Y`` in
 the representation above.
+If you prefer the representation as `X` you can use the [`Stiefel`](@ref)`(n, n, ℂ)` manifold.
 
 # Constructor
 
@@ -130,6 +131,36 @@ function Random.rand(
     end
 end
 
+@doc raw"""
+    rand(::Unitary; vector_at=nothing, σ::Real=1.0)
+
+Generate a random point on the [`UnitaryMatrices`](@ref) manifold,
+if `vector_at` is nothing, by computing the QR decomposition of
+a ``n×x`` matrix.
+
+Generate a tangent vector at `vector_at` by projecting a normally
+distributed matrix onto the tangent space.
+"""
+rand(::UnitaryMatrices; σ::Real=1.0)
+
+function Random.rand!(
+    rng::AbstractRNG,
+    M::UnitaryMatrices,
+    pX;
+    vector_at=nothing,
+    σ::Real=one(real(eltype(pX))),
+)
+    n = get_parameter(M.size)[1]
+    if vector_at === nothing
+        A = σ * randn(rng, eltype(pX), n, n)
+        pX .= Matrix(qr(A).Q)
+    else
+        Z = σ * randn(rng, eltype(pX), size(pX))
+        project!(M, pX, vector_at, Z)
+    end
+    return pX
+end
+
 function Base.show(io::IO, ::UnitaryMatrices{TypeParameter{Tuple{n}},ℂ}) where {n}
     return print(io, "UnitaryMatrices($(n))")
 end

src/utils.jl

@@ -210,7 +210,7 @@ unrealify!(Y, X, ::typeof(ℝ), args...) = copyto!(Y, X)
 @doc raw"""
     symmetrize!(Y, X)
 
-Given a quare matrix `X` compute `1/2 .* (X' + X)` in place of `Y`
+Given a square matrix `X` compute `1/2 .* (X' + X)` in place of `Y`.
 """
 function symmetrize!(Y, X)
     Y .= (X' .+ X) ./ 2
@@ -220,7 +220,7 @@ end
 @doc raw"""
     symmetrize(X)
 
-Given a quare matrix `X` compute `1/2 .* (X' + X)`.
+Given a square matrix `X` compute `1/2 .* (X' + X)`.
 """
 function symmetrize(X)
     return (X' .+ X) ./ 2
@@ -229,7 +229,7 @@ end
 @doc raw"""
     vec2skew!(X, v, k)
 
-create a skew symmetric matrix inplace in `X` of size ``k×k`` from a vector `v`,
+Create a skew symmetric matrix in-place in `X` of size ``k×k`` from a vector `v`,
 for example for `v=[1,2,3]` and `k=3` this
 yields
 ````julia

test/manifolds/probability_simplex.jl

@@ -78,7 +78,7 @@ include("../header.jl")
                 test_vee_hat=false,
                 is_tangent_atol_multiplier=40.0,
                 default_inverse_retraction_method=nothing,
-                test_inplace=false,
+                test_inplace=true,
                 test_rand_point=true,
                 test_rand_tvector=true,
                 rand_tvector_atol_multiplier=40.0,

test/manifolds/unitary_matrices.jl

@@ -59,6 +59,22 @@ end
     r1 = exp(M, p, X, 1.0)
     @test isapprox(M, r, r1; atol=1e-10)
 
+    @testset "Projection" begin
+        M = UnitaryMatrices(2)
+        pE = [2im 0.0; 0.0 2im]
+        p = project(M, pE)
+        @test is_point(M, p; error=:error)
+        pE[2, 1] = 1.0
+        X = project(M, p, pE)
+        @test is_vector(M, p, X; error=:error)
+    end
+    @testset "Random" begin
+        M = UnitaryMatrices(2)
+        Random.seed!(23)
+        p = rand(M)
+        @test is_point(M, p; error=:error)
+        @test is_vector(M, p, rand(M; vector_at=p); error=:error)
+    end
     @testset "Riemannian Hessian" begin
         p = Matrix{Float64}(I, 2, 2)
         X = [0.0 3.0; -3.0 0.0]