JSoC'19: Trace of Matrix Inverse

Introduction

This article will focus on estimation of trace of large sparse matrix inverse, that is, where is a Symmetric Positive Definite matrix andThis method is based on Estimating the Trace of the Matrix Inverse by Interpolating from the Diagonal of an Approximate Inverse by Lingfei Wu, et al.

Diagonal Approximation Algorithm

Diagonal Approximation Algorithm makes use of a preconditioner to generate an initial guess of the distribution pattern for the diagonal elements. This initial guess is used to find the a set of points that best describe the characteristics of the diagonal. With help of these points, we interpolate the diagonal using a linear or a cubic hermite polynomial model. These steps can be summarized in the following high-level description of the algorithm:

Findas an approximate inverse tousing some preconditioner such as IncompleteLU or Incomplete Cholesky (or any other Preconditioner). Take
Compute fitting sample a set of indices that best describes the distribution pattern for diagonal of
Solve linear system
Obtain a fitting model using linear or polynomial regression methods.
Use the fitting model to calculate diagonal elements and trace of the matrix inverse.

Usage

Setup the environment.

# Ingredients
import Pkg;
Pkg.add(PackageSpec(url="https://github.com/luca-aki/traceEstimation.jl", rev = "dapp"))
Pkg.add("MatrixMarket")
using LinearAlgebra, SparseArrays, MatrixMarket, TraceEstimation

80.7s

Julia

Let us import a few matrices from SuiteSparse to work on.

Nasa2146 Matrix (Structural Problem)

Size: 2146x2146
Nonzero Entries: 72,250
Condition Number: 1.724336e+03
Trace of Matrix Inverse: 0.003418

nasa2146.mtx

A = MatrixMarket.mmread(nasa2146.mtx
)

# We will use tr_inv from the diagonalapprox.jl in TraceEstimation package

# tr_inv(A::SPD-Matrix, Probing-Vector-Count, Sample-Point-Count; Preconditioner, Fitting-Model)
# Here, we are not passing values for Preconditioner or Fitting model. Default Incomplete Cholesky and Linear Regression will be used.
tr_inv(A, 6, 40)

0.7s

Julia

0.00348673

Kuu Matrix (Structural Problem)

Size: 7102x7102

Nonzero Entries: 340,200
Condition Number: 1.575800e+04
Trace of Matrix Inverse: 3618.675849

Kuu.mtx

A = MatrixMarket.mmread(Kuu.mtx
)

# Here, we are passing IncompleteLU and pchip as Preconditioner and fitting model resp.
tr_inv(A, 8, 80, pc = "ilu", model = "pchip")

4.9s

Julia

3562.28

Trefethen_20000 (Combinatorial Problem)

Size:20,000x20,000
Nonzero Entries: 554,466
Condition Number: 2.005593e+05
Trace of Matrix Inverse: 3.2153278

Trefethen_20000.mtx

A = MatrixMarket.mmread(Trefethen_20000.mtx
)
A = SparseMatrixCSC{Float64, Int64}(A)
# We will increase Sample Point count as the matrix is large.
tr_inv(A, 8, 250, pc = "ilu")

3.5s

Julia

3.2471

Scope for Improvement and Insights

As evident from the Step 1 and Step 4 of the algorithm given above, we can plug different preconditioners and fitting models to obtain a better approximation of diagonal elements and ultimately the trace of the matrix inverse. These can be modified according to the type of problem one is trying to solve.

It is also evident that some methods would be faster than others. For example, calculating IncompleteLU is faster than calculating Incomplete Cholesky, but the same can be less accurate sometimes. Step 3 of the algorithm uses CG to solve the linear systems, which also increases the execution time as the number of sample points increase.

Thank you for reading!