diff --git a/NEWS.md b/NEWS.md index d0e50a50f9..ff94e315fc 100644 --- a/NEWS.md +++ b/NEWS.md @@ -14,12 +14,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * introduce `rand(:HamiltonianMatrices)` * extend `rand` to also `rand!` for `HamiltonianMatrices`, `SymplecticMatrices` and `SymplecticStiefel` * implement `riemannian_gradient` conversion for `SymplecticMatrices` and `SymplecticGrassmann` +* the new manifold of `MultinomialSymmetricPositiveDefinite` matrices +* `rand!` for `MultinomialDoublyStochastic` and `MultinomialSymmetric` ### Deprecated * Rename `Symplectic` to `SimplecticMatrices` in order to have a `Symplectic` wrapper for such matrices as well in the future for the next breaking change. * Rename `SymplecticMatrix` to `SymplecticElement` to clarify that it is the special matrix ``J_{2n}`` and not an arbitrary symplectic matrix. +### Fixed + +* a bug that cause `project` for tangent vectors to return wrong results on `MultinomialDoublyStochastic` + ## [0.9.12] – 2024-01-21 ### Fixed diff --git a/docs/make.jl b/docs/make.jl index c67ad99467..1268cf3af2 100755 --- a/docs/make.jl +++ b/docs/make.jl @@ -119,6 +119,7 @@ makedocs(; "Multinomial doubly stochastic matrices" => "manifolds/multinomialdoublystochastic.md", "Multinomial matrices" => "manifolds/multinomial.md", "Multinomial symmetric matrices" => "manifolds/multinomialsymmetric.md", + "Multinomial symmetric positive definite matrices" => "manifolds/multinomialsymmetricpositivedefinite.md", "Oblique manifold" => "manifolds/oblique.md", "Probability simplex" => "manifolds/probabilitysimplex.md", "Positive numbers" => "manifolds/positivenumbers.md", diff --git a/docs/src/manifolds/multinomialdoublystochastic.md b/docs/src/manifolds/multinomialdoublystochastic.md index 8cbaba852c..2b8cfcddc0 100644 --- a/docs/src/manifolds/multinomialdoublystochastic.md +++ b/docs/src/manifolds/multinomialdoublystochastic.md @@ -7,3 +7,8 @@ Order = [:type, :function] ``` ## Literature + +```@bibliography +Pages = ["multinomialdoublystochastic.md"] +Canonical=false +``` diff --git a/docs/src/manifolds/multinomialsymmetric.md b/docs/src/manifolds/multinomialsymmetric.md index d895045cdb..d8de6978ff 100644 --- a/docs/src/manifolds/multinomialsymmetric.md +++ b/docs/src/manifolds/multinomialsymmetric.md @@ -7,3 +7,8 @@ Order = [:type, :function] ``` ## Literature + +```@bibliography +Pages = ["multinomialsymmetric.md"] +Canonical=false +``` diff --git a/docs/src/manifolds/multinomialsymmetricpositivedefinite.md b/docs/src/manifolds/multinomialsymmetricpositivedefinite.md new file mode 100644 index 0000000000..14089f8594 --- /dev/null +++ b/docs/src/manifolds/multinomialsymmetricpositivedefinite.md @@ -0,0 +1,14 @@ +# Multinomial symmetric positive definite matrices + +```@autodocs +Modules = [Manifolds] +Pages = ["manifolds/MultinomialSymmetricPositiveDefinite.jl"] +Order = [:type, :function] +``` + +## Literature + +```@bibliography +Pages = ["multinomialsymmetricpositivedefinite.md"] +Canonical=false +``` diff --git a/docs/src/references.bib b/docs/src/references.bib index 98a24f12f6..3f2cc96041 100644 --- a/docs/src/references.bib +++ b/docs/src/references.bib @@ -275,6 +275,7 @@ @article{DouikHassibi:2019 TITLE = {Manifold Optimization Over the Set of Doubly Stochastic Matrices: A Second-Order Geometry}, JOURNAL = {IEEE Transactions on Signal Processing} } + # # E # ---------------------------------------------------------------------------------------- diff --git a/src/Manifolds.jl b/src/Manifolds.jl index 0ab1a05787..c0d5b3e746 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -252,6 +252,7 @@ using ManifoldsBase: IsEmbeddedSubmanifold, IsExplicitDecorator, LogarithmicInverseRetraction, + ManifoldDomainError, ManifoldsBase, NestedPowerRepresentation, NestedReplacingPowerRepresentation, @@ -426,6 +427,7 @@ include("manifolds/GeneralizedStiefel.jl") include("manifolds/Hyperbolic.jl") include("manifolds/MultinomialDoublyStochastic.jl") include("manifolds/MultinomialSymmetric.jl") +include("manifolds/MultinomialSymmetricPositiveDefinite.jl") include("manifolds/PositiveNumbers.jl") include("manifolds/ProjectiveSpace.jl") include("manifolds/SkewHermitian.jl") @@ -646,6 +648,7 @@ export Euclidean, MultinomialDoubleStochastic, MultinomialMatrices, MultinomialSymmetric, + MultinomialSymmetricPositiveDefinite, Oblique, OrthogonalMatrices, PositiveArrays, @@ -781,7 +784,7 @@ export CachedBasis, InducedBasis, ProjectedOrthonormalBasis # Errors on Manifolds -export ComponentManifoldError, CompositeManifoldError +export ComponentManifoldError, CompositeManifoldError, ManifoldDomainError # Functions on Manifolds export ×, ^, diff --git a/src/manifolds/Multinomial.jl b/src/manifolds/Multinomial.jl index 5d0dbc42b4..b62a658c25 100644 --- a/src/manifolds/Multinomial.jl +++ b/src/manifolds/Multinomial.jl @@ -5,8 +5,11 @@ The multinomial manifold consists of `m` column vectors, where each column is of `n` and unit norm, i.e. ````math -\mathcal{MN}(n,m) \coloneqq \bigl\{ p ∈ ℝ^{n×m}\ \big|\ p_{i,j} > 0 \text{ for all } i=1,…,n, j=1,…,m \text{ and } p^{\mathrm{T}}\mathbb{1}_m = \mathbb{1}_n\bigr\}, +\mathcal{MN}(n,m) \coloneqq \bigl\{ + p ∈ ℝ^{n×m}\ \big|\ p_{i,j} > 0 \text{ for all } i=1,…,n, j=1,…,m + \text{ and } p^{\mathrm{T}}\mathbb{1}_m = \mathbb{1}_n\bigr\}, ```` + where ``\mathbb{1}_k`` is the vector of length ``k`` containing ones. This yields exactly the same metric as @@ -85,6 +88,28 @@ function power_dimensions(M::MultinomialMatrices) return (m,) end +@doc raw""" + riemannian_gradient(M::MultinomialMatrices, p, Y; kwargs...) + +Let ``Y`` denote the Euclidean gradient of a function ``\tilde f`` defined in the +embedding neighborhood of `M`, then the Riemannian gradient is given by +Equation 5 of [DouikHassibi:2019](@cite) as + +```math + \operatorname{grad} f(p) = \proj_{T_p\mathcal M}(Y⊙p) +``` + +where ``⊙`` denotes the Hadamard or elementwise product. + +""" +riemannian_gradient(M::MultinomialMatrices, p, Y; kwargs...) + +function riemannian_gradient!(M::MultinomialMatrices, X, p, Y; kwargs...) + X .= p .* Y + project!(M, X, p, X) + return X +end + representation_size(M::MultinomialMatrices) = get_parameter(M.size) function Base.show(io::IO, ::MultinomialMatrices{TypeParameter{Tuple{n,m}}}) where {n,m} diff --git a/src/manifolds/MultinomialDoublyStochastic.jl b/src/manifolds/MultinomialDoublyStochastic.jl index 6ff2769a94..cd9de91962 100644 --- a/src/manifolds/MultinomialDoublyStochastic.jl +++ b/src/manifolds/MultinomialDoublyStochastic.jl @@ -3,7 +3,10 @@ A common type for manifolds that are doubly stochastic, for example by direct constraint [`MultinomialDoubleStochastic`](@ref) or by symmetry [`MultinomialSymmetric`](@ref), +or additionally by symmetric positive definiteness [`MultinomialSymmetricPositiveDefinite`](@ref) as long as they are also modeled as [`IsIsometricEmbeddedManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/decorator.html#ManifoldsBase.IsIsometricEmbeddedManifold). + +That way they share the inner product (just by restriction), and even the Riemannian gradient """ abstract type AbstractMultinomialDoublyStochastic <: AbstractDecoratorManifold{ℝ} end @@ -11,11 +14,46 @@ function active_traits(f, ::AbstractMultinomialDoublyStochastic, args...) return merge_traits(IsIsometricEmbeddedManifold()) end +@doc raw""" + representation_size(M::AbstractMultinomialDoublyStochastic) + +return the representation size of doubly stochastic matrices, whic are embedded +in the ``ℝ^{n×n}`` matrices and hence the answer here is `` +""" +function representation_size(M::AbstractMultinomialDoublyStochastic) + n = get_parameter(M.size)[1] + return (n, n) +end + +@doc raw""" + riemannian_gradient(M::AbstractMultinomialDoublyStochastic, p, Y; kwargs...) + +Let ``Y`` denote the Euclidean gradient of a function ``\tilde f`` defined in the +embedding neighborhood of `M`, then the Riemannian gradient is given by +Lemma 1 [DouikHassibi:2019](@cite) as + +```math + \operatorname{grad} f(p) = \proj_{T_p\mathcal M}(Y⊙p) +``` + +where ``⊙`` denotes the Hadamard or elementwise product, and the projection +is the projection onto the tangent space of the corresponding manifold. + +""" +riemannian_gradient(M::AbstractMultinomialDoublyStochastic, p, Y; kwargs...) + +function riemannian_gradient!(M::AbstractMultinomialDoublyStochastic, X, p, Y; kwargs...) + X .= p .* Y + project!(M, X, p, X) + return X +end + @doc raw""" MultinomialDoublyStochastic{T} <: AbstractMultinomialDoublyStochastic The set of doubly stochastic multinomial matrices consists of all ``n×n`` matrices with stochastic columns and rows, i.e. + ````math \begin{aligned} \mathcal{DP}(n) \coloneqq \bigl\{p ∈ ℝ^{n×n}\ \big|\ &p_{i,j} > 0 \text{ for all } i=1,…,n, j=1,…,m,\\ @@ -60,22 +98,14 @@ end Checks whether `p` is a valid point on the [`MultinomialDoubleStochastic`](@ref)`(n)` `M`, i.e. is a matrix with positive entries whose rows and columns sum to one. """ -function check_point( - M::MultinomialDoubleStochastic, - p::T; - atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), - kwargs..., -) where {T} +function check_point(M::MultinomialDoubleStochastic, p; kwargs...) n = get_parameter(M.size)[1] - r = sum(p, dims=2) - if !isapprox(norm(r - ones(n, 1)), 0.0; atol=atol, kwargs...) - return DomainError( - r, - "The point $(p) does not lie on $M, since its rows do not sum up to one.", - ) - end - return nothing + s = check_point(MultinomialMatrices(n, n), p; kwargs...) + !isnothing(s) && return s + s2 = check_point(MultinomialMatrices(n, n), p'; kwargs...) + return s2 end + @doc raw""" check_vector(M::MultinomialDoubleStochastic p, X; kwargs...) @@ -83,21 +113,12 @@ Checks whether `X` is a valid tangent vector to `p` on the [`MultinomialDoubleSt This means, that `p` is valid, that `X` is of correct dimension and sums to zero along any column or row. """ -function check_vector( - M::MultinomialDoubleStochastic, - p, - X::T; - atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), - kwargs..., -) where {T} - r = sum(X, dims=2) # check for stochastic rows - if !isapprox(norm(r), 0.0; atol=atol, kwargs...) - return DomainError( - r, - "The matrix $(X) is not a tangent vector to $(p) on $(M), since its rows do not sum up to zero.", - ) - end - return nothing +function check_vector(M::MultinomialDoubleStochastic, p, X; kwargs...) + n = get_parameter(M.size)[1] + s = check_vector(MultinomialMatrices(n, n), p, X; kwargs...) + !isnothing(s) && return s + s2 = check_vector(MultinomialMatrices(n, n), p, X'; kwargs...) + return s2 end function get_embedding(::MultinomialDoubleStochastic{TypeParameter{Tuple{n}}}) where {n} @@ -134,11 +155,14 @@ end Project `Y` onto the tangent space at `p` on the [`MultinomialDoubleStochastic`](@ref) `M`, return the result in `X`. The formula reads + ````math \operatorname{proj}_p(Y) = Y - (α\mathbf{1}_n^{\mathrm{T}} + \mathbf{1}_nβ^{\mathrm{T}}) ⊙ p, ```` + where ``⊙`` denotes the Hadamard or elementwise product and ``\mathbb{1}_n`` is the vector of length ``n`` containing ones. The two vectors ``α,β ∈ ℝ^{n×n}`` are computed as a solution (typically using the left pseudo inverse) of + ````math \begin{pmatrix} I_n & p\\p^{\mathrm{T}} & I_n \end{pmatrix} \begin{pmatrix} α\\ β\end{pmatrix} @@ -152,8 +176,8 @@ project(::MultinomialDoubleStochastic, ::Any, ::Any) function project!(M::MultinomialDoubleStochastic, X, p, Y) n = get_parameter(M.size)[1] - ζ = [I p; p I] \ [sum(Y, dims=2); sum(Y, dims=1)'] # Formula (25) from 1802.02628 - return X .= Y .- (repeat(ζ[1:n], 1, 3) .+ repeat(ζ[(n + 1):end]', 3, 1)) .* p + ζ = [I p; p' I] \ [sum(Y, dims=2); sum(Y, dims=1)'] # Formula (25) from 1802.02628 + return X .= Y .- (repeat(ζ[1:n], 1, n) .+ repeat(ζ[(n + 1):end]', n, 1)) .* p end @doc raw""" @@ -201,15 +225,45 @@ function project!( return q end -function representation_size(M::MultinomialDoubleStochastic) - n = get_parameter(M.size)[1] - return (n, n) +@doc raw""" + rand(::MultinomialDoubleStochastic; vector_at=nothing, σ::Real=1.0, kwargs...) + +Generate random points on the [`MultinomialDoubleStochastic`](@ref) manifold +or tangent vectors at the point `vector_at` if that is not `nothing`. + +Let ``n×n`` denote the matrix dimension of the [`MultinomialDoubleStochastic`](@ref). + +When `vector_at` is nothing, this is done by generating a random matrix`rand(n,n)` +with positive entries and projecting it onto the manifold. The `kwargs...` are +passed to this projection. + +When `vector_at` is not `nothing`, a random matrix in the ambient space is generated +and projected onto the tangent space +""" +rand(::MultinomialDoubleStochastic; σ::Real=1.0) + +function Random.rand!( + rng::AbstractRNG, + M::MultinomialDoubleStochastic, + pX; + vector_at=nothing, + σ::Real=one(real(eltype(pX))), + kwargs..., +) + rand!(rng, pX) + pX .*= σ + if vector_at === nothing + project!(M, pX, pX; kwargs...) + else + project!(M, pX, vector_at, pX) + end + return pX end @doc raw""" retract(M::MultinomialDoubleStochastic, p, X, ::ProjectionRetraction) -compute a projection based retraction by projecting $p\odot\exp(X⨸p)$ back onto the manifold, +compute a projection based retraction by projecting ``p\odot\exp(X⨸p)`` back onto the manifold, where ``⊙,⨸`` are elementwise multiplication and division, respectively. Similarly, ``\exp`` refers to the elementwise exponentiation. """ diff --git a/src/manifolds/MultinomialSymmetric.jl b/src/manifolds/MultinomialSymmetric.jl index 470f6bb779..c0d11bfc6f 100644 --- a/src/manifolds/MultinomialSymmetric.jl +++ b/src/manifolds/MultinomialSymmetric.jl @@ -1,5 +1,5 @@ @doc raw""" - MultinomialSymmetric{T} <: AbstractMultinomialDoublyStochastic{N} + MultinomialSymmetric{T} <: AbstractMultinomialDoublyStochastic The multinomial symmetric matrices manifold consists of all symmetric ``n×n`` matrices with positive entries such that each column sums to one, i.e. @@ -56,7 +56,10 @@ i.e. is a symmetric matrix with positive entries whose rows sum to one. """ function check_point(M::MultinomialSymmetric, p; kwargs...) n = get_parameter(M.size)[1] - return check_point(SymmetricMatrices(n, ℝ), p) + s = check_point(SymmetricMatrices(n, ℝ), p; kwargs...) + !isnothing(s) && return s + s2 = check_point(MultinomialMatrices(n, n), p; kwargs...) + return s2 end @doc raw""" check_vector(M::MultinomialSymmetric p, X; kwargs...) @@ -67,7 +70,10 @@ along any row. """ function check_vector(M::MultinomialSymmetric, p, X; kwargs...) n = get_parameter(M.size)[1] - return check_vector(SymmetricMatrices(n, ℝ), p, X; kwargs...) + s = check_vector(SymmetricMatrices(n, ℝ), p, X; kwargs...) + !isnothing(s) && return s + s2 = check_vector(MultinomialMatrices(n, n), p, X) + return s2 end embed!(::MultinomialSymmetric, q, p) = copyto!(q, p) @@ -106,7 +112,9 @@ end project(M::MultinomialSymmetric, p, Y) Project `Y` onto the tangent space at `p` on the [`MultinomialSymmetric`](@ref) `M`, return the result in `X`. -The formula reads + +The formula from [DouikHassibi:2019](@cite), Sec. VI reads + ````math \operatorname{proj}_p(Y) = Y - (α\mathbf{1}_n^{\mathrm{T}} + \mathbf{1}_n α^{\mathrm{T}}) ⊙ p, ```` @@ -116,13 +124,47 @@ The two vector ``α ∈ ℝ^{n×n}`` is given by solving (I_n+p)α = Y\mathbf{1}, ```` where ``I_n`` is teh ``n×n`` unit matrix and ``\mathbf{1}_n`` is the vector of length ``n`` containing ones. - """ project(::MultinomialSymmetric, ::Any, ::Any) -function project!(::MultinomialSymmetric, X, p, Y) +function project!(M::MultinomialSymmetric, X, p, Y) + n = get_parameter(M.size)[1] α = (I + p) \ sum(Y, dims=2) # Formula (49) from 1802.02628 - return X .= Y .- (repeat(α, 1, 3) .+ repeat(α', 3, 1)) .* p + return X .= Y .- (repeat(α, 1, n) .+ repeat(α', n, 1)) .* p +end + +@doc raw""" + rand(::MultinomialSymmetric; vector_at=nothing, σ::Real=1.0, kwargs...) + +Generate random points on the [`MultinomialSymmetric`](@ref) manifold +or tangent vectors at the point `vector_at` if that is not `nothing`. + +Let ``n×n`` denote the matrix dimension of the [`MultinomialSymmetric`](@ref). + +When `vector_at` is nothing, this is done by generating a random matrix `rand(n, n)` +with positive entries and projecting it onto the manifold. The `kwargs...` are +passed to this projection. + +When `vector_at` is not `nothing`, a random matrix in the ambient space is generated +and projected onto the tangent space +""" +rand(::MultinomialSymmetric; σ::Real=1.0) + +function Random.rand!( + rng::AbstractRNG, + M::MultinomialSymmetric, + pX; + vector_at=nothing, + kwargs..., +) + n = get_parameter(M.size)[1] + rand!(rng, SymmetricMatrices(n), pX; kwargs...) + if vector_at === nothing + project!(M, pX, pX) + else + project!(M, pX, vector_at, pX) + end + return pX end function representation_size(M::MultinomialSymmetric) @@ -143,6 +185,35 @@ function retract_project!(M::MultinomialSymmetric, q, p, X, t::Number) return project!(M, q, p .* exp.(t .* X ./ p)) end +@doc raw""" + Y = riemannian_Hessian(M::MultinomialSymmetric, p, G, H, X) + riemannian_Hessian!(M::MultinomialSymmetric, Y, p, G, H, X) + +Compute the Riemannian Hessian ``\operatorname{Hess} f(p)[X]`` given the +Euclidean gradient ``∇ f(\tilde p)`` in `G` and the Euclidean Hessian ``∇^2 f(\tilde p)[\tilde X]`` in `H`, +where ``\tilde p, \tilde X`` are the representations of ``p,X`` in the embedding,. + +The Riemannian Hessian can be computed as stated in Corollary 3 [DouikHassibi:2019](@cite). +""" +riemannian_Hessian(M::MultinomialSymmetric, p, G, H, X) + +function riemannian_Hessian!(M::MultinomialSymmetric, Y, p, G, H, X) + # The notation here is the same as in (53) DouikHassibi:2019 + # with the small change their X is our p their ξ_X is our X , Hessf is H, Gradf is G + n = get_parameter(M.size)[1] + ov = ones(n) # \bf 1 + I_p = lu(I + p) + γ = G .* p + α = I_p \ (γ * ov) + α_sq = (repeat(α, 1, n) .+ repeat(α', n, 1)) + δ = γ .- α_sq .* p + γ_dot = H .* p + G .* X + α_dot = (I_p \ γ_dot .- (I_p \ X) * (I_p \ γ)) * ov + δ_dot = γ_dot .- (repeat(α_dot, 1, n) .+ repeat(α_dot', n, 1)) .* p .- α_sq .* X + project!(M, Y, p, δ_dot .- 0.5 * ((δ .* X) ./ p)) + return Y +end + function Base.show(io::IO, ::MultinomialSymmetric{TypeParameter{Tuple{n}}}) where {n} return print(io, "MultinomialSymmetric($(n))") end diff --git a/src/manifolds/MultinomialSymmetricPositiveDefinite.jl b/src/manifolds/MultinomialSymmetricPositiveDefinite.jl new file mode 100644 index 0000000000..faf8198650 --- /dev/null +++ b/src/manifolds/MultinomialSymmetricPositiveDefinite.jl @@ -0,0 +1,123 @@ +@doc raw""" + MultinomialSymmetricPositiveDefinite <: AbstractMultinomialDoublyStochastic + +The symmetric positive definite multinomial matrices manifold consists of all +symmetric ``n×n`` matrices with positive eigenvalues, and +positive entries such that each column sums to one, i.e. + +````math +\begin{aligned} +\mathcal{SP}^+(n) \coloneqq \bigl\{ + p ∈ ℝ^{n×n}\ \big|\ &p_{i,j} > 0 \text{ for all } i=1,…,n, j=1,…,m,\\ +& p^\mathrm{T} = p,\\ +& p\mathbf{1}_n = \mathbf{1}_n\\ +a^\mathrm{T}pa > 0 \text{ for all } a ∈ ℝ^{n}\backslash\{\mathbf{0}_n\} +\bigr\}, +\end{aligned} +```` + +where ``\mathbf{1}_n`` and ``\mathbr{0}_n`` are the vectors of length ``n`` +containing ones and zeros, respectively. More details about this manifold can be found in +[DouikHassibi:2019](@cite). + +# Constructor + + MultinomialSymmetricPositiveDefinite(n) + +Generate the manifold of matrices ``\mathbb R^{n×n}`` that are symmetric, positive definite, and doubly stochastic. +""" +struct MultinomialSymmetricPositiveDefinite{T} <: AbstractMultinomialDoublyStochastic + size::T +end + +function MultinomialSymmetricPositiveDefinite(n::Int; parameter::Symbol=:type) + size = wrap_type_parameter(parameter, (n,)) + return MultinomialSymmetricPositiveDefinite{typeof(size)}(size) +end + +function check_point(M::MultinomialSymmetricPositiveDefinite, p; kwargs...) + # Multinomial checked first via embedding + n = get_parameter(M.size)[1] + s = check_point(SymmetricPositiveDefinite(n), p; kwargs...) + !isnothing(s) && + return ManifoldDomainError("The point $(p) does not lie on the $(M).", s) + return nothing +end + +function check_vector(M::MultinomialSymmetricPositiveDefinite, p, X; kwargs...) + # Multinomial checked first via embedding + n = get_parameter(M.size)[1] + s = check_vector(SymmetricPositiveDefinite(n), p, X; kwargs...) + !isnothing(s) && return ManifoldDomainError( + "The vector $(X) is not a tangent vector to $(p) on $(M)", + s, + ) + return nothing +end + +function get_embedding( + ::MultinomialSymmetricPositiveDefinite{TypeParameter{Tuple{n}}}, +) where {n} + return MultinomialMatrices(n, n) +end +function get_embedding(M::MultinomialSymmetricPositiveDefinite{Tuple{Int}}) + n = get_parameter(M.size)[1] + return MultinomialMatrices(n, n; parameter=:field) +end + +""" + Random.rand!( + rng::AbstractRNG, + M::MultinomialSymmetricPositiveDefinite, + p::AbstractMatrix, + ) + +Generate a random point on [`MultinomialSymmetricPositiveDefinite`](@ref) manifold. +The steps are as follows: +1. Generate a random [totally positive matrix](https://en.wikipedia.org/wiki/Totally_positive_matrix) + a. Construct a vector `L` of `n` random positive increasing real numbers. + b. Construct the [Vandermonde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix) + `V` based on the sequence `L`. + c. Perform LU factorization of `V` in such way that both L and U components have + positive elements. + d. Convert the LU factorization into LDU factorization by taking the diagonal of U + and dividing U by it, `V=LDU`. + e. Construct a new matrix `R = UDL` which is totally positive. +2. Project the totally positive matrix `R` onto the manifold of [`MultinomialDoubleStochastic`](@ref) + matrices. +3. Symmetrize the projected matrix and return the result. + +This method roughly follows the procedure described in https://math.stackexchange.com/questions/2773460/how-to-generate-a-totally-positive-matrix-randomly-using-software-like-maple +""" +function Random.rand!( + rng::AbstractRNG, + M::MultinomialSymmetricPositiveDefinite, + p::AbstractMatrix, +) + n = get_parameter(M.size)[1] + L = sort(exp.(randn(rng, n))) + V = reduce(hcat, map(xi -> [xi^k for k in 0:(n - 1)], L))' + @static if VERSION < v"1.7" + Vlu = lu(V, Val(false)) + else + Vlu = lu(V, LinearAlgebra.RowNonZero()) + end + dm = Diagonal(Vlu.U) + uutd = dm \ Vlu.U + random_totally_positive = uutd * dm * Vlu.L + MMDS = MultinomialDoubleStochastic(n) + ds = project(MMDS, random_totally_positive) + p .= (ds .+ ds') ./ 2 + return p +end + +function Base.show( + io::IO, + ::MultinomialSymmetricPositiveDefinite{TypeParameter{Tuple{n}}}, +) where {n} + return print(io, "MultinomialSymmetricPositiveDefinite($(n))") +end +function Base.show(io::IO, M::MultinomialSymmetricPositiveDefinite{Tuple{Int}}) + n = get_parameter(M.size)[1] + return print(io, "MultinomialSymmetricPositiveDefinite($(n); parameter=:field)") +end diff --git a/src/manifolds/Symmetric.jl b/src/manifolds/Symmetric.jl index 5f0a982e24..686ea65f6d 100644 --- a/src/manifolds/Symmetric.jl +++ b/src/manifolds/Symmetric.jl @@ -226,6 +226,18 @@ project(::SymmetricMatrices, ::Any, ::Any) project!(M::SymmetricMatrices, Y, p, X) = (Y .= (X .+ transpose(X)) ./ 2) +function Random.rand!( + rng::AbstractRNG, + M::SymmetricMatrices, + pX; + σ::Real=one(real(eltype(pX))), + kwargs..., +) + rand!(rng, pX) + pX .= (σ / (2 * norm(pX))) .* (pX + pX') + return pX +end + function Base.show(io::IO, ::SymmetricMatrices{TypeParameter{Tuple{n}},F}) where {n,F} return print(io, "SymmetricMatrices($(n), $(F))") end diff --git a/test/manifolds/multinomial_doubly_stochastic.jl b/test/manifolds/multinomial_doubly_stochastic.jl index 22e55e77d9..fb1177413b 100644 --- a/test/manifolds/multinomial_doubly_stochastic.jl +++ b/test/manifolds/multinomial_doubly_stochastic.jl @@ -10,14 +10,14 @@ include("../header.jl") @test is_vector(M, p, X) pf1 = [0.1 0.9 0.1; 0.1 0.9 0.1; 0.1 0.1 0.9] #not sum 1 @test_throws ManifoldDomainError is_point(M, pf1; error=:error) - pf2r = [0.1 0.9 0.1; 0.8 0.05 0.15; 0.1 0.05 0.75] - @test_throws DomainError is_point(M, pf2r; error=:error) + pf2r = [0.1 0.9 0.1; 0.8 0.05 0.15; 0.1 0.05 0.75] # sum(pf2r,dims=2)[3] is 0.9 + @test_throws CompositeManifoldError is_point(M, pf2r; error=:error) @test_throws ManifoldDomainError is_point(M, pf2r'; error=:error) pf3 = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] # contains nonpositive entries @test_throws ManifoldDomainError is_point(M, pf3; error=:error) Xf2c = [-0.1 0.0 0.1; -0.2 0.1 0.1; 0.2 -0.1 -0.1] #nonzero columns @test_throws ManifoldDomainError is_vector(M, p, Xf2c; error=:error) - @test_throws DomainError is_vector(M, p, Xf2c'; error=:error) + @test_throws CompositeManifoldError is_vector(M, p, Xf2c'; error=:error) @test representation_size(M) == (3, 3) @test !is_flat(M) pE = similar(p) @@ -68,4 +68,19 @@ include("../header.jl") @test repr(M) == "MultinomialDoubleStochastic(3; parameter=:field)" @test get_embedding(M) === MultinomialMatrices(3, 3; parameter=:field) end + @testset "random" begin + Random.seed!(42) + p = rand(M) + @test is_point(M, p) + X = rand(M; vector_at=p) + @test is_vector(M, p, X; atol=1e-9) + end + @testset "Riemannian Gradient" begin + M = MultinomialDoubleStochastic(3) + p = ones(3, 3) ./ 3 + Y = [1.0 -1.0 0.0; 0.0 0.0 1.0; -1.0 1.0 1.1] + G = project(M, p, p .* Y) + X = riemannian_gradient(M, p, Y) + @test isapprox(M, p, G, X) + end end diff --git a/test/manifolds/multinomial_matrices.jl b/test/manifolds/multinomial_matrices.jl index f23fa4e760..4fd20a46c9 100644 --- a/test/manifolds/multinomial_matrices.jl +++ b/test/manifolds/multinomial_matrices.jl @@ -55,4 +55,12 @@ include("../header.jl") M = MultinomialMatrices(3, 2; parameter=:field) @test repr(M) == "MultinomialMatrices(3, 2; parameter=:field)" end + @testset "Riemannian Gradient" begin + M = MultinomialMatrices(3, 2) + p = [0.5 0.4 0.1; 0.5 0.4 0.1]' + Y = [1.0 -1.0; 0.0 0.0; -1.0 1.0] + G = project(M, p, p .* Y) + X = riemannian_gradient(M, p, Y) + @test isapprox(M, p, G, X) + end end diff --git a/test/manifolds/multinomial_spd.jl b/test/manifolds/multinomial_spd.jl new file mode 100644 index 0000000000..70a6027a22 --- /dev/null +++ b/test/manifolds/multinomial_spd.jl @@ -0,0 +1,39 @@ +include("../header.jl") + +@testset "Multinomial symmetric positive definite matrices" begin + @testset "Basics" begin + M = MultinomialSymmetricPositiveDefinite(3) + Mf = MultinomialSymmetricPositiveDefinite(3; parameter=:field) + @test repr(M) == "MultinomialSymmetricPositiveDefinite(3)" + @test repr(Mf) == "MultinomialSymmetricPositiveDefinite(3; parameter=:field)" + @test get_embedding(M) == MultinomialMatrices(3, 3) + @test get_embedding(Mf) == MultinomialMatrices(3, 3; parameter=:field) + # + # Checks + # (a) Points + p = [0.6 0.2 0.2; 0.2 0.6 0.2; 0.2 0.2 0.6] + @test is_point(M, p; error=:error) + # Symmetric but does not sum to 1 + pf1 = zeros(3, 3) + @test_throws ManifoldDomainError is_point(M, pf1; error=:error) + # in theory this is not spd since it has an EV 0 but numerically it is + pf2 = [0.3 0.4 0.3; 0.4 0.2 0.4; 0.3 0.4 0.3] + # Multinomial but not symmetric + pf3 = [0.2 0.3 0.5; 0.4 0.2 0.4; 0.4 0.5 0.1] + @test_throws ManifoldDomainError is_point(M, pf3; error=:error) + # (b) Tangent vectors + X = [1.0 -0.5 -0.5; -0.5 1.0 -0.5; -0.5 -0.5 1.0] + @test is_vector(M, p, X; error=:error) + Xf1 = ones(3, 3) # Symmetric but does not sum to zero + @test_throws ManifoldDomainError is_vector(M, p, Xf1; error=:error) + #sums to zero but not symmetric + Xf2 = [-0.5 0.3 0.2; 0.2 -0.5 0.3; 0.3 0.2 -0.5] + @test_throws ManifoldDomainError is_vector(M, p, Xf2; error=:error) + end + + @testset "Random" begin + q = zeros(3, 3) + M = MultinomialSymmetricPositiveDefinite(3) + @test is_point(M, rand!(MersenneTwister(), M, q)) + end +end diff --git a/test/manifolds/multinomial_symmetric.jl b/test/manifolds/multinomial_symmetric.jl index fcf0a2360e..fe7961ec57 100644 --- a/test/manifolds/multinomial_symmetric.jl +++ b/test/manifolds/multinomial_symmetric.jl @@ -73,4 +73,22 @@ include("../header.jl") @test repr(M) == "MultinomialSymmetric(3; parameter=:field)" @test get_embedding(M) === MultinomialMatrices(3, 3; parameter=:field) end + @testset "random" begin + Random.seed!(42) + p = rand(M) + @test is_point(M, p) + X = rand(M; vector_at=p) + @test is_vector(M, p, X) + end + @testset "Hessian call" begin + p = ones(3, 3) ./ 3 + Y = one(p) + G = zero(p) + H = 0.5 * one(p) + X = riemannian_Hessian(M, p, G, H, Y) + X2 = similar(X) + riemannian_Hessian!(M, X2, p, G, H, Y) + @test isapprox(M, p, X, X2) + @test is_vector(M, p, X) + end end diff --git a/test/runtests.jl b/test/runtests.jl index bc740b3112..1764273c24 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -154,6 +154,7 @@ end include_test("manifolds/lorentz.jl") include_test("manifolds/multinomial_doubly_stochastic.jl") include_test("manifolds/multinomial_symmetric.jl") + include_test("manifolds/multinomial_spd.jl") include_test("manifolds/positive_numbers.jl") include_test("manifolds/probability_simplex.jl") include_test("manifolds/projective_space.jl")