PaperlayerView exact revision on arXiv

Source-linked AI summary

A fast direct solver for structured linear systems by recursive skeletonization

Kenneth L. Ho, Leslie Greengard

arXiv:1110.3105v2math.NA

TL;DR

Large dense systems are expensive for classical direct methods, while iterative methods can struggle with ill-conditioning and repeated right-hand sides. The paper addresses this with recursive interpolative-decomposition compression embedded in a sparse system, then factors and applies the inverse in two phases. For non-oscillatory boundary integral equations, the abstract reports O(N) two-dimensional phases and three-dimensional costs of O(N^{3/2}) precomputation and O(N log N) solution, with about 100x solution-phase speedup over fast multipole methods.

  • Problem

    Iterative methods can become costly or unpredictable for ill-conditioned systems and do not effectively reuse work across multiple right-hand sides.

  • Method

    Recursive skeletonization uses interpolative decomposition and auxiliary variables to embed a compressed approximation of a dense structured matrix into a sparse system for factorization and inverse application.

  • Results

    For non-oscillatory boundary integral equations, the abstract reports O(N) complexity for both phases in 2D, while the implementation has O(N^{3/2}) precomputation and O(N log N) solution costs in 3D.

  • Takeaways & Limitations

    The solver is particularly suited to problems requiring many iterations or multiple right-hand sides, where its rapid inverse application can be reused.

  • Takeaways & Limitations

    A principal limitation is growth in factorization cost in 3D or higher, preventing optimal O(N) or nearly optimal O(N log N) complexity.

Abstract

from arXiv · show

We present a fast direct solver for structured linear systems based on multilevel matrix compression. Using the recently developed interpolative decomposition of a low-rank matrix in a recursive manner, we embed an approximation of the original matrix into a larger, but highly structured sparse one that allows fast factorization and application of the inverse. The algorithm extends the Martinsson/Rokhlin method developed for 2D boundary integral equations and proceeds in two phases: a precomputation phase, consisting of matrix compression and factorization, followed by a solution phase to apply the matrix inverse. For boundary integral equations which are not too oscillatory, e.g., based on the Green's functions for the Laplace or low-frequency Helmholtz equations, both phases typically have complexity O(N) in two dimensions, where $N$ is the number of discretization points. In our current implementation, the corresponding costs in three dimensions are $O(N^{3/2})$ and $O(N \log N)$ for precomputation and solution, respectively. Extensive numerical experiments show a speedup of $\sim 100$ for the solution phase over modern fast multipole methods; however, the cost of precomputation remains high. Thus, the solver is particularly suited to problems where large numbers of iterations would be required. Such is the case with ill-conditioned linear systems or when the same system is to be solved with multiple right-hand sides. Our algorithm is implemented in Fortran and freely available.

1. Introduction.

The paper presents a recursive-skeletonization direct solver for structured dense systems, targeting iterative-method weaknesses in ill-conditioned problems and repeated solves. It compresses the matrix, embeds it in a sparse framework, and applies the inverse rapidly after factorization.

  • Motivation: Ill-conditioning can make iterative solution times unpredictable, whereas direct methods remain relatively insensitive to system conditioning.The paper cites near-resonance problems, close-to-touching geometries, and large material contrasts as examples.
  • Motivation: Direct methods efficiently reuse a factorization for multiple right-hand sides, while many iterative methods treat each right-hand side as a new problem.The paper identifies scattering, optimization, and time-dependent processes in fixed geometries as examples.
  • Contribution: The algorithm constructs a compressed, data-sparse representation of the inverse for matrices with block low-rank structure, such as those arising from integral equations.The procedure assumes block low-rank structure similar to that used by fast matrix-vector product methods, without making a priori rank assumptions.
  • Contribution: Compression and auxiliary variables embed an approximation of the dense matrix into a sparse matrix framework compatible with standard sparse matrix technology.This sparse embedding is the central mechanism used to obtain a fast direct solver.
  • Method: Recursive skeletonization applies interpolative decomposition through a multilevel procedure to compute the compressed matrix representation.The compressed representation supports fast matrix-vector multiplication and matrix inversion.
  • Method: The direct solver separates precomputation, which constructs and factors the compressed representation, from solution, which rapidly applies the inverse to a right-hand side.The paper reports that the solution phase can often beat a single fast multipole method call by several orders of magnitude.

