Akshay Jain / Aug 02 2019
Remix of Julia by Nextjournal

JSoC'19: Trace of Matrix Inverse


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:

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


Setup the environment.

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

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
A = MatrixMarket.mmread(
) # 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)

Kuu Matrix (Structural Problem)

Size: 7102x7102

  • Nonzero Entries: 340,200
  • Condition Number: 1.575800e+04
  • Trace of Matrix Inverse: 3618.675849
A = MatrixMarket.mmread(
) # Here, we are passing IncompleteLU and pchip as Preconditioner and fitting model resp. tr_inv(A, 8, 80, pc = "ilu", model = "pchip")

Trefethen_20000 (Combinatorial Problem)

  • Size:20,000x20,000
  • Nonzero Entries: 554,466
  • Condition Number: 2.005593e+05
  • Trace of Matrix Inverse: 3.2153278
A = MatrixMarket.mmread(
) A = SparseMatrixCSC{Float64, Int64}(A) # We will increase Sample Point count as the matrix is large. tr_inv(A, 8, 250, pc = "ilu")

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!