# 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. $$Tr(A) = \mathbb{E}(\textbf{z}^T A\textbf{z})$$ Here, $\textbf{z}$is a random variable usually following either Gaussian or Rademacher distribution, i.e., $z_i$are independent vectors of zero-mean, unit-variance and uncorrelated entries. ## How to solve the given bilinear form? Let $M = A^{-1},$general case being $M = f(A),$then $tr(M) = \mathbb{E}(\textbf{z}^T M\textbf{z}) =$ $\frac{1}{N}\sum_{i=1}^N(z_i^TMz_i).$One prominent question that arises from the use of this technique is *what method to employ while evaluating*$Mz_i = A^{-1}z_i.$A quick resolution is using some iterative solver such as CG to evaluate$Ax_i = z_i.$We can have a good estimate of $tr(M)$via this procedure in most cases, but if$A$is so ill-conditioned that even a stable iterative solver leaves a substantial backward error in evaluation of$Mz_i$then we will have a significant departure from true value of $tr(A).$While such a scenario does occur in practice, it is fairly uncommon. On the other hand, the opposite scenario where it is possible to have a good estimates for$Mz_i,$given$A$is not highly ill-conditioned for our solver. It is also possible to estimate these values via polynomial interpolation. # 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 $tr(f(A)) = \sum(f(\lambda_i))$for trace estimation. Now, we take a different but fundamentally similar approach to the problem by interpolating eigenvalues of the matrix for a given function $f:[a, b] \rightarrow\mathbb{R}$via chebyshev interpolation. $$f(x) \equiv p_n(x) = \sum_{j=0}^n(c_jT_j(\frac{2}{b-a}x - \frac{b+a}{b-a}))$$ where, $p_n$is the interpolation polynomial of degree $n,$$c_jis the interpolation coefficient and T_j(x)is the j^{th}chebyshev polynomial. Chebyshev polynomial of first kind defined recursively as$$ T_o(x) = 1,~ T_1(x) = x,  T_{j+1}(x) = 2xT_j(x) - T_{j-1}(x)$$Where,$x_k = cos(\frac{\pi(k+0.5)}{n+1}).$## Cheby-Hutch Trace Estimator With the interpolation power of chebyshev interpolation, we can feed the values of$Mz_i$into out traditional Hutchinson estimator. A high-level description of the algorithm for a symmetric positive definite matrix$A\in\mathbb{R}^{n\times n}$is given below: * Find$\lambda_1$and$\lambda_n,$the extreme eigenvalues of the$A.$* Calculate the values for interpolation coefficient$c_j.$* Calculate$p_n$to approximate$f$for the eigenvalues of$A.$* Form$Mz_i,$and feed it to the Hutchinson estimator for calculation of$tr(f(A)).\$ # Usage and Comparison ## Code Snippet Add the required modules julia id=a66a486d-f0d7-4f68-8be0-e0da6a92cf18 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 julia id=07a7d88b-2616-42fd-a098-820c611db8c2 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 julia id=692637e8-45b5-444b-86aa-6f70da609cdb @benchmark tr(inv(A))  julia id=a34b09bf-015e-40a6-989f-0a395b2d20c2 @benchmark chebyhutch(A, 4, 6, fn = inv)  For sparse matrices, the performance increases by manyfold. julia id=090e1be7-e583-465a-8170-74babfb415fd Bm = Matrix(B) @benchmark tr(inv(Bm))  julia id=250caf7c-14b2-4150-bee8-f3ae764bf62b @benchmark 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) ![speedtrinv2.png][nextjournal#file#963eae0e-770f-49a0-a629-9cc85c99a158] Test details: * Default: tr(sqrt(A)) * SLQ: slq(A, fn = sqrt, mtol = 0.005) * Cheby-Hutch: chebyhutch(A, 6, 25, fn = sqrt) ![speedtrsqrt.png][nextjournal#file#a9f29e63-989f-4bec-91c9-d7790addb598] Test details: * SLQ: slq(A, fn = inv, mtol = 0.005) * Cheby-Hutch: chebyhutch(A, 6, 25, fn = inv) ![error.PNG][nextjournal#file#758d998e-3806-4600-9c1b-0818baa9e1b1] 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) ![error3.PNG][nextjournal#file#59b65b3c-7d6b-4393-b5be-61145da88db5] ![error2.PNG][nextjournal#file#4fd99d40-c470-4147-8ac3-6cdc8d8c3cf7] # 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! [nextjournal#file#963eae0e-770f-49a0-a629-9cc85c99a158]: [nextjournal#file#a9f29e63-989f-4bec-91c9-d7790addb598]: [nextjournal#file#758d998e-3806-4600-9c1b-0818baa9e1b1]: [nextjournal#file#59b65b3c-7d6b-4393-b5be-61145da88db5]: [nextjournal#file#4fd99d40-c470-4147-8ac3-6cdc8d8c3cf7]:
This notebook was exported from https://nextjournal.com/a/LKaM6wyaUsmM5ap6GW2HG?change-id=CVeXQdaqk7AiyncQQSTEiy edn nextjournal-metadata {:article {:settings nil, :nodes {"07a7d88b-2616-42fd-a098-820c611db8c2" {:compute-ref #uuid "a6d6f9d5-ab71-4b59-a32c-3e60fdf553c7", :exec-duration 2318, :id "07a7d88b-2616-42fd-a098-820c611db8c2", :kind "code", :output-log-lines {:stdout 0}, :runtime [:runtime "46ea9893-c4ef-4ed5-b145-e68e1fc45d6a"]}, "090e1be7-e583-465a-8170-74babfb415fd" {:compute-ref #uuid "fa28ed18-05bd-4601-b477-ccad6623b42f", :exec-duration 32408, :id "090e1be7-e583-465a-8170-74babfb415fd", :kind "code", :output-log-lines {:stdout 0}, :runtime [:runtime "46ea9893-c4ef-4ed5-b145-e68e1fc45d6a"]}, "250caf7c-14b2-4150-bee8-f3ae764bf62b" {:compute-ref #uuid "00a8d56c-9c0f-4400-bc0a-3e81d68b8be9", :exec-duration 12408, :id "250caf7c-14b2-4150-bee8-f3ae764bf62b", :kind "code", :output-log-lines {:stdout 0}, :runtime [:runtime "46ea9893-c4ef-4ed5-b145-e68e1fc45d6a"]}, "46ea9893-c4ef-4ed5-b145-e68e1fc45d6a" {:environment [:environment {:article/nextjournal.id #uuid "5b460d39-8c57-43a6-8b13-e217642b0146", :change/nextjournal.id #uuid "5cc9b90b-0bd4-4e38-bb95-834151b2dc86", :node/id "39e3f06d-60bf-4003-ae1a-62e835085aef"}], :id "46ea9893-c4ef-4ed5-b145-e68e1fc45d6a", :kind "runtime", :language "julia", :resources {:machine-type "n1-standard-4"}, :type :nextjournal}, "4fd99d40-c470-4147-8ac3-6cdc8d8c3cf7" {:id "4fd99d40-c470-4147-8ac3-6cdc8d8c3cf7", :kind "file"}, "59b65b3c-7d6b-4393-b5be-61145da88db5" {:id "59b65b3c-7d6b-4393-b5be-61145da88db5", :kind "file"}, "692637e8-45b5-444b-86aa-6f70da609cdb" {:compute-ref #uuid "993a60a8-ebf3-48b4-b55d-6616150fa8a1", :exec-duration 30193, :id "692637e8-45b5-444b-86aa-6f70da609cdb", :kind "code", :output-log-lines {:stdout 0}, :runtime [:runtime "46ea9893-c4ef-4ed5-b145-e68e1fc45d6a"]}, "758d998e-3806-4600-9c1b-0818baa9e1b1" {:id "758d998e-3806-4600-9c1b-0818baa9e1b1", :kind "file"}, "963eae0e-770f-49a0-a629-9cc85c99a158" {:id "963eae0e-770f-49a0-a629-9cc85c99a158", :kind "file"}, "a34b09bf-015e-40a6-989f-0a395b2d20c2" {:compute-ref #uuid "21319384-4a37-4ef0-8dab-08c7114074c3", :exec-duration 15203, :id "a34b09bf-015e-40a6-989f-0a395b2d20c2", :kind "code", :output-log-lines {:stdout 0}, :runtime [:runtime "46ea9893-c4ef-4ed5-b145-e68e1fc45d6a"]}, "a66a486d-f0d7-4f68-8be0-e0da6a92cf18" {:compute-ref #uuid "a2980825-36ae-4685-86ee-fb6b548496b5", :exec-duration 1255, :id "a66a486d-f0d7-4f68-8be0-e0da6a92cf18", :kind "code", :output-log-lines {:stdout 6}, :runtime [:runtime "46ea9893-c4ef-4ed5-b145-e68e1fc45d6a"]}, "a9f29e63-989f-4bec-91c9-d7790addb598" {:id "a9f29e63-989f-4bec-91c9-d7790addb598", :kind "file"}}, :nextjournal/id #uuid "02b2942f-23e5-4000-8679-7cd9bf9bf3ab", :article/change {:nextjournal/id #uuid "5d1479d9-32d4-4e03-8934-3165a034fc96"}}}