JSoC'19: Logdet via Chebyhutch

Log-determinant

Logarithm determinant, aka log-det, is useful in numerous statistical applications and related fields, including machine learning, for example, it is used in calculating log probability density function (logpdf) and maximum log-likelihood function.

The log-determinant of a positive definite matrix can be calculated by finding the trace of

logdet(A) = tr(log(A))

We have previously seen ways to estimate trace of logarithm by using Stochastic Lanczos Quadrature and Chebyshev-Hutchinson Polynomial Interpolation (or Chebyhutch). In this article, we are going to focus on using Chebyhutch to calculate logdet for sparse matrices.

Chebyhutch for log-determinant

(Article on Chebyhutch)

Here is a matrix that I grabbed from SuiteSparse Matrix Collection website. We will be using this for our testing purpose.

Kuu.jld2

Kuu matrix is a Symmetric Positive Definite matrix of size 7102, nonzero count of 340,200 and condition number 1.575800e+04.

using LinearAlgebra, SparseArrays, BenchmarkTools, TraceEstimation, FileIO, JLD2, Distributions, IterativeSolvers

11.3s

Julia

@load Kuu.jld2
 A

9.3s

Julia

1-element Array{Symbol,1}: :A

That's all we need to start our analysis. Let's first take a look at the available method logdet in the LinearAlgebra library.

M = Matrix(A)
@benchmark logdet(M)

37.5s

Julia

BenchmarkTools.Trial: memory estimate: 384.87 MiB allocs estimate: 6 -------------- minimum time: 9.142 s (0.00% GC) median time: 9.142 s (0.00% GC) mean time: 9.142 s (0.00% GC) maximum time: 9.142 s (0.00% GC) -------------- samples: 1 evals/sample: 1

Now, we will try running Chebyhutch on the same matrix for log-determinant. Chebyhutch requires 4 parameters:

A Symmetric Positive Definite Matrix.
Iteration number (higher iteration number improves trace precision).
Polynomial Degree (higher degree improves trace accuracy).
An analytical function such as log or sqrt.

@benchmark chebyhutch(A, 4, 6, fn = log)

12.0s

Julia

BenchmarkTools.Trial: memory estimate: 6.84 MiB allocs estimate: 342 -------------- minimum time: 82.443 ms (0.00% GC) median time: 97.026 ms (0.00% GC) mean time: 98.233 ms (0.40% GC) maximum time: 129.210 ms (1.05% GC) -------------- samples: 51 evals/sample: 1

So, it is evident that Chebyhutch can provide a significant boost in the calculation of sparse log-determinant but what about the accuracy? Here is a simple relative error analysis.

act = logdet(M)
obv = chebyhutch(A, 4, 6, fn = log)
abs((act - obv)/act)

10.4s

Julia

0.00225348

A relative error in the magnitude of 10is good enough for most of the numerical work.

Application - SparseMvNormal

I have already setup the environment for SparseMvNormal, which is internally using Chebyhutch for calculating the logdet. We are going to see a simple application where we find logpdf of a sparse covariance matrix with a random mean vector.

# Start by importing the code from environment
include("/SparseMvNormalCode")

# Generating the vectors for MvNormal example
v = rand(size(A, 1))
x = rand(-100.0:2:100.0, size(A, 1))
println("")

0.9s

Julia

Using the above generated v and x, let us first calculate logpdf for MvNormal.

d = MvNormal(v, M)
logpdf(d, x)

6.6s

Julia

-5.70833e6

Now, let us do the same for SparseMvNormal using the original sparse matrix

include("/SparseMvNormalCode")
ds = SparseMvNormal(v, A)
logpdf(ds, x)

2.0s

Julia

-5.7083e6

As we can see, the results are good for any numerical purposes.

Since MvNormal uses Cholesky decomposition to calculate logdet we should to do another execution time analysis.

# Wrapper for complete logpdf using MvNormal
function clogpdf(A::AbstractMatrix, v::AbstractVector, x::AbstractVector)
  d = MvNormal(v, A)
  logpdf(d, x)
end
# Wrapper for complete logpdf using SparseMvNormal
function clogpdf(A::AbstractSparseMatrix, v::AbstractVector, x::AbstractVector)
  d = SparseMvNormal(v, A)
  logpdf(d, x)
end

0.7s

Julia

clogpdf (generic function with 2 methods)

@benchmark clogpdf($M, v, x)

21.7s

Julia

BenchmarkTools.Trial: memory estimate: 384.92 MiB allocs estimate: 13 -------------- minimum time: 5.112 s (0.01% GC) median time: 5.112 s (0.01% GC) mean time: 5.112 s (0.01% GC) maximum time: 5.112 s (0.01% GC) -------------- samples: 1 evals/sample: 1

@benchmark clogpdf($A, v, x)

15.5s

Julia

BenchmarkTools.Trial: memory estimate: 7.17 MiB allocs estimate: 367 -------------- minimum time: 727.398 ms (0.17% GC) median time: 733.421 ms (0.00% GC) mean time: 734.705 ms (0.05% GC) maximum time: 746.142 ms (0.00% GC) -------------- samples: 7 evals/sample: 1

With this, I will end the SparseMvNormal example.

Thank you for reading!