# JSoC'19: Stochastic Trace Estimation

## Hutchinson's estimator

Hutchinson in 1989, came up with a general technique to transform algebraic problem of calculating trace into statistical problem of computing expectation of a quadratic function.

Here,

### How to solve the given bilinear form?

Let *what method to employ while evaluating*A quick resolution is using some iterative solver such as CG to evaluate

## Chebyshev Polynomial Interpolation

When working with Stochastic Lanczos Quadrature method, we built a *tridiagonal* matrix, to extract Ritz-values (a close approximation to eigenvalues of the original matrix) and used the identity

where,

Chebyshev polynomial of first kind defined recursively as

Where,

### Cheby-Hutch Trace Estimator

With the interpolation power of chebyshev interpolation, we can feed the values of

- Find
and the extreme eigenvalues of the - Calculate the values for interpolation coefficient
- Calculate
to approximate for the eigenvalues of - Form
and feed it to the Hutchinson estimator for calculation of

## Usage and Comparison

### Code Snippet

Add the required modules

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

Create dummy random matrices dense and sparse

A = rand(5000,5000) # Make sure it is symmetric positive definite A = A + A' + 250I while(!isposdef(A)) A = A + 100I end B = sprand(5000,5000,0.001) B = B + B' + 250I while(!isposdef(B)) B = B + 100I end

Now we benchmark Cheby-Hutch comparing it with default methods

tr(inv(A))

chebyhutch(A, 4, 6, fn = inv)

For sparse matrices, the performance increases by manyfold.

Bm = Matrix(B) tr(inv(Bm))

chebyhutch(B, 4, 6, fn = inv)

### Because Everyone Loves Graphs!

Following tests are performed on random dense matrices generated inside julia.

Test details:

- Default: tr(inv(A))
- SLQ: slq(A, fn = inv, mtol = 0.005)
- Cheby-Hutch: chebyhutch(A, 6, 25, fn = inv)

Test details:

- Default: tr(sqrt(A))
- SLQ: slq(A, fn = sqrt, mtol = 0.005)
- Cheby-Hutch: chebyhutch(A, 6, 25, fn = sqrt)

Test details:

- SLQ: slq(A, fn = inv, mtol = 0.005)
- Cheby-Hutch: chebyhutch(A, 6, 25, fn = inv)

Following tests are performed on a Sparse Matrix downloaded form SuiteSparse Matrix Collection and imported using MatrixMarket.jl.

Test details:

Matrix: Kuu.mtx (7102 x 7102), Condition number = 1.57e+04

- SLQ: slq(A, mtol = X, eps = X, fn = inv)
- Cheby-Hutch: chebyhutch(A, 15, X, fn = inv)

## Up Next

Continuing the journey with trace estimation method, I will be working on Diagonal Approximation method. Along with that, I will also be working on a Sparse Approximate Inverse to serve as preconditioner for Diagonal Approximation.

Thank you for reading!