Bayesian computations and efficient algorithms for computing functions of large, sparse matrices

Abstract

The need for computing functions of large, sparse matrices arises in Bayesian spatial models where the computations using Gaussian Markov random fields require the evaluation of G -1 and G -1/2 for the precision matrix G and in the geostatistical approach where approximations of R -1 and R 1/2 are needed for the covariance matrix R . In both cases, good approximations to the desired matrix functions are required over a range of probable values of a vector v drawn randomly from a given population, as occurs in simulation techniques for finding posterior distributions such as Markov chain Monte Carlo. Consequently, it is preferable that the complete matrix function approximation be determined rather than for its action on a given v . The aim of this work is to find low degree polynomial approximations p( A ) such that e = ? f( A ) - p( A ) ? 2 is small in some sense on the spectral interval [a,b], where the extreme eigenvalues a and b are calculated using Krylov subspace approximation. Algorithms based on low order near-minimax polynomial approximations are proposed for the required matrix functions for a typical case study in computational Bayesian statistics, where a good balance between accuracy and computationally efficiency is achieved.