SparseMvNormal Runtime

using Pkg; Pkg.add(PackageSpec(url="https://github.com/luca-aki/traceEstimation.jl"))
Pkg.add("Distributions")
Pkg.add("IterativeSolvers")
Pkg.add("JLD2")
using LinearAlgebra, SparseArrays, Distributions, IterativeSolvers

80.0s

SparseMvNormal (Julia)

Following is an example application of Chebyhutch log-determinant. In this example, we are creating a struct SparseMvNormal that inherits from AbstractMvNormal. We will be using this to calculate log probability density function.

The implementation is fairly straightforward. The major change that we have done here is to include chebyhutch as a replacement for logdet.

export SparseMvNormal

struct SparseMvNormal{T<:Real, Cov<:AbstractSparseMatrix, Mean<:AbstractVector{T}} <: AbstractMvNormal
    μ::Mean
    Σ::Cov
end

# Constructor
function SparseMvNormal(μ::AbstractVector{T}, Σ::AbstractSparseMatrix{T}) where {T<:Real}
    size(Σ, 1) == size(μ, 1) || throw(DimensionMismatch("The dimensions of mu and Sigma are inconsistent."))
    SparseMvNormal{T, typeof(Σ), typeof(μ)}(μ, Σ)
end

# Basic stats
Base.length(d::SparseMvNormal) = length(d.μ)
Distributions.mean(d::SparseMvNormal) = d.μ

Distributions.var(d::SparseMvNormal) = diag(d.Σ)
Distributions.cov(d::SparseMvNormal) = d.Σ

struct Inv{T, Cov<:AbstractSparseMatrix}
    cov::Cov
end
Base.:\(a::Inv, x::AbstractVector) = cg(a.cov, x)
Distributions.invcov(d::SparseMvNormal) = Inv(d.cov)
Distributions.logdetcov(d::SparseMvNormal) = chebyhutch(d.Σ, 4, 6, fn = log)

# Evaluation
invquad(a::AbstractSparseMatrix, x::StridedVector) = dot(x, cg(a, x))
Distributions.sqmahal(d::SparseMvNormal, x::AbstractVector) = invquad(d.Σ, x .- d.μ)
Distributions.sqmahal!(r::AbstractVector, d::SparseMvNormal, x::AbstractMatrix) = invquad!(r, d.Σ, x .- d.μ)
gradlogpdf(d::SparseMvNormal, x::AbstractVector) = -(d.Σ \ (x .- d.μ))

SparseMvNormalCode