Point Spread Function Approximation of High-Rank Hessians with Locally Supported Nonnegative Integral Kernels

SIAM Journal on Scientific Computing, Volume 46, Issue 3, Page A1658-A1689, June 2024.
Abstract. We present an efficient matrix-free point spread function (PSF) method for approximating operators that have locally supported nonnegative integral kernels. The PSF-based method computes impulse responses of the operator at scattered points and interpolates these impulse responses to approximate entries of the integral kernel. To compute impulse responses efficiently, we apply the operator to Dirac combs associated with batches of point sources, which are chosen by solving an ellipsoid packing problem. The ability to rapidly evaluate kernel entries allows us to construct a hierarchical matrix (H-matrix) approximation of the operator. Further matrix computations are then performed with fast H-matrix methods. This end-to-end procedure is illustrated on a blur problem. We demonstrate the PSF-based method’s effectiveness by using it to build preconditioners for the Hessian operator arising in two inverse problems governed by PDEs: inversion for the basal friction coefficient in an ice sheet flow problem and for the initial condition in an advective-diffusive transport problem. While for many ill-posed inverse problems the Hessian of the data misfit term exhibits a low-rank structure, and hence a low-rank approximation is suitable, for many problems of practical interest, the numerical rank of the Hessian is still large. The Hessian impulse responses, on the other hand, typically become more local as the numerical rank increases, which benefits the PSF-based method. Numerical results reveal that the preconditioner clusters the spectrum of the preconditioned Hessian near one, yielding roughly [math]–[math] reductions in the required number of PDE solves, as compared to classical regularization-based preconditioning and no preconditioning. We also present a comprehensive numerical study for the influence of various parameters (that control the shape of the impulse responses and the rank of the Hessian) on the effectiveness of the advection-diffusion Hessian approximation. The results show that the PSF-based method is able to form good approximations of high-rank Hessians using only a small number of operator applications.