Multinomial doubly stochastic matrices

AbstractMultinomialDoublyStochastic <: AbstractDecoratorManifold{ℝ}

A common type for manifolds that are doubly stochastic. This can be by direct constraint MultinomialDoubleStochastic, by symmetry MultinomialSymmetric, or additionally by symmetric positive definiteness MultinomialSymmetricPositiveDefinite, as long as they are also modeled as IsometricallyEmbeddedManifoldType.

That way they share the inner product (just by restriction), and even the Riemannian gradient

MultinomialDoubleStochastic{T} <: AbstractMultinomialDoublyStochastic

The set of doubly stochastic multinomial matrices consists of all $n×n$ matrices with stochastic columns and rows, i.e.

\[\begin{aligned} \mathcal{DP}(n) \coloneqq \bigl\{p ∈ ℝ^{n×n}\ \big|\ &p_{i,j} > 0 \text{ for all } i=1,…,n, j=1,…,n,\\ & p\mathbf{1}_n = p^{\mathrm{T}}\mathbf{1}_n = \mathbf{1}_n \bigr\}, \end{aligned}\]

where $\mathbf{1}_n$ is the vector of length $n$ containing ones.

The tangent space is given in Proposition 3.1 of [Dou20] as

\[T_p\mathcal{DP}(n) \coloneqq \bigl\{ X ∈ ℝ^{n×n}\ \big|\ X\mathbf{1}_n = X^{\mathrm{T}}\mathbf{1}_n = \mathbf{0}_n \bigr\},\]

where $\mathbf{0}_n$ is the vector of length $n$ containing zeros.

More details can be found in Section III [DH19].

Constructor

MultinomialDoubleStochastic(n::Int; parameter::Symbol=:type)

Generate the manifold of matrices $ℝ^{n×n}$ that are doubly stochastic.

rand(::MultinomialDoubleStochastic; vector_at=nothing, σ::Real=1.0, kwargs...)

Generate random points on the MultinomialDoubleStochastic manifold or tangent vectors at the point vector_at if that is not nothing.

Let $n×n$ denote the matrix dimension of the MultinomialDoubleStochastic.

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; σ has no effect in this case, since the projection normalizes rows and columns. The kwargs... are passed to this projection.

When vector_at is not nothing, a random matrix in the ambient space is generated, multiplied by σ and projected onto the tangent space.

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 [DH19] as

\[ \operatorname{grad} f(p) = \operatorname{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.

check_point(M::MultinomialDoubleStochastic, p)

Checks whether p is a valid point on the MultinomialDoubleStochastic(n) M, i.e. is a matrix with positive entries whose rows and columns sum to one.

check_vector(M::MultinomialDoubleStochastic p, X; kwargs...)

Checks whether X is a valid tangent vector to p on the MultinomialDoubleStochastic M. This means, that p is valid, that X is of correct dimension and sums to zero along any column or row.

default_retraction_method(M::AbstractMultinomialDoublyStochastic)

Return ProjectionRetraction, the retraction that the doubly stochastic multinomial manifolds implement.

is_flat(::MultinomialDoubleStochastic)

Return false. MultinomialDoubleStochastic is not a flat manifold.

manifold_dimension(M::MultinomialDoubleStochastic)

returns the dimension of the MultinomialDoubleStochastic manifold namely

\[\operatorname{dim}_{\mathcal{DP}(n)} = (n-1)^2.\]

project(
    M::AbstractMultinomialDoublyStochastic,
    p;
    maxiter = 1000,
    tolerance = eps(real(float(eltype(p)))),
    warn_nonconvergence = true,
)

project a matrix p with positive entries applying Sinkhorn's algorithm. Note that this project method – different from the usual case, accepts keywords. If the algorithm does not reach tolerance within maxiter iterations, a warning is issued unless warn_nonconvergence is set to false.

project(M::MultinomialDoubleStochastic, p, Y)

Project Y onto the tangent space at p on the MultinomialDoubleStochastic M, return the result in X. The formula reads

\[ \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$ are computed with an LU decomposition, or with the pseudo inverse if that fails, as a solution of

\[ \begin{pmatrix} I_n & p\\p^{\mathrm{T}} & I_n \end{pmatrix} \begin{pmatrix} α\\ β\end{pmatrix} = \begin{pmatrix} Y\mathbf{1}\\Y^{\mathrm{T}}\mathbf{1}\end{pmatrix},\]

where $I_n$ is the $n×n$ unit matrix and $\mathbf{1}_n$ is the vector of length $n$ containing ones. This system has infinitely many solutions, and all of them yield the same projection, see [DH19] and Theorem 3.1 in [Dou20].

representation_size(M::AbstractMultinomialDoublyStochastic)

return the representation size of doubly stochastic matrices, which are embedded in the $ℝ^{n×n}$ matrices and hence the answer here is $(n,n)$.

retract(M::MultinomialDoubleStochastic, p, X, ::ProjectionRetraction)

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.