Krylov-aware stochastic trace estimation

We introduce an algorithm for estimating the trace of a matrix function $f(\mathbf{A})$ using implicit products with a symmetric matrix $\mathbf{A}$. Existing methods for implicit trace estimation of a matrix function tend to treat matrix-vector products with $f(\mathbf{A})$ as a black-box to be computed by a Krylov subspace method. Like other recent algorithms for implicit trace estimation, our approach is based on a combination of deflation and stochastic trace estimation. However, we take a closer look at how products with $f(\mathbf{A})$ are integrated into these approaches which enables several efficiencies not present in previously studied methods. In particular, we describe a Krylov subspace method for computing a low-rank approximation of a matrix function by a computationally efficient projection onto Krylov subspace.Comment: Figure 5.1 differs somewhat from the published version due to a clerical error made when uploading the images to the journa

Scalable iterative methods for sampling from massive Gaussian random vectors

Sampling from Gaussian Markov random fields (GMRFs), that is multivariate Gaussian ran- dom vectors that are parameterised by the inverse of their covariance matrix, is a fundamental problem in computational statistics. In this paper, we show how we can exploit arbitrarily accu- rate approximations to a GMRF to speed up Krylov subspace sampling methods. We also show that these methods can be used when computing the normalising constant of a large multivariate Gaussian distribution, which is needed for both any likelihood-based inference method. The method we derive is also applicable to other structured Gaussian random vectors and, in particu- lar, we show that when the precision matrix is a perturbation of a (block) circulant matrix, it is still possible to derive O(n log n) sampling schemes.Comment: 17 Pages, 4 Figure

Suboptimal subspace construction for log-determinant approximation

Variance reduction is a crucial idea for Monte Carlo simulation and the stochastic Lanczos quadrature method is a dedicated method to approximate the trace of a matrix function. Inspired by their advantages, we combine these two techniques to approximate the log-determinant of large-scale symmetric positive definite matrices. Key questions to be answered for such a method are how to construct or choose an appropriate projection subspace and derive guaranteed theoretical analysis. This paper applies some probabilistic approaches including the projection-cost-preserving sketch and matrix concentration inequalities to construct a suboptimal subspace. Furthermore, we provide some insights on choosing design parameters in the underlying algorithm by deriving corresponding approximation error and probabilistic error estimations. Numerical experiments demonstrate our method's effectiveness and illustrate the quality of the derived error bounds

Faster randomized partial trace estimation

We develop randomized matrix-free algorithms for estimating partial traces. Our algorithm improves on the typicality-based approach used in [T. Chen and Y-C. Cheng, Numerical computation of the equilibrium-reduced density matrix for strongly coupled open quantum systems, J. Chem. Phys. 157, 064106 (2022)] by deflating important subspaces (e.g. corresponding to the low-energy eigenstates) explicitly. This results in a significant variance reduction for matrices with quickly decaying singular values. We then apply our algorithm to study the thermodynamics of several Heisenberg spin systems, particularly the entanglement spectrum and ergotropy

Randomized matrix-free quadrature for spectrum and spectral sum approximation

We study randomized matrix-free quadrature algorithms for spectrum and spectral sum approximation. The algorithms studied are characterized by the use of a Krylov subspace method to approximate independent and identically distributed samples of $\mathbf{v}^{\sf H}f[\mathbf{A}]\mathbf{v}$, where $\mathbf{v}$ is an isotropic random vector, $\mathbf{A}$ is a Hermitian matrix, and $f[\mathbf{A}]$ is a matrix function. This class of algorithms includes the kernel polynomial method and stochastic Lanczos quadrature, two widely used methods for approximating spectra and spectral sums. Our analysis, discussion, and numerical examples provide a unified framework for understanding randomized matrix-free quadrature and shed light on the commonalities and tradeoffs between them. Moreover, this framework provides new insights into the practical implementation and use of these algorithms, particularly with regards to parameter selection in the kernel polynomial method

Randomized low-rank approximation of monotone matrix functions

This work is concerned with computing low-rank approximations of a matrix function $f(A)$ for a large symmetric positive semi-definite matrix $A$, a task that arises in, e.g., statistical learning and inverse problems. The application of popular randomized methods, such as the randomized singular value decomposition or the Nystr\"om approximation, to $f(A)$ requires multiplying $f(A)$ with a few random vectors. A significant disadvantage of such an approach, matrix-vector products with $f(A)$ are considerably more expensive than matrix-vector products with $A$, even when carried out only approximately via, e.g., the Lanczos method. In this work, we present and analyze funNystr\"om, a simple and inexpensive method that constructs a low-rank approximation of $f(A)$ directly from a Nystr\"om approximation of $A$, completely bypassing the need for matrix-vector products with $f(A)$. It is sensible to use funNystr\"om whenever $f$ is monotone and satisfies $f(0) = 0$. Under the stronger assumption that $f$ is operator monotone, which includes the matrix square root $A^{1/2}$ and the matrix logarithm $\log(I+A)$, we derive probabilistic bounds for the error in the Frobenius, nuclear, and operator norms. These bounds confirm the numerical observation that funNystr\"om tends to return an approximation that compares well with the best low-rank approximation of $f(A)$. Our method is also of interest when estimating quantities associated with $f(A)$, such as the trace or the diagonal entries of $f(A)$. In particular, we propose and analyze funNystr\"om++, a combination of funNystr\"om with the recently developed Hutch++ method for trace estimation

Accuracy of the finite-temperature Lanczos method compared to simple typicality-based estimates

We study trace estimators for equilibrium thermodynamic observables that rely on the idea of typicality and derivatives thereof such as the finite-temperature Lanczos method (FTLM). As numerical examples quantum spin systems are studied. Our initial aim was to identify pathological examples or circumstances, such as strong frustration or unusual densities of states, where these methods could fail. It turned out that all investigated systems allow such approximations. Only at temperatures of the order of the lowest energy gap is the convergence somewhat slower in the number of random vectors over which observables are averaged