Hierarchical Probing for Estimating the Trace of the Matrix Inverse on Toroidal Lattices

The standard approach for computing the trace of the inverse of a very large, sparse matrix A is to view the trace as the mean value of matrix quadratures, and use the Monte Carlo algorithm to estimate it. This approach is heavily used in our motivating application of Lattice QCD. Often, the elements of A−1 display certain decay properties away from the non zero structure of A, but random vectors cannot exploit this induced structure of A−1. Probing is a technique that, given a sparsity pattern of A, discovers elements of A through matrix-vector multiplications with specially designed vectors. In the case of A−1, the pattern is obtained by distance-k coloring of the graph of A. For sufficiently large k, the method produces accurate trace estimates but the cost of producing the colorings becomes prohibitively expensive. More importantly, it is difficult to search for an optimal k value, since none of the work for prior choices of k can be reused. First, we introduce the idea of hierarchical probing that produces distance-2i colorings for a sequence of distances 20, 21, . . . , 2m, up to the diameter of the graph. To achieve this, we do not color the entire graph, but at each level, i, we compute the distance-1 coloring independently for each of the node-groups associated with a color of the distance-(2i−1) coloring. Second, based on this idea, we develop an algorithm for uniform, toroidal lattices that simply applies bit-arithmetic on local coordinates to produce the hierarchical permutation. Third, we provide an algorithm for choosing an appropriate sequence of Hadamard and Fourier vectors, so that earlier vectors in the sequence correspond to hierarchical probing vectors of smaller distances. This allows us to increase the number of probing vectors until the required accuracy is achieved. Several experiments show that when a decay structure exists in the matrix, our algorithm finds it and approximates the trace incrementally, starting with the most important contributions. We have observed up to an order of magnitude speedup over the standard Monte Carlo.