Merge pull request #1097 from johnczito/add_singular_wishart

Add singular branch of the Wishart

src/matrix/wishart.jl

@@ -1,58 +1,70 @@
-# Wishart distribution
-#
-#   following the Wikipedia parameterization
-#
 """
     Wishart(ν, S)
 ```julia
-ν::Real           degrees of freedom (greater than p - 1)
+ν::Real           degrees of freedom (whole number or a real number greater than p - 1)
 S::AbstractPDMat  p x p scale matrix
 ```
 The [Wishart distribution](http://en.wikipedia.org/wiki/Wishart_distribution)
-generalizes the gamma distribution to ``p\\times p`` real, positive definite
-matrices ``\\mathbf{H}``. If ``\\mathbf{H}\\sim \\textrm{W}_p(\\nu,\\mathbf{S})``,
-then its probability density function is
+generalizes the gamma distribution to ``p\\times p`` real, positive semidefinite
+matrices ``\\mathbf{H}``.
+
+If ``\\nu>p-1``, then ``\\mathbf{H}\\sim \\textrm{W}_p(\\nu, \\mathbf{S})``
+has rank ``p`` and its probability density function is
 
 ```math
 f(\\mathbf{H};\\nu,\\mathbf{S}) = \\frac{1}{2^{\\nu p/2} \\left|\\mathbf{S}\\right|^{\\nu/2} \\Gamma_p\\left(\\frac {\\nu}{2}\\right ) }{\\left|\\mathbf{H}\\right|}^{(\\nu-p-1)/2} e^{-(1/2)\\operatorname{tr}(\\mathbf{S}^{-1}\\mathbf{H})}.
 ```
 
-If ``\\nu`` is an integer, then a random matrix ``\\mathbf{H}`` given by
+If ``\\nu\\leq p-1``, then ``\\mathbf{H}`` is rank ``\\nu`` and it has
+a density with respect to a suitably chosen volume element on the space of
+positive semidefinite matrices. See [here](https://doi.org/10.1214/aos/1176325375).
+
+For integer ``\\nu``, a random matrix given by
 
 ```math
-\\mathbf{H} = \\mathbf{X}\\mathbf{X}^{\\rm{T}}, \\quad\\mathbf{X} \\sim \\textrm{MN}_{p,\\nu}(\\mathbf{0}, \\mathbf{S}, \\mathbf{I}_{\\nu})
+\\mathbf{H} = \\mathbf{X}\\mathbf{X}^{\\rm{T}},
+\\quad\\mathbf{X} \\sim \\textrm{MN}_{p,\\nu}(\\mathbf{0}, \\mathbf{S}, \\mathbf{I}_{\\nu})
 ```
 
-has ``\\mathbf{H}\\sim \\textrm{W}_p(\\nu, \\mathbf{S})``. For non-integer degrees of freedom,
-Wishart matrices can be generated via the [Bartlett decomposition](https://en.wikipedia.org/wiki/Wishart_distribution#Bartlett_decomposition).
+has ``\\mathbf{H}\\sim \\textrm{W}_p(\\nu, \\mathbf{S})``.
+For non-integer ``\\nu``, Wishart matrices can be generated via the
+[Bartlett decomposition](https://en.wikipedia.org/wiki/Wishart_distribution#Bartlett_decomposition).
 """
-struct Wishart{T<:Real, ST<:AbstractPDMat} <: ContinuousMatrixDistribution
-    df::T     # degree of freedom
-    S::ST     # the scale matrix
-    logc0::T  # the logarithm of normalizing constant in pdf
+struct Wishart{T<:Real, ST<:AbstractPDMat, R<:Integer} <: ContinuousMatrixDistribution
+    df::T          # degree of freedom
+    S::ST          # the scale matrix
+    logc0::T       # the logarithm of normalizing constant in pdf
+    rank::R        # rank of a sample
+    singular::Bool # singular of nonsingular wishart?
 end
 
 #  -----------------------------------------------------------------------------
 #  Constructors
 #  -----------------------------------------------------------------------------
 
-function Wishart(df::T, S::AbstractPDMat{T}) where T<:Real
+function Wishart(df::T, S::AbstractPDMat{T}, warn::Bool = true) where T<:Real
+    df > 0 || throw(ArgumentError("df must be positive. got $(df)."))
     p = dim(S)
