analysis of H-matrix inversion

Data-sparse matrix formats are an important tool for the numerical treatment of partial differential equations. Especially discretizations of integral equations such as the boundary element method have only become competitive through the incorporation of data-sparse matrix formats, since these permit to reduce the computational complexity from quadratic (in the problem size) to log-linear complexity. The fast multipole method is a prominent example and counted among the top 10 algorithms of the 20th century. Early work in this field has focussed primarily on identifying suitable data formats and designing fast matrix-vector multiplication algorithms for use in iterative solvers. More recently, the quest for algorithms that realize more complex functions such the inversion or factorization has emerged as a highly active research area. Applications of (approximate) inverses include their use as preconditioners, especially for problem classes for which the well-established ones struggle. The theme of the project is the analysis of the matrix inversion (and related problems of factorizations) in the data-sparse matrix format of H-matrices, which are due to Hackbusch. H-matrices are blockwise low rank matrices with the additional feature that the blocks feature a tree structure, which is exploited algorithmically and is a key ingredient for a fast (approximate) arithmetic that includes multiplication and inversion. The block structure is, of course, problem- dependent. A fundamental question is to identify the block structure which can accommodate the inverse at high accuracy. We study stiffness matrices arising from discretizations of partial differential equations. For standard block structures and for scalar second order elliptic equations it is known that the block structure can, in principle, accommodate the inverse of the stiffness matrix. The project extends these results in several directions: First, we aim at corresponding results for systems such as the Lame system, the Stokes system, the Maxwell equations as well as multifield problems. The unifying property of these systems is their ellipticity so that similar block structures can be employed. Wave propgation problems require completely different block structures.We extend the approximability results for the inverse to high wavenumber Helmholtz problems, which arise in accoustics. Finally, the approximability question discussed so far is only the first important step towards an analysis of the H-matrix inversion algorithms. Therefore, we will attempt an actual error analysis of an inversion algorithm in the setting of finite element discretizations of symmetric positive definite systems.

Many problems in mathematical physics or engineering are described by non-local operators. A classical example are many-particle systems with N particles occuring, e.g., in astrophysics where the particles are coupled gravitatively. The computational complexity is quadratic in the number N of particles. The 'fast multipole algorithm' of Rokhlin and Greengard, one of the top-10-algorithms of the 20th century reduces the work to an essentially linear one in N. Mathematically, the algorithm realizes a data-sparse representation of an N by N matrix with O(N) data. In 1999 W. Hackbusch generalized this class of matrices to the class of H-matrices, which are blockwise low-rank matrices but have additionally a hierarchical structure that allowed him to endow the class with a fast arithmetic (additional, multiplication, inversion). The arithmetic is necessarily approxiamte, so that the mathematical question arises, for which matrices the inverse can be approximated well from the class of H-matrix. We showed in the project that (Galerkin) discretizations of elliptic problems lead to matrices whose inverses can be approximated very well by H-matrices. A number of elliptic boundary value problems was studied in more detail in the project. A classical problem class are elliptic partial differential equations (PDEs) posed in full space. Numerically, this problem is treated by discretizing the PDE in a bounded computational domain, realizing the PDE in the remaining unbounded region by an integral equation, and finally coupling the two discretizations. For discretizations of the coupled problem based on quasi-uniform meshes the inverse of the system matrix can be approximated very well by H-matrices. The situation is the same for discretizations of Maxwell's equations, which describe certain electromagnetic phenomen. For scalar elliptic boundary value problems, the project could generalize the existing theory based on quasi-uniform meshes to a large class of highly non-uniform meshes. The ability to deal with such general meshes is essential for the numerical treatment of elliptic problems as it permits to resolve local features of the solution or the prescribed geometry. Another very different class of problems arises from the question to interpolate "smoothly" between given data. An important class of interpolating functions are radial basis functions (RBFs). For the popular class of 'polyharmonic splines', the project showed that the inverse of the interpolation matrix can be represented in the H-matrix format.