2. Preliminaries.

The paper formalizes hierarchical block-separable matrices and explains how interpolative decomposition compresses their low-rank off-diagonal interactions. A tree organizes contiguous index subsets, enabling recursive compression while diagonal blocks may remain full-rank.

  • Block structure: Classical Gaussian elimination requires O(N^3) work for a dense linear system.The system is partitioned into blocks with dimensions summing to N.
  • Block structure: A matrix is block separable when each off-diagonal submatrix factors through three low-rank matrices.The corresponding off-diagonal block rows and columns have low rank, while the diagonal blocks need not be compressed.
  • Hierarchical structure: The hierarchy partitions contiguous index blocks into children, with each child corresponding to a contiguous subset of indices.Nodes at finer levels are concatenated to form their parent index sets.
  • Hierarchical structure: Hierarchical block separability requires low-rank off-diagonal block rows and columns at every level of an index hierarchy.This structure arises, for example, in integral-equation discretizations with non-oscillatory kernels.
  • Interpolative decomposition: The interpolative decomposition approximates A as BP, retaining a subset of A's columns in B and encoding interpolation coefficients in P.Its approximation error is comparable to the next singular value, and P contains an identity submatrix.
  • Interpolative decomposition: Applying the interpolative decomposition to the transpose compresses the row space, producing analogous row and column representations.The decomposition can adaptively determine rank k to achieve a specified relative precision.
  • Interpolative decomposition: Random-sampling compression typically requires O(mn log k + k^2n) operations, replacing the earlier deterministic O(kmn) approach.The paper uses this newer compression technology in its implementation.

3. Algorithm.

The algorithm recursively compresses off-diagonal matrix blocks using interpolative decompositions, regrouping skeleton matrices up an index tree while extracting diagonal blocks. Proxy surfaces localize compression for suitable potential-field kernels, and the resulting representation supports fast matrix-vector multiplication and sparse direct inversion.

  • 3.1. Hierarchical matrix compression.: The one-level scheme compresses each off-diagonal block row and column with interpolative decomposition, then approximates A_ij as L_iS_ijR_j.The skeleton submatrix S_ij is selected by the row and column projections.
  • 3.1. Hierarchical matrix compression.: Recursive skeletonization repeatedly extracts diagonal blocks, recompresses off-diagonal blocks, and regroups them while ascending the index tree.The process terminates when regrouping reaches the root.
  • 3.2. Accelerated compression via proxy surfaces.: Global compression steps dominate cost when no additional matrix knowledge is available, motivating local acceleration through proxy-surface constructions.The proxy-surface approach is applicable to several potential-theory kernels, with extra care required for oscillatory problems.
  • 3.2. Accelerated compression via proxy surfaces.: For potential-field matrices, proxy surfaces replace distant sources with a constant-size representation, reducing compression to neighboring points plus proxy-surface points.This follows because the distant-source field is harmonic inside the proxy surface and can be replicated there; the interaction rank depends on the desired precision.
  • 3.3. Compressed matrix-vector multiplication.: The compressed representation applies matrix-vector products through an upward pass of column projections followed by a downward pass of row projections.These passes hierarchically compress input data to outgoing skeletons and project incoming skeletons back to the output.
  • 3.4. Compressed matrix inversion.: The compressed system is embedded in a larger sparse representation whose recursive factorization enables direct inversion with limited fill-in.The skeleton system is expanded recursively, and the resulting inverse representation preserves the compressed structure with minimal fill-in.

4. Complexity analysis.

