Source-linked AI summary
A Fast Block Low-Rank Dense Solver with Applications to Finite-Element Matrices
Amirhossein Aminfar, Sivaram Ambikasaran, Eric Darve
TL;DR
Large finite-element systems can be difficult to solve, and dense frontal matrices are a costly stage of multifrontal elimination. The paper combines an HODLR solver with graph-guided BDLR approximations and iterative preconditioning; its benchmarks report faster near-machine-accuracy solves than conventional LU in almost all cases.
Problem
Method
The paper approximates frontal matrices with HODLR structure, uses graph-based boundary-distance row and column selection for BDLR, and applies a low-accuracy HODLR solver as a preconditioner.
Source: §2.2 · p. 6, §5.3 · pp. 12–13Results
BDLR and ACA preconditioners reached near-machine accuracy faster than conventional LU in almost all cases; BDLR achieved higher speedup on the largest structured and unstructured cases.
Source: §7.3 · pp. 17–18Takeaways & Limitations
The HODLR block-low-rank representation significantly reduces multifrontal solver memory consumption, a stated bottleneck on existing hardware.
Source: §8 · p. 20In 3D, off-diagonal rank grows slowly like n^1/2 rather than remaining constant, although the authors report that HODLR remains competitive.
Source: §8 · p. 20
Abstract
from arXiv · showhide
This article presents a fast solver for the dense "frontal" matrices that arise from the multifrontal sparse elimination process of 3D elliptic PDEs. The solver relies on the fact that these matrices can be efficiently represented as a hierarchically off-diagonal low-rank (HODLR) matrix. To construct the low-rank approximation of the off-diagonal blocks, we propose a new pseudo-skeleton scheme, the boundary distance low-rank approximation, that picks rows and columns based on the location of their corresponding vertices in the sparse matrix graph. We compare this new low-rank approximation method to the adaptive cross approximation (ACA) algorithm and show that it achieves betters speedup specially for unstructured meshes. Using the HODLR direct solver as a preconditioner (with a low tolerance) to the GMRES iterative scheme, we can reach machine accuracy much faster than a conventional LU solver. Numerical benchmarks are provided for frontal matrices arising from 3D finite element problems corresponding to a wide range of applications.
1. Introduction
Large finite-element systems can be difficult for iterative methods when ill-conditioned, while conventional direct solves are costly, especially in 3D. The paper targets dense frontal matrices with a fast direct solver used as a preconditioner for an iterative solve.
1. Introduction
- Dense frontal matrices are a computational bottleneck in multifrontal solves, while conventional direct solvers scale as O(N^2) for 3D meshes.The motivating systems come from finite-element discretizations of elliptic PDEs and may be ill-conditioned, challenging iterative methods.Source: §1 · p. 1
- The proposed strategy uses a fast direct solver as a preconditioner to combine iterative and direct solvers’ speed, accuracy, and robustness for ill-conditioned matrices.The direct solver is intended to act as a highly accurate preconditioner.Source: §1 · pp. 1–2
- The paper introduces an HODLR direct solver, develops BDLR for approximating frontal-matrix off-diagonal blocks, and benchmarks iterative preconditioning and FETI-DP applications.The benchmarks cover frontal matrices from varied 3D finite-element problems.Source: §1 · p. 2
2. Previous Work
The review contrasts hierarchical dense-matrix structures and their solvers, then surveys approaches to reducing multifrontal costs. It positions the paper’s algebraic BDLR approximation and HODLR frontal solver as applicable to varied mesh structures.
2.1. Fast Direct Solvers for Dense Hierarchical Matrices
- Hierarchical matrices represent certain dense matrices through recursively partitioned blocks, with admissible subblocks stored at low rank; H-matrices also subdivide off-diagonal blocks into low- and full-rank blocks.The H-matrix factorization literature reports O(n log^2(n)) complexity; other reviewed HSS methods use factorization, Woodbury-based solves, or recursive constructions.Source: §2.1 · pp. 2–4
- HODLR uses low-rank off-diagonal blocks without nested bases and has O(n log^2 n) solver cost; HSS adds nesting and can achieve O(n) complexity.The simpler HODLR structure can have lower average off-diagonal rank than HSS, while HSS’s nested structure is more complicated to work with.Source: §2.1 · pp. 3–4
- The paper proposes BDLR, an algebraic scheme that selects rows and columns from the sparse matrix graph, and compares it with ACA for off-diagonal low-rank approximation.Unlike analytic methods such as Chebyshev approximation or multipole expansions, the algebraic approach does not require a matrix defined by an analytical kernel.Source: §2.1 · p. 4
2.2. Fast Direct Solvers for Sparse Matrices
- Sparse Cholesky and multifrontal direct methods have high costs for large finite-element systems, with nested dissection scaling as O(N^1.5) in 2D and O(N^2) in 3D.The efficiency of sparse Cholesky depends strongly on mesh-node ordering.Source: §2.2 · p. 4
- Earlier multifrontal compression work includes HSS-based methods and BLR frontal formats, with reported scope or implementation gaps for some approaches.Some reviewed HSS work considered regular, well-shaped 2D meshes; the cited BLR discussion does not cover extend-add operations or a full multifrontal solver.Source: §2.2 · pp. 5–6
- The paper observes HODLR structure in multifrontal frontal matrices and presents a dense HODLR solver that can be applied without constructing HSS trees and across mesh structures.The benchmarks focus on challenging frontal matrices from large mesh problems; unstructured meshes are harder to precondition than structured meshes, and 3D problems have higher off-diagonal ranks than 2D problems.Source: §2.2 · p. 6
3. An Iterative Solver with Direct Solver Preconditioning
The paper investigates low-accuracy HODLR direct solves as preconditioners for GMRES. It reports faster convergence to high accuracy than conventional LU and high-accuracy HODLR solves.
3. An Iterative Solver with Direct Solver Preconditioning
- The low-accuracy HODLR direct solver is used as a preconditioner for GMRES.The authors state that the preconditioning approach can also be applied to other iterative solvers, including conjugate gradient.Source: §3 · pp. 6–7
- The approach is framed as a way to use a relatively low-accuracy direct solve within an iterative scheme.The supplied table passage identifies a summary of fast sparse direct solvers but provides no table contents.Source: §3 · p. 6, Table 2 · p. 7
- The preconditioned approach is reported to be much faster than both conventional LU and a high-accuracy direct HODLR solver.The passage presents this as the paper’s result for the proposed approach.Source: §3 · p. 7
4. A Fast Direct Solver for HODLR Matrices
The section presents a recursive direct solver that exploits HODLR structure in dense frontal matrices. Its cost depends on off-diagonal ranks: decaying ranks can make the solver efficient enough to compete with HSS solvers.
4. A Fast Direct Solver for HODLR Matrices
- Approximating dense frontal matrices as HODLR matrices reduces factorization cost from conventional dense O(n^3) to O(n log^2(n)).Source: §4 · p. 7
4.1. HODLR Matrices
- HODLR matrices have low-rank off-diagonal blocks across multiple levels, with separate low-rank bases that are not nested across levels.Source: §4.1 · pp. 7–8
4.2. Solver Derivation and Algorithm
- The solver recursively factors diagonal blocks independently, using Woodbury-equivalent operations; this requires both diagonal blocks to be nonsingular and suits parallel implementation.Source: §4.2 · pp. 8–9
4.3. Algorithm Summary
- The algorithm summary applies the same recursive procedure at each HODLR tree level and block.Source: §4.3 · p. 9
4.3.1. Factorize
- Factorization approximates off-diagonal blocks and recursively forms Schur-complement terms and the corresponding right-hand sides across levels.Source: §4.3.1 · p. 9
4.3.2. Solve
- The solve phase uses the previously computed factorization, so its work is mainly matrix-matrix products and is small relative to factorization.Source: §4.3.2 · p. 10
4.4. Solver Computational Cost
- At each level, solving costs include two child solves with s+r right-hand sides, Schur-complement construction, and right-hand-side construction.Source: §4.4 · p. 10
- With decreasing off-diagonal rank, HODLR solvers can become very efficient and compete with HSS solvers; constant rank gives O(r^2n log^2(n)) cost.Source: §4.4 · pp. 10–11
5. Low-Rank Approximation Schemes
The section surveys low-rank approximation methods for HODLR off-diagonal blocks and introduces BDLR, which selects rows and columns using graph-boundary distance. This selection addresses limitations of existing methods and is motivated by the observed relationship between graph proximity and large matrix entries.
5. Low-Rank Approximation Schemes
- Chebyshev low-rank approximation is suited only to interactions mediated by smooth kernels.The section identifies this as one drawback among the approximation schemes considered.Source: §5 · p. 11
5.1. ACA Low-Rank Approximation
- Partial-pivoting ACA costs O(r(m + n)) for an m × n matrix but may miss critical rows or columns when leverage scores are nonuniform.It selects pivots by searching for a large entry in one row or column at each step.Source: §5 · p. 11, §5.1 · p. 12
5.3. Pseudo-Skeleton and Boundary Distance Low-Rank Approximation
- BDLR adds a row-and-column selection method to pseudo-skeleton approximation, whose accuracy depends strongly on which rows and columns are chosen.Pseudo-skeleton forms an approximation from selected rows and columns; the paper increases their number until a relative Frobenius norm error meets a tolerance.Source: §5 · p. 11, §5.3 · pp. 12–13
5.2. Randomized Algorithms
- Randomized low-rank methods cost O(n^2) for an n × n matrix without a fast matrix-vector product, which is unavailable in the benchmark settings.The paper therefore limits its practical choices to methods such as Chebyshev, partial-pivoting ACA, and pseudo-skeleton approximation.Source: §5 · p. 11, §5.2 · p. 12
5.3. Pseudo-Skeleton and Boundary Distance Low-Rank Approximation
- For FEM or BEM blocks related to a Green’s function, nearby points correspond to large entries, motivating graph distance as a row-and-column selection heuristic.The stated rationale assumes the Green’s function is smooth away from its singularity at the origin.Source: §5.3 · pp. 12–13
- BDLR selects boundary rows and columns first, then adds rows and columns by increasing graph distance from the boundary until the desired accuracy is reached.The paper associates boundary distance with interactions across row and column sets in the sparse matrix graph.Source: §5.3 · p. 13
6. Application for Multifrontal Solve Process
The section applies the fast dense solver to frontal matrices from sparse multifrontal finite-element problems, including elasticity and FETI-DP settings. It presents the benchmarks as a proof of concept, with structured and unstructured mesh cases.
6. Application for Multifrontal Solve Process
- The solver is applied to frontal matrices from sparse multifrontal finite-element problems, with the goal of speeding the multifrontal process to approximately O(N^4/3).The authors characterize the results as a proof of concept and state that SCOTCH supplies separators handled by the algorithm without changes.Source: §6 · pp. 14–15
6.1. Elasticity Equation for a 3D Beam and a Cylinder Head Geometry
- Both beam and cylinder-head elasticity problems use unstructured meshes, with 10-node tetrahedra for the beam and mixed solid and shell elements for the cylinder head.The cylinder-head solid mesh combines 8-node hexahedral, 6-node pentahedral, and 4-node tetrahedral elements, alongside 3-node shell elements.Source: §6.1 · p. 15
6.2. FETI-DP Solver for a 3D Elasticity Problem
- The FETI-DP benchmarks use two local 3D elasticity matrices: one from a structured cube mesh and one from an unstructured engine mesh.The cube uses 8-node hexahedral elements, while the engine uses 10-node tetrahedral elements.Source: §6.2 · pp. 15–16
7. Numerical Benchmarks
The benchmarks evaluate BDLR and ACA preconditioners on frontal matrices from 3D elasticity problems, including structured and unstructured meshes. Both can reach near-machine accuracy faster than conventional LU in almost all cases, while BDLR has higher speedup than ACA on the largest cases.
7.1. Elasticity Equation for a 3D Beam and a Cylinder Head Geometry
- Beam and cylinder-head fronts show low-rank off-diagonal blocks, with larger pivots generally associated with rows and columns nearer the interaction-graph boundary.For the cylinder head, good approximation at low ϵ requires more rows and columns, corresponding to greater BDLR depth.Source: §7.1 · pp. 16–17
- BDLR off-diagonal rank decays from root to leaf, and the benchmarks include level-by-level factorization, solve, and approximation timings.The study also compares BDLR and ACA convergence as GMRES preconditioners for beam and cylinder-head fronts.Source: §7.1 · p. 17
7.2. FETI-DP Solver for a 3D Elasticity Problem
- The FETI-DP benchmarks compare structured cube and unstructured engine meshes using BDLR and ACA preconditioners on frontal matrices.The cube and engine convergence comparisons are shown for sample fronts.Source: §7.2 · p. 17
- For the structured cube, BDLR error is comparable to optimal SVD error; the unstructured engine mesh has important pivots beyond the points nearest the boundary.For ϵ = 10^-10, the engine approximation needs more points; insufficient sampling can yield small LU pivots and large relative error.Source: §7.2 · p. 17
7.3. Summary
- BDLR accuracy and convergence depend on choosing depth alongside tolerance; lowering ϵ without sufficient depth can degrade the preconditioner.For ACA, the authors note that decreasing tolerance can improve results.Source: §7.3 · p. 18
- BDLR rows and columns are known a priori from the separator-graph structure, a property the authors say could accelerate extend-add in future work.They propose tracking only the required rows and columns to avoid constructing large dense frontal and update matrices.Source: §7.3 · p. 18
- In almost all cases, BDLR- and ACA-preconditioned GMRES reached near-machine accuracy faster than conventional LU; BDLR had higher speedup on the largest cases.The largest cases were the 1.5M structured cube and 2.3M unstructured cylinder head; convergence depends on method parameters.Source: §7.3 · pp. 17–18, Table 3 · p. 18
8. Conclusion and Future Work
The conclusion presents BDLR as a robust low-rank scheme for algebraic Schur complements, tracking selected rows and columns during multifrontal extend-add. It also discusses rank growth, computational cost, and parallel scalability of the HODLR solver.
8. Conclusion and Future Work
- Using block low-rank structures like HODLR significantly reduces memory consumption in multifrontal sparse solvers, where memory is a major hardware bottleneck.Source: §8 · p. 20
- BDLR is presented as a robust low-rank approximation scheme for algebraic Schur complements that tracks selected rows and columns during multifrontal extend-add.The authors say a follow-up paper will show runtime and memory improvements from this tracking.Source: §8 · p. 20
- In 3D, off-diagonal rank grows slowly like n^1/2, but HODLR remains competitive because its cost is O(r^2 n log^2 n) with a relatively small constant.Source: §8 · p. 20
- Recursive solves on independent diagonal blocks give the direct algorithm parallel scalability, potentially favoring it on large parallel computers despite higher flop counts than an optimal O(n) method.Source: §8 · p. 20
- The authors estimate the direct solver starts being faster at r ∼0.4n than an (2/3)n^3 LU factorization algorithm.Source: §8 · p. 20