Matrix-free Krylov iteration for implicit convolution of numerically low-rank data

Evaluating the response of a linear shift-invariant system is a problem that occurs frequently in a wide variety of science and engineering problems. Calculating the system response via a convolution may be done efficiently with Fourier transforms. When one must compute the response of one system to mm input signals, or the response of mm systems to one signal, it may be the case that one may approximate all system responses without having to compute all mm Fourier transforms. This can lead to substantial computational savings. Rather than process each point individually, one may only process basis vectors that span the output data space. However, to get a low-error approximation, it is necessary that the output vectors have low numerical rank if they were assembled into a matrix. We develop theory that shows how the singular value decay of a matrix ΦAΦA that is a product of a convolution operator ΦΦ and an arbitrary matrix AA depends in a linear fashion on the singular value decays of ΦΦ and AA. We propose gap-rank, a measure of the relative numerical rank of a matrix. We show that convolution cannot destroy the numerical low-rank-ness of ΦAΦA data with only modest assumptions. We then develop a new method that exploits low-rank problems with block Golub–Kahan iteration in a Krylov subspace to approximate the low-rank problem. Our method can exploit parallelism in both the individual convolutions and the linear algebra operations in the block Golub–Kahan algorithm. We present numerical examples from signal and image processing that show the low error and scalability of our method.