The analysis derives costs for recursive compression, matrix-vector products, factorization, and inverse application under block low-rank structure. It also notes that UMFPACK improves practical stability but complicates detailed fill-in and complexity analysis when pivoting is used.

  • Hierarchy: λ ∼ (1/d) log N levels arise because each subdivision reduces the number of blocks by roughly 2^d, reaching O(1) blocks at the coarsest level.The finest level has p1 ∼ N blocks, while pλ = O(1).
  • Rank structure: The interaction rank k_l is determined by recursively exposing well-separated source subblocks, assuming separated interactions have constant rank at fixed precision.This assumption controls the compressed block dimensions across levels.
  • Matrix-vector multiplication: N for d = 1, N log N for d = 2, and N^2(1−1/d) for d > 2 are the total costs for matrix-vector multiplication.The estimate combines levelwise applications of D^(l), L^(l), R^(l), and the coarsest sparse operation.
  • Factorization: The factorization cost has the same complexity as equation (4.1), obtained by summing levelwise construction costs and the final inversion of S.At each level, constructing D^−1 and Λ costs O(p_l n_l^3), while the final solve costs O((p_λ k_λ)^3).
  • Inverse application: N for d = 1, N log N for d = 2, and N^2(1−1/d) for d > 2 also describe the cost of applying the inverse.The inverse has the same dimensions as the corresponding compressed factors and therefore matches matrix-vector multiplication complexity.
  • Practical considerations: Detailed UMFPACK cost analysis is omitted; pivoting can introduce fill-in, making complexity and storage analysis more involved.The sparse approach may provide backward stability when unpivoted Gaussian elimination is unstable.

5. Error analysis.

The error analysis bounds solution error for compressed matrices in terms of compression accuracy and conditioning. It also notes mild observed error propagation across levels and focuses subsequent examples on representative applications.

  • Error propagation: The local ID precision need not equal the final global precision because errors may accumulate across levels, although observed propagation is mild.The paper states this behavior based on its experiments and prior work.
  • Error bound: 2ϵκ(A)/(1−ϵκ(A)) bounds the relative solution error when the compressed system is applied under the stated condition-number requirement.The bound explicitly depends on the condition number κ(A).
  • Well-conditioned case: For a well-conditioned matrix with κ(A) = O(1), such as a second-kind integral-equation discretization, the error estimate simplifies accordingly.The supplied passage identifies second-kind discretizations as an example of this regime.
  • Examples: The numerical study examines Laplace and Helmholtz kernels in two and three dimensions, including timing, molecular electrostatics, and multiple scattering.These experiments evaluate the algorithm's efficiency and flexibility.

6. Numerical examples.

The numerical examples use tree-structured discretizations and proxy compression to evaluate recursive skeletonization for matrix-vector products. They report dimensions, timings, errors, storage, and comparisons with LAPACK/ATLAS and FMM.

  • Experimental setup: Quadtrees in 2D and octrees in 3D uniformly subdivide blocks to size O(1), while empty boxes are truncated adaptively.Proxy compression is used throughout the described setup.
  • Reported metrics: Table 6.1 reports N, row and column skeleton dimensions, compression and matrix-vector timings, relative error, and compressed-matrix storage for the 2D Laplace surface case at ϵ = 10^-9.The table defines Tcm, Tmv, E, and M as measured quantities.
  • Experimental setup: Proxy surfaces use a constant number of discretization points independent of N for Laplace kernels, with Helmholtz counts also depending on frequency.The Helmholtz choice follows Nyquist-Shannon sampling consistency.
  • Baselines: The comparisons use LAPACK/ATLAS and FMM, with FMM timings obtained from the open-source FMMLIB package.The cited FMMLIB implementation excludes plane-wave and diagonal-translation optimizations.
  • Evaluation goals: Recursive skeletonization is first evaluated as a generalized FMM for rapid matrix-vector products.The direct-solver experiments are introduced separately from this matrix-vector application study.

6.1. Generalized fast multipole method.

Recursive skeletonization compresses low-rank kernel matrices through multilevel interpolative decomposition, enabling fast matrix-vector products. For Laplace and low-frequency Helmholtz problems, its observed scaling is favorable, though compression costs can dominate.

  • Recursive skeletonization uses multilevel interpolative decomposition to construct a compressed representation of structured kernel matrices.
  • For 3D Laplace surface geometries, compression and application are expected to cost O(N^3/2) and O(N log N), respectively.
  • The algorithm scales as predicted and can outperform uncompressed matrix-vector products, while remaining competitive with the O(N) FMM over tested 3D sizes.
  • Compression takes one to two orders of magnitude longer than FMM application, but compressed matrix-vector products are consistently small.
  • At N = 131072 in the 2D volume case, Tcm increases from 41 s at ϵ = 10^-3 to 339 s at ϵ = 10^-9, while Tmv increases from 0.09 s to 0.27 s.
  • At low frequency, Helmholtz performance is similar to Laplace performance because both kernels are essentially non-oscillatory.
  • As frequency increases, compression efficiency deteriorates because matrix-block ranks grow; no asymptotic gain remains in the high-frequency regime.

6.2. Recursive skeletonization as a direct solver.

As a direct solver, recursive skeletonization compresses and factors boundary-integral systems, producing rapid inverse applications. It is competitive with FMM/GMRES in tested 2D and medium-scale 3D cases, with especially large solve-time advantages after precomputation.

  • The direct-solver experiments use second-kind boundary-integral discretizations for Laplace and Helmholtz equations in two and three dimensions.
  • In 2D Laplace tests, recursive skeletonization has linear complexity and approaches the performance of the O(N) FMM/GMRES solver.
  • At N = 131072, total recursive-skeletonization time is TRS = 8.5 s versus TFMM = 6.9 s, while inverse application takes Tsv = 0.07 s.
  • At N = 20480 in 3D, recursive skeletonization remains competitive despite O(N^3/2) complexity, with TRS = 850 s and TFMM = 839 s at ϵ = 10^-9.
  • The 3D solve-time advantage remains large because each inverse application costs only O(N log N), yielding TFMM/Tsv ∼2500 at N = 20480 and ϵ = 10^-9.
  • For low-frequency Helmholtz problems, recursive skeletonization can be faster than FMM/GMRES at low to moderate frequencies when memory requirements are not excessive.
  • In 3D Helmholtz tests, the compressed solve remains rapid, with TFMM/Tsv ∼2000 at N = 20480.

6.3. Molecular electrostatics.

The molecular-electrostatics example applies recursive skeletonization to a boundary-integral formulation of polarization response. For a DNA surface, precomputation is costly but subsequent inverse applications are very fast.

  • The molecular model decomposes the electrostatic potential into source and piecewise harmonic polarization components.
  • The polarization response is represented as a single-layer potential, leading to a boundary integral equation for the polarization charge.
  • The DNA surface contains N = 19752 triangles and uses partial charges assigned to n = 486 heavy atoms.
  • At precision ϵ = 10^-3, the recursive-skeletonization computation takes TRS = 592 s, while each inverse application takes Tsv = 0.08 s.
  • The corresponding FMM/GMRES solve takes TFMM = 27 s for this DNA example.

6.4. Multiple scattering.

For multiple scattering, recursive skeletonization supplies block-diagonal inverse operators as preconditioners for FMM/GMRES. This sharply reduces iteration counts and total solution time across configurations with varying scatterer separation.

  • The scattering problem uses a second-kind integral equation for sound-hard obstacles and solves it with FMM/GMRES.
  • The preconditioner is block diagonal, with each isolated-scatterer inverse computed using recursive skeletonization.
  • The test uses two identical smooth scatterers, each ten wavelengths in size, at varied horizontal separation distances.
  • After 0.76 s of precomputation, preconditioning reduces the iteration count from nFMM ∼700 to nRS ∼10 for each configuration.
  • The total solution-time ratio is ∼60 for unpreconditioned versus preconditioned FMM/GMRES solves.
  • The broader conclusion is that recursive skeletonization is competitive when kernels are not too oscillatory or 3D systems are not too large, and is especially useful for ill-conditioning or multiple right-hand sides.

7. Generalizations and conclusions.

The approach remains effective but its factorization cost grows in three or more dimensions, limiting optimal complexity. The algorithm is naturally parallelizable, while reducing higher-dimensional costs remains an active research direction.

  • 7. Generalizations and conclusions.: Hierarchical compression approaches such as HSS matrices, H-matrices, and skeletonization can overcome the higher-dimensional factorization obstacle.
  • 7. Generalizations and conclusions.: The algorithm is naturally parallelizable because many block-sweep computations can be processed independently.This applies clearly to recursive skeletonization using proxy surfaces, with a possible O(log N) performance loss across hierarchy levels.
Loading 1110.3105v2…