A Fast Randomized Algorithm for Computing a Hierarchically Semiseparable Representation of a Matrix
SIAM Journal on Matrix Analysis and Applications
| Hi-index | 0.00 |
In 1985, Vladimir Rokhlin introduced the fast multipole method (FMM) for the solution of the Laplace equation. The method has since been named, by Computing in Science & Engineering magazine, one of the 20th century's ten most influential algorithms. In 1990, Rokhlin introduced a version of FMM for the Helmholtz equation that has greatly influenced algorithm development in the computational electromagnetics community. Unfortunately, scattering FMM is numerically unstable. Its instability can be traced to the asymptotic behavior of the multipole basis functions. After demonstrating the instability, I introduce measures that eliminate it without sacrificing the method's efficiency. Relative accuracies approaching 10−16 are possible, even for scattering obstacles with diameters smaller than 10−200 wavelengths. When applied to the solution of PDEs, multipole methods are embedded into iterative solvers such as GMRES. In many circumstances, the discretization produces a poorly conditioned algebraic system, and iterative solvers are slow to converge. I introduce a fast direct solver which avoids that difficulty. The multipole structure is embedded into a large sparse system, to which a standard sparse solver is applied. I demonstrate scattering problems for which the new solver is clearly superior to either dense Gaussian elimination or iterative FMM. To apply FMM, the PDE is first reformulated as a weakly singular integral equation. Multipole methods constrain the integral discretization in a way that makes a high-order numerical solution difficult. The existing high-order discretizations are numerically unstable. I present a stable discretization of arbitrarily high order that satisfies the multipole constraints. I routinely solve scattering problems using rules with order 32, and I have constructed stable discretizations with orders as high as 288. To illustrate these techniques, I exhibit numerical solutions in two space dimensions. All computations have been carried out in MATLAB . Scattering obstacles with diameters greater than 1000 wavelengths can be comfortably treated with a personal computer.