-    df > p - 1 || throw(ArgumentError("df should be greater than dim - 1."))
-    logc0 = wishart_logc0(df, S)
+    rnk = p
+    singular = df <= p - 1
+    if singular
+        isinteger(df) || throw(ArgumentError("singular df must be an integer. got $(df)."))
+        rnk = convert(Integer, df)
+        warn && @warn("got df <= dim - 1; returning a singular Wishart")
+    end
+    logc0 = wishart_logc0(df, S, rnk)
     R = Base.promote_eltype(T, logc0)
     prom_S = convert(AbstractArray{T}, S)
-    Wishart{R, typeof(prom_S)}(R(df), prom_S, R(logc0))
+    Wishart{R, typeof(prom_S), typeof(rnk)}(R(df), prom_S, R(logc0), rnk, singular)
 end
 
-function Wishart(df::Real, S::AbstractPDMat)
+function Wishart(df::Real, S::AbstractPDMat, warn::Bool = true)
     T = Base.promote_eltype(df, S)
-    Wishart(T(df), convert(AbstractArray{T}, S))
+    Wishart(T(df), convert(AbstractArray{T}, S), warn)
 end
 
-Wishart(df::Real, S::Matrix) = Wishart(df, PDMat(S))
-
-Wishart(df::Real, S::Cholesky) = Wishart(df, PDMat(S))
+Wishart(df::Real, S::Matrix, warn::Bool = true) = Wishart(df, PDMat(S), warn)
+Wishart(df::Real, S::Cholesky, warn::Bool = true) = Wishart(df, PDMat(S), warn)
 
 #  -----------------------------------------------------------------------------
 #  REPL display
@@ -66,23 +78,30 @@ show(io::IO, d::Wishart) = show_multline(io, d, [(:df, d.df), (:S, Matrix(d.S))]
 
 function convert(::Type{Wishart{T}}, d::Wishart) where T<:Real
     P = convert(AbstractArray{T}, d.S)
-    Wishart{T, typeof(P)}(T(d.df), P, T(d.logc0))
+    Wishart{T, typeof(P), typeof(d.rank)}(T(d.df), P, T(d.logc0), d.rank, d.singular)
 end
-function convert(::Type{Wishart{T}}, df, S::AbstractPDMat, logc0) where T<:Real
+function convert(::Type{Wishart{T}}, df, S::AbstractPDMat, logc0, rnk, singular) where T<:Real
     P = convert(AbstractArray{T}, S)
-    Wishart{T, typeof(P)}(T(df), P, T(logc0))
+    Wishart{T, typeof(P), typeof(rnk)}(T(df), P, T(logc0), rnk, singular)
 end
 
 #  -----------------------------------------------------------------------------
 #  Properties
 #  -----------------------------------------------------------------------------
 
-insupport(::Type{Wishart}, X::Matrix) = isposdef(X)
-insupport(d::Wishart, X::Matrix) = size(X) == size(d) && isposdef(X)
+insupport(::Type{Wishart}, X::AbstractMatrix) = ispossemdef(X)
+function insupport(d::Wishart, X::AbstractMatrix)
+    size(X) == size(d) || return false
+    if d.singular
+        return ispossemdef(X, rank(d))
+    else
+        return isposdef(X)
+    end
+end
 
 dim(d::Wishart) = dim(d.S)
 size(d::Wishart) = (p = dim(d); (p, p))
-rank(d::Wishart) = dim(d)
+rank(d::Wishart) = d.rank
 params(d::Wishart) = (d.df, d.S)
 @inline partype(d::Wishart{T}) where {T<:Real} = T
 
@@ -95,6 +114,7 @@ function mode(d::Wishart)
 end
 
 function meanlogdet(d::Wishart)
+    d.singular && return -Inf
     p = dim(d)
     df = d.df
     v = logdet(d.S) + p * logtwo
@@ -105,6 +125,7 @@ function meanlogdet(d::Wishart)
 end
 
 function entropy(d::Wishart)
+    d.singular && throw(ArgumentError("entropy not defined for singular Wishart."))
     p = dim(d)
     df = d.df
     -d.logc0 - 0.5 * (df - p - 1) * meanlogdet(d) + 0.5 * df * p
@@ -125,13 +146,44 @@ end
 #  Evaluation
 #  -----------------------------------------------------------------------------
 
-function wishart_logc0(df::Real, S::AbstractPDMat)
-    h_df = df / 2
+function wishart_logc0(df::Real, S::AbstractPDMat, rnk::Integer)
     p = dim(S)
-    -h_df * (logdet(S) + p * typeof(df)(logtwo)) - logmvgamma(p, h_df)
+    if df <= p - 1
+        return singular_wishart_logc0(p, df, S, rnk)
+    else
+        return nonsingular_wishart_logc0(p, df, S)
+    end
 end
 
 function logkernel(d::Wishart, X::AbstractMatrix)
+    if d.singular
+        return singular_wishart_logkernel(d, X)
+    else
+        return nonsingular_wishart_logkernel(d, X)
+    end
+end
+
+#  Singular Wishart pdf: Theorem 6 in Uhlig (1994 AoS)
+function singular_wishart_logc0(p::Integer, df::Real, S::AbstractPDMat, rnk::Integer)
+    h_df = df / 2
+    -h_df * (logdet(S) + p * typeof(df)(logtwo)) - logmvgamma(rnk, h_df) + (rnk*(rnk - p) / 2)*typeof(df)(logπ)
+end
+
+function singular_wishart_logkernel(d::Wishart, X::AbstractMatrix)
+    df = d.df
+    p = dim(d)
+    r = rank(d)
+    L = eigvals(Hermitian(X), (p - r + 1):p)
+    0.5 * ((df - (p + 1)) * sum(log.(L)) - tr(d.S \ X))
+end
+
+#  Nonsingular Wishart pdf
+function nonsingular_wishart_logc0(p::Integer, df::Real, S::AbstractPDMat)
+    h_df = df / 2
+    -h_df * (logdet(S) + p * typeof(df)(logtwo)) - logmvgamma(p, h_df)
+end
+
+function nonsingular_wishart_logkernel(d::Wishart, X::AbstractMatrix)
     df = d.df
     p = dim(d)
     Xcf = cholesky(X)
@@ -143,7 +195,12 @@ end
 #  -----------------------------------------------------------------------------
 
 function _rand!(rng::AbstractRNG, d::Wishart, A::AbstractMatrix)
-    _wishart_genA!(rng, dim(d), d.df, A)
+    if d.singular
+        A .= zero(eltype(A))
+        A[:, 1:rank(d)] = randn(rng, dim(d), rank(d))
+    else
+        _wishart_genA!(rng, dim(d), d.df, A)
+    end
     unwhiten!(d.S, A)
     A .= A * A'
 end
@@ -179,7 +236,7 @@ end
 
 function _rand_params(::Type{Wishart}, elty, n::Int, p::Int)
     n == p || throw(ArgumentError("dims must be equal for Wishart"))
-    ν = elty( n + 1 + abs(10randn()) )
+    ν = elty( n - 1 + abs(10randn()) )
     S = (X = 2rand(elty, n, n) .- 1; X * X')
     return ν, S
 end

src/utils.jl

@@ -64,3 +64,58 @@ function trycholesky(a::Matrix{Float64})
         return e
     end
 end
+
+"""
+    ispossemdef(A, k) -> Bool
+Test whether a matrix is positive semi-definite with specified rank `k` by
+checking that `k` of its eigenvalues are positive and the rest are zero.
+# Examples
+```jldoctest
+julia> A = [1 0; 0 0]
+2×2 Array{Int64,2}:
+ 1  0
+ 0  0
+julia> ispossemdef(A, 1)
+true
+julia> ispossemdef(A, 2)
+false
+```
+"""
+function ispossemdef(X::AbstractMatrix, k::Int;
+                     atol::Real=0.0,
+                     rtol::Real=(minimum(size(X))*eps(real(float(one(eltype(X))))))*iszero(atol))
+    _check_rank_range(k, minimum(size(X)))
+    ishermitian(X) || return false
+    dp, dz, dn = eigsigns(Hermitian(X), atol, rtol)
+    return dn == 0 && dp == k
+end
+function ispossemdef(X::AbstractMatrix;
+                     atol::Real=0.0,
+                     rtol::Real=(minimum(size(X))*eps(real(float(one(eltype(X))))))*iszero(atol))
+    ishermitian(X) || return false
+    dp, dz, dn = eigsigns(Hermitian(X), atol, rtol)
+    return dn == 0
+end
+
+function _check_rank_range(k::Int, n::Int)
+    0 <= k <= n || throw(ArgumentError("rank must be between 0 and $(n) (inclusive)"))
+    nothing
+end
+
+#  return counts of the number of positive, zero, and negative eigenvalues
+function eigsigns(X::AbstractMatrix,
+                  atol::Real=0.0,
+                  rtol::Real=(minimum(size(X))*eps(real(float(one(eltype(X))))))*iszero(atol))
+    eigs = eigvals(X)
+    eigsigns(eigs, atol, rtol)
+end
+function eigsigns(eigs::Vector{<: Real}, atol::Real, rtol::Real)
+    tol = max(atol, rtol * eigs[end])
+    eigsigns(eigs, tol)
+end
+function eigsigns(eigs::Vector{<: Real}, tol::Real)
+    dp = count(x -> tol < x, eigs)        #  number of positive eigenvalues
+    dz = count(x -> -tol < x < tol, eigs) #  number of numerically zero eigenvalues
+    dn = count(x -> x < -tol, eigs)       #  number of negative eigenvalues
+    return dp, dz, dn
+end

test/matrixvariates.jl

@@ -28,7 +28,6 @@ import Distributions: _univariate, _multivariate, _rand_params
 
 function test_draw(d::MatrixDistribution, X::AbstractMatrix)
     @test size(d) == size(X)
-    @test size(d) == size(mean(d))
     @test size(d, 1) == size(X, 1)
     @test size(d, 2) == size(X, 2)
     @test length(d) == length(X)
@@ -331,18 +330,32 @@ function test_special(dist::Type{Wishart})
         @test pvalue(kstest) >= α
     end
     @testset "H ~ W(ν, I) ⟹ H[i, i] ~ χ²(ν)" begin
-        ν = n + 1
-        ρ = Chisq(ν)
-        d = Wishart(ν, ScalMat(n, 1))
+        κ = n + 1
+        ρ = Chisq(κ)
+        g = Wishart(κ, ScalMat(n, 1))
         mymats = zeros(n, n, M)
         for m in 1:M
-            mymats[:, :, m] = rand(d)
+            mymats[:, :, m] = rand(g)
         end
         for i in 1:n
             kstest = ExactOneSampleKSTest(mymats[i, i, :], ρ)
             @test pvalue(kstest) >= α / n
         end
     end
+    @testset "Check Singular Branch" begin
+        X = H[1]
+        rank1 = Wishart(n - 2, Σ, false)
+        rank2 = Wishart(n - 1, Σ, false)
+        test_draw(rank1)
+        test_draw(rank2)
+        test_draws(rank1, rand(rank1, 10^6))
+        test_draws(rank2, rand(rank2, 10^6))
+        test_cov(rank1)
+        test_cov(rank2)
+        @test Distributions.singular_wishart_logkernel(d, X) ≈ Distributions.nonsingular_wishart_logkernel(d, X)
+        @test Distributions.singular_wishart_logc0(n, ν, d.S, rank(d)) ≈ Distributions.nonsingular_wishart_logc0(n, ν, d.S)
+        @test logpdf(d, X) ≈ Distributions.singular_wishart_logkernel(d, X) + Distributions.singular_wishart_logc0(n, ν, d.S, rank(d))
+    end
     nothing
 end
 

test/utils.jl

@@ -1,6 +1,6 @@
 using Distributions, PDMats
 using Test, LinearAlgebra
-
+import Distributions: ispossemdef
 
 # RealInterval
 r = RealInterval(1.5, 4.0)
@@ -36,3 +36,33 @@ N = GenericArray([1.0 0.0; 1.0 0.0])
 
 @test Distributions.isApproxSymmmetric(N) == false
 @test Distributions.isApproxSymmmetric(M)
+
+
+n = 10
+areal = randn(n,n)/2
+aimg  = randn(n,n)/2
+@testset "For A containing $eltya" for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int)
+    ainit = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
+    @testset "Positive semi-definiteness" begin
+        notsymmetric = ainit
+        notsquare    = [ainit ainit]
+        @test !ispossemdef(notsymmetric)
+        @test !ispossemdef(notsquare)
+        for truerank in 0:n
+            X = ainit[:, 1:truerank]
+            A = truerank == 0 ? zeros(eltya, n, n) : X * X'
+            @test ispossemdef(A)
+            for testrank in 0:n
+                if testrank == truerank
+                    @test ispossemdef(A, testrank)
+                else
+                    @test !ispossemdef(A, testrank)
+                end
+            end
+            @test !ispossemdef(notsymmetric, truerank)
+            @test !ispossemdef(notsquare, truerank)
+            @test_throws ArgumentError ispossemdef(A, -1)
+            @test_throws ArgumentError ispossemdef(A, n + 1)
+        end
+    end
+end