Scalable Matrix-Free Adaptive Product-Convolution Approximation for Locally Translation-Invariant Operators

Abstract

We present an adaptive grid matrix-free operator approximation scheme based on a "product-convolution" interpolation of convolution operators. This scheme is appropriate for operators that are locally translation-invariant, even if these operators are high-rank or full-rank. Such operators arise in Schur complement methods for solving partial differential equations (PDEs), as Hessians in PDE-constrained optimization and inverse problems, as integral operators, as covariance operators, and as Dirichlet-to-Neumann maps. Constructing the approximation requires computing the impulse responses of the operator to point sources centered on nodes in an adaptively refined grid of sample points. A randomized a-posteriori error estimator drives the adaptivity. Once constructed, the approximation can be efficiently applied to vectors using the fast Fourier transform. The approximation can be efficiently converted to hierarchical matrix (H-matrix) format, then inverted or factorized using scalable H-matrix arithmetic. The quality of the approximation degrades gracefully as fewer sample points are used, allowing cheap lower quality approximations to be used as preconditioners. This yields an automated method to construct preconditioners for locally translationinvariant Schur complements. We directly address issues related to boundaries and prove that our scheme eliminates boundary artifacts. We test the scheme on a spatially varying blurring kernel, on the non-local component of an interface Schur complement for the Poisson operator, and on the data misfit Hessian for an advection dominated advection-diffusion inverse problem. Numerical results show that the scheme outperforms existing methods.

hIPPYlib-MUQ: A Bayesian Inference Software Framework for Integration of Data with Complex Predictive Models under Uncertainty

“…For example, the results of section 5.1.5 for the Poisson coefficient inverse problem indicate that 𝑂 (10 6 ) PDE solves may still be required even with the most efficient MCMC methods. In such cases, hIPPYlib-MUQ can be used as a prototyping environment to study new methods that further exploit problem structure, for example through the use of various reduced models (e.g., [20]) or via advanced Hessian approximations that go beyond low rank [2,3].…”

Fast wavelet decomposition of linear operators through product-convolution expansions

“…In the context of PDEs, product-convolution expansions are encountered in Schur complement methods for solving PDEs, Dirichlet-to-Neumann maps and in PDE-constrained optimization as Hessians. In [1], the functions (u k , v k ) m k=1 are chosen as above, but instead of being defined on a Euclidean grid, locations (y k ) m k=1 are adaptively sampled to best approximate the operator with a fixed number m.…”

“…Product-convolution expansions have appeared at least 3 decades ago and have been given different names in various fields, see e.g. [5,25,15,17,10,1]. A remarkable aspect of these expansions is their essential role in the interpolation of linear operators from scattered impulse responses.…”

Hierarchical Matrix Approximations of Hessians Arising in Inverse Problems Governed by PDEs

“…Also exploiting the pseudodifferential structure of seismic inversion Hessians is the matrix probing method of [31], which approximates the Hessian (and its inverse) with basis matrices stemming from the Hessian's symbol and find their coefficients by probing the Hessian in random directions. Recently, methods that exploit the local translation invariance of Hessians have been introduced [74,3]. The adaptive product-convolution approximation in particular is demonstrated to be robust to the Peclet number for advection-dominated transport and the frequency for an auxiliary operator that arises in connection with KKT preconditioning [3] for a wave inverse problem [2].…”

“…Recently, methods that exploit the local translation invariance of Hessians have been introduced [74,3]. The adaptive product-convolution approximation in particular is demonstrated to be robust to the Peclet number for advection-dominated transport and the frequency for an auxiliary operator that arises in connection with KKT preconditioning [3] for a wave inverse problem [2]. Here we focus on comparisons with low rank approximation, and defer comparison to these other methods in appropriate contexts for future work.…”