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.
Let us import a few matrices from SuiteSparse to work on.
Nasa2146 Matrix (Structural Problem)
Nonzero Entries: 72,250
Condition Number: 1.724336e+03
Trace of Matrix Inverse: 0.003418
# 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)
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.