Approximating Spectral Sums of Large-Scale Matrices using Stochastic Chebyshev Approximations

Abstract

Computation of the trace of a matrix function plays an important role in many scientific computing applications, including applications in machine learning, computational physics (e.g., lattice quantum chromodynamics), network analysis and computational biology (e.g., protein folding), just to name a few application areas. We propose a linear-time randomized algorithm for approximating the trace of matrix functions of large symmetric matrices. Our algorithm is based on coupling function approximation using Chebyshev interpolation with stochastic trace estimators (Hutchinson's method), and as such requires only implicit access to the matrix, in the form of a function that maps a vector to the product of the matrix and the vector. We provide rigorous approximation error in terms of the extremal eigenvalue of the input matrix, and the Bernstein ellipse that corresponds to the function at hand. Based on our general scheme, we provide algorithms with provable guarantees for important matrix computations, including log-determinant, trace of matrix inverse, Estrada index, Schatten p-norm, and testing positive definiteness. We experimentally evaluate our algorithm and demonstrate its effectiveness on matrices with tens of millions dimensions. * This article is partially based on preliminary results published in the proceeding of the 32nd International Conference on Machine Learning (ICML 2015).

Genome-wide association analysis of age-at-onset traits using Cox mixed-effects models

“…The nodes and weights in the Gaussian quadrature can be elegantly computed by the eigenvalues and first elements in the eigenvectors of the tridiagonal matrix obtained from the Lanczos algorithm. We found in our simulation study that the accuracy of the SLQ method for approximating a log-determinant was generally higher by one order of magnitude than that using Chebyshev orthogonal polynomials (Han et al, 2016;Pace and LeSage, 2004) , both of which were more accurate than Martin's Taylor expansion (Barry and Kelley Pace, 1999;Martin, 1992) , consistent with previous reports (Han et al, 2016;Ubaru et al, 2017) . More details about the SLQ approximation were described in Supplementary materials A.4.…”

Spectrum-adapted Polynomial Approximation for Matrix Functions

“…Efficiently computing f (A)b, a function of a large, sparse Hermitian matrix times a vector, is an important component in numerous signal processing, machine learning, applied mathematics, and computer science tasks. Application examples include graph-based semi-supervised learning methods [2]- [4]; graph spectral filtering in graph signal processing [5]; convolutional neural networks / deep learning [6,7]; clustering [8,9]; approximating the spectral density of a large matrix [10]; estimating the numerical rank of a matrix [11,12]; approximating spectral sums such as the log-determinant of a matrix [13] or the trace of a matrix inverse for applications in physics, biology, information theory, and other disciplines [14]; solving semidefinite programs [15]; simulating random walks [16,Chapter 8]; and solving ordinary and partial differential equations [17]- [19].…”

A posteriori error estimate for computing tr(f(A)) by using the Lanczos method

“…Computing the trace in this manner, however, requires the computation of all the eigenvalues, which is also often prohibitively expensive. Hence, various methods proposed for approximately computing tr( f ( A )) consist of the following two ingredients: Approximate the trace of f ( A ) by using the average of unbiased samples , i =1,…, N , where the u i are independent random vectors of some nature. Approximately compute the bilinear form by using some numerical technique. …”

“…Computing the trace in this manner, however, requires the computation of all the eigenvalues, which is also often prohibitively expensive. Hence, various methods proposed for approximately computing tr( f(A)) consist of the following two ingredients 1,10,[13][14][15][16][17][18][19] The various methods differ in the random mechanism of selecting the u i and the numerical technique for computing the bilinear form. Several variants of these ingredients exist (e.g., computing deterministically tr( (A)) = ∑ n i=1 e T i (A)e i rather than using random vectors u i , or even using block vectors to replace the canonical vectors e i 20 ; or using moment extrapolation for, particularly, f(t) = t with real value 21,22 ), but they are not the focus of this work.…”

“…In the two ingredients, the convergence of the approximation is generally gauged through a combination of Monte Carlo convergence of the sample average and the convergence of the numerical technique. In practice, however, the convergence results obtained for these methods 10,15,17,23,24 rarely translate into practical schemes to monitor convergence. Monitoring the accuracy of a given approximation to tr( f(A)) can be very challenging for certain functions f. This is in complete contrast to the situation prevalent when solving linear systems, where simple residual norms provide a computable measure of the backward error.…”