  We present an efficient matrix-free point spread function (PSF) method for
approximating operators that have locally supported non-negative integral
kernels. The method computes impulse responses at scattered points, and
interpolates these impulse responses to approximate integral kernel entries.
Impulse responses are computed by applying the operator to Dirac comb batches
of point sources, which are chosen via an ellipsoid packing procedure.
Evaluation of kernel entries allows us to construct a hierarchical matrix
approximation of the operator, which is used for further matrix computations.
We illustrate the end-to-end method on a blur problem, then use the method to
build preconditioners for the Hessian in two inverse problems governed by
partial differential equations (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. But Hessian impulse
responses typically become more local as the numerical rank increases, which
benefits the PSF method. Numerical results reveal that the PSF preconditioner
clusters the spectrum of the preconditioned Hessian near one, yielding roughly
5x-10x reductions in the required number of PDE solves, as compared to
regularization preconditioning and no preconditioning. We also present a
numerical study for the influence of various parameters (that control the shape
of the impulse responses) on the effectiveness of the advection-diffusion
Hessian approximation. The results show that the PSF-based preconditioners are
able to form good approximations of high-rank Hessians using a small number of
operator applications.
