Sparse Hierarchical Preconditioners Using Piecewise Smooth Approximations of Eigenvectors

TL;DR

This paper introduces a hierarchical preconditioning method that improves the efficiency of iterative solvers for PDE discretizations by focusing on accurately approximating smooth eigenvectors, leading to near-optimal solution times.

Contribution

It presents a novel hierarchical approximate factorization approach that constructs sparse preconditioners with low complexity, tailored to smooth eigenvectors in elliptic PDE problems.

Findings

Achieves $O(n)$ or $O(n \, \log n)$ construction complexity.

Demonstrates optimal $O(n)$ solution times on large elliptic problems.

Preconditioners are exact on near-kernel modes for elasticity equations.

Abstract

When solving linear systems arising from PDE discretizations, iterative methods (such as Conjugate Gradient, GMRES, or MINRES) are often the only practical choice. To converge in a small number of iterations, however, they have to be coupled with an efficient preconditioner. The efficiency of the preconditioner depends largely on its accuracy on the eigenvectors corresponding to small eigenvalues, and unfortunately, black-box methods typically cannot guarantee sufficient accuracy on these eigenvectors. Thus, constructing the preconditioner becomes a problem-dependent task. However, for a large class of problems, including many elliptic equations, the eigenvectors corresponding to small eigenvalues are smooth functions of the PDE grid. In this paper, we describe a hierarchical approximate factorization approach which focuses on improving accuracy on the smooth eigenvectors. The improved accuracy is achieved by preserving the action of the factorized matrix on piecewise polynomial functions of the grid. Based on the factorization, we propose a family of sparse preconditioners with $O(n)$ or $O(n \log{n})$ construction complexities. Our methods exhibit the optimal $O(n)$ solution times in benchmarks run on large elliptic problems of different types, arising for example in flow or mechanical simulations. In the case of the linear elasticity equation the preconditioners are exact on the near-kernel rigid body modes.

Tables

Table 1. Table 7.1: Accuracy on the unit-length eigenvectors corresponding to the largest eigenvalue ( E 1 subscript 𝐸 1 E_{1} ), and the smallest eigenvalue ( E n subscript 𝐸 𝑛 E_{n} ), for Nest-All-All with polynomial compression, and its standard low-rank equivalent, applied to Poisson equation Eq. (25) .

	$E_1=\frac{1}{{\lambda }_1}{\Vert \left(𝐀-{𝐀}_{\ell }\right){𝐯}_1\Vert }_2$		$E_n=\frac{1}{{\lambda }_n}{\Vert \left(𝐀-{𝐀}_{\ell }\right){𝐯}_n\Vert }_2$
$n$	Pcw. quad.	Low-rank	Pcw. quad.	Low-rank
4.0M	$4.5\cdot {10}^{-2}$	$4.4\cdot {10}^{-2}$	$1.2\cdot {10}^0$	$4.3\cdot {10}^2$
6.0M	$4.5\cdot {10}^{-2}$	$4.4\cdot {10}^{-2}$	$1.3\cdot {10}^0$	$5.6\cdot {10}^2$
8.0M	$4.5\cdot {10}^{-2}$	$4.4\cdot {10}^{-2}$	$1.2\cdot {10}^0$	$6.8\cdot {10}^2$
12.0M	$4.5\cdot {10}^{-2}$	$4.4\cdot {10}^{-2}$	$1.1\cdot {10}^0$	$8.9\cdot {10}^2$
16.0M	$4.7\cdot {10}^{-2}$	$4.4\cdot {10}^{-2}$	$1.4\cdot {10}^0$	$1.1\cdot {10}^3$

Equations

$$\mathbf{Ax}=\mathbf{b}$$

$$\Vert {\mathbf{A}}_{\ell }-\mathbf{A}{\Vert }_2\le \epsilon \Vert \mathbf{A}{\Vert }_2$$

$$\Vert \mathbf{I}-{\mathbf{A}}_{\ell }^{-1/2}\mathbf{A}{\mathbf{A}}_{\ell }^{-1/2}{\Vert }_2\le \epsilon$$

$$\mathbf{A}\approx {\mathbf{A}}_{\ell }={\mathbf{G}}_1{\mathbf{B}}_1{\mathbf{G}}_2{\mathbf{B}}_2\cdots {\mathbf{G}}_{\ell -1}{\mathbf{B}}_{\ell -1}{\mathbf{G}}_{\ell }{\mathbf{G}}_{\ell }^T{\mathbf{B}}_{\ell -1}^T{\mathbf{G}}_{\ell -1}^T\cdots {\mathbf{B}}_2^T{\mathbf{G}}_2^T{\mathbf{B}}_1^T{\mathbf{G}}_1^T$$

$$\begin{array}{cc}{\mathbf{\Pi }}^1:=\left(\begin{array}{cccc}1& x_1& y_1& z_1\\ 1& x_2& y_2& z_2\\ ⋮& ⋮& ⋮& ⋮\\ 1& x_n& y_n& z_n\end{array}\right),& {\mathbf{\Pi }}^2:=\left(\begin{array}{cc}\begin{array}{c}{\mathbf{\Pi }}^1\end{array}& \begin{array}{cccccc}{x}_1^2& y_1^2& z_1^2& x_1{y}_1& y_1{z}_1& z_1{x}_1\\ x_2^2& y_1^2& z_1^2& x_2{y}_2& y_2{z}_2& z_2{x}_2\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ x_n^2& y_n^2& z_n^2& x_n{y}_n& y_n{z}_n& z_n{x}_n\end{array}\end{array}\right)\end{array}$$

$$\mathbf{M}=\left(\begin{array}{ccc}{\mathbf{M}}_{\mathcal{BB}}& {\mathbf{M}}_{\mathcal{B}\mathcal{N}}& \\ {\mathbf{M}}_{\mathcal{N}\mathcal{B}}& {\mathbf{M}}_{\mathcal{NN}}& {\mathbf{M}}_{\mathcal{NW}}\\ & {\mathbf{M}}_{\mathcal{WN}}& {\mathbf{M}}_{\mathcal{WW}}\end{array}\right)$$

$$\mathbf{M}=\left(\begin{array}{cc}\mathbf{L}& \\ & \mathbf{I}\end{array}\right)\left(\begin{array}{ccc}\mathbf{I}& {\hat{\mathbf{M}}}_{\mathcal{B}{\mathcal{N}}_1}& {\hat{\mathbf{M}}}_{\mathcal{B}{\mathcal{N}}_2}\\ {\hat{\mathbf{M}}}_{{\mathcal{N}}_1\mathcal{B}}& {\mathbf{M}}_{{\mathcal{N}}_1{\mathcal{N}}_1}& {\mathbf{M}}_{{\mathcal{N}}_1{\mathcal{N}}_2}\\ {\hat{\mathbf{M}}}_{{\mathcal{N}}_2\mathcal{B}}& {\mathbf{M}}_{{\mathcal{N}}_2{\mathcal{N}}_1}& {\mathbf{M}}_{{\mathcal{N}}_2{\mathcal{N}}_2}\end{array}\right)\left(\begin{array}{cc}{\mathbf{L}}^T& \\ & \mathbf{I}\end{array}\right)$$

$$\mathbf{\Phi }=\left(\begin{array}{ccc}{\mathbf{\Phi }}_B& & \\ & {\mathbf{\Phi }}_{N_1}& \\ & & {\mathbf{\Phi }}_{N_2}\end{array}\right)$$

$$\mathbf{N}=\left(\begin{array}{ccc}{\mathbf{L}}^T{\mathbf{\Phi }}_B& {\hat{\mathbf{M}}}_{\mathcal{B}{\mathcal{N}}_1}{\mathbf{\Phi }}_{N_1}& {\hat{\mathbf{M}}}_{\mathcal{B}{\mathcal{N}}_2}{\mathbf{\Phi }}_{N_2}\end{array}\right)$$

$$\mathbf{NP}=\mathbf{QR}=\left(\begin{array}{cc}{\mathbf{Q}}_1& {\mathbf{Q}}_2\end{array}\right)\left(\begin{array}{c}{\mathbf{R}}_1\\ \end{array}\right)$$

$$\mathbf{B}=\left(\begin{array}{cc}\mathbf{LQ}& \\ & \mathbf{I}\end{array}\right)$$

$$\mathbf{M\Phi }$$

$$=$$

$$\mathbf{M}\approx \mathbf{B}\left(\begin{array}{cccc}\mathbf{I}& & {\mathbf{Q}}_1^T{\hat{\mathbf{M}}}_{\mathcal{B}{N}_1}& {\mathbf{Q}}_1^T{\hat{\mathbf{M}}}_{\mathcal{B}{\mathcal{N}}_2}\\ & \mathbf{I}& & \\ {\hat{\mathbf{M}}}_{{\mathcal{N}}_1\mathcal{B}}{\mathbf{Q}}_1& & {\mathbf{M}}_{N_1{\mathcal{N}}_1}& {\mathbf{M}}_{{\mathcal{N}}_1{\mathcal{N}}_2}\\ {\hat{\mathbf{M}}}_{{\mathcal{N}}_2\mathcal{B}}{\mathbf{Q}}_1& & {\mathbf{M}}_{{\mathcal{N}}_2{\mathcal{N}}_1}& {\mathbf{M}}_{{\mathcal{N}}_2{\mathcal{N}}_2}\end{array}\right){\mathbf{B}}^T:=\mathbf{B}{\mathbf{M}}_{(+)}{\mathbf{B}}^T$$

$${\mathbf{\Phi }}_B\leftarrow {\mathbf{Q}}_1^T{\mathbf{L}}^T{\mathbf{\Phi }}_B$$

$$\tilde{\mathbf{M}}:=\mathbf{B}\left(\begin{array}{ccccccc}\mathbf{I}& & {\mathbf{Q}}_1^T{\hat{\mathbf{M}}}_{\mathcal{B}{\mathcal{N}}_1}& {\mathbf{Q}}_1^T{\hat{\mathbf{M}}}_{\mathcal{B}{\mathcal{N}}_2}& \cdots & {\mathbf{Q}}_1^T{\hat{\mathbf{M}}}_{\mathcal{B}{\mathcal{N}}_g}& \\ & \mathbf{I}& & & & & \\ {\hat{\mathbf{M}}}_{{\mathcal{N}}_1\mathcal{B}}{\mathbf{Q}}_1& & {\mathbf{M}}_{{\mathcal{N}}_1{\mathcal{N}}_1}& {\mathbf{M}}_{{\mathcal{N}}_1{\mathcal{N}}_2}& \cdots & {\mathbf{M}}_{{\mathcal{N}}_1{\mathcal{N}}_g}& {\mathbf{M}}_{{\mathcal{N}}_1\mathcal{W}}\\ {\hat{\mathbf{M}}}_{{\mathcal{N}}_2\mathcal{B}}{\mathbf{Q}}_1& & {\mathbf{M}}_{{\mathcal{N}}_2{\mathcal{N}}_1}& {\mathbf{M}}_{{\mathcal{N}}_2{\mathcal{N}}_2}& \cdots & {\mathbf{M}}_{{\mathcal{N}}_2{\mathcal{N}}_g}& {\mathbf{M}}_{{\mathcal{N}}_2\mathcal{W}}\\ ⋮& & ⋮& ⋮& \ddots & ⋮& ⋮\\ {\hat{\mathbf{M}}}_{{\mathcal{N}}_g\mathcal{B}}{\mathbf{Q}}_1& & {\mathbf{M}}_{{\mathcal{N}}_g{\mathcal{N}}_1}& {\mathbf{M}}_{{\mathcal{N}}_g{\mathcal{N}}_2}& \cdots & {\mathbf{M}}_{{\mathcal{N}}_g{\mathcal{N}}_g}& {\mathbf{M}}_{{\mathcal{N}}_g\mathcal{W}}\\ & & {\mathbf{M}}_{\mathcal{W}{\mathcal{N}}_1}& {\mathbf{M}}_{\mathcal{W}{\mathcal{N}}_2}& \cdots & {\mathbf{M}}_{\mathcal{W}{\mathcal{N}}_g}& {\mathbf{M}}_{\mathcal{WW}}\end{array}\right){\mathbf{B}}^T$$

$$\mathbf{M}{\mathbf{\Pi }}_B=\tilde{\mathbf{M}}{\mathbf{\Pi }}_B$$

$$({\mathbf{\Pi }}_B{)}_{ij}=\{\begin{array}{ll}{\mathbf{\Pi }}_{ij}& \text{if }i\in B\\ 0& \text{otherwise. }\end{array}$$

$${\mathbf{A}}_k\mathbf{Y}=\mathbf{AY}$$

$${\mathbf{A}}_{(i-1)}=\left(\prod _{j=1}^{|I_i|}{\mathbf{G}}_{(j)}\right){\mathbf{A}}_{(i-\frac{1}{2})}\left(\prod _{j=|I_i|}^{j=1}{\mathbf{G}}_{(j)}^T\right)={\mathbf{G}}_i{\mathbf{A}}_{(i-\frac{1}{2})}{\mathbf{G}}_i^T$$

$${\mathbf{A}}_{(i-\frac{1}{2})}\approx \left(\prod _{j=1}^{|S_i|}{\mathbf{B}}_{(i)}\right){\mathbf{A}}_{(i)}\left(\prod _{j=|S_i|}^1{\mathbf{B}}_{(i)}^T\right)={\mathbf{B}}_i{\mathbf{A}}_{(i)}{\mathbf{B}}_i^T$$

$$\mathbf{\Phi }=\left(\begin{array}{cccc}{\mathbf{\Phi }}_{B_0}& & & \\ & {\mathbf{\Phi }}_{B_1}& & \\ & & \ddots & \\ & & & {\mathbf{\Phi }}_{B_g}\end{array}\right)$$

$$\mathbf{N}=\left(\begin{array}{ccccc}{\mathbf{L}}^T{\mathbf{\Phi }}_{B_0}& {\hat{\mathbf{A}}}_{{\mathcal{B}}_0{\mathcal{B}}_1}{\mathbf{\Phi }}_{B_1}& {\hat{\mathbf{A}}}_{{\mathcal{B}}_0{\mathcal{B}}_2}{\mathbf{\Phi }}_{B_2}& \cdots & {\hat{\mathbf{A}}}_{{\mathcal{B}}_0{\mathcal{B}}_g}{\mathbf{\Phi }}_{B_g}\end{array}\right)$$

$${\mathbf{\Phi }}_{B_0}\leftarrow {\mathbf{Q}}_1^T{\mathbf{L}}^T{\mathbf{\Phi }}_{B_0}$$

$${\mathbf{A}}_{(i-1)}\mathbf{\Phi }={\mathbf{G}}_i{\mathbf{A}}_{(i-\frac{1}{2})}{\mathbf{G}}_i^T\mathbf{\Phi }$$

$${\mathbf{\Phi }}_D\leftarrow \left(\begin{array}{c}{\mathbf{\Phi }}_{B_1}\\ {\mathbf{\Phi }}_{B_2}\\ ⋮\\ {\mathbf{\Phi }}_{B_c}\end{array}\right)$$

$${\mathbf{\Phi }}_D={\tilde{\mathbf{\Pi }}}_D(\mathcal{D},:)$$

$${\mathbf{A}}^{-1}\approx {\mathbf{A}}_{\ell }^{-1}$$

$${\mathbf{G}}_1^{-T}{\mathbf{B}}_1^{-T}{\mathbf{G}}_2^{-T}{\mathbf{B}}_3^{-T}\cdots {\mathbf{G}}_{\ell -1}^{-T}{\mathbf{B}}_{\ell -1}^{-T}{\mathbf{G}}_{\ell }^{-T}$$

$$D_{rsu}=\left((rd{h}_1,(r+1)d{h}_1)\times (sd{h}_2,(s+1)d{h}_2)\times (ud{h}_3,(u+1)d{h}_3)\right)$$

Full text

{AMS}

65F05, 65F08, 65F50, 65Y20

1 Introduction

A significant class of problems in engineering lead to large and sparse symmetric positive definite (SPD) systems:

[TABLE]

where $\mathbf{A}\in\mathbb{R}^{n\times n}$ is a sparse SPD matrix, $\mathbf{b}\in\mathbb{R}^{n}$, and $\mathbf{x}\in\mathbb{R}^{n}$ is the desired unknown solution. In particular, we are interested in the discretizations of second-order elliptic partial differential equations (PDEs) obtained using finite stencil, finite volume, or finite element methods. Examples include the Laplace, elasticity, or (some cases of) Maxwell equations. A substantial effort in scientific computing has been devoted to efficiently solve Eq. (1) arising from such PDEs.

1.1 Previous work

The most reliable method for solving Eq. (1) is the (exact) Cholesky factorization. A naive implementation has $\mathcal{O}(n^{3})$ computational cost but sparsity can be exploited to limit the fill-ins and reduce the cost. Many methods have been designed to limit the fill-ins based on appropriately ordering the variables [22, 41]. In the context of PDEs, an efficient method is nested dissection [21, 37], which can reduce the costs to $\mathcal{O}(n^{3/2})$ in 2D, and $\mathcal{O}(n^{2})$ in 3D, under mild assumptions. In fact, nested dissection is at the heart of many state-of-the-art direct solvers [2, 16, 29], which are useful for small to middle-size systems. Still, the $\mathcal{O}(n^{2})$ complexity becomes impractical for large-scale problems.

Another group of approaches are the iterative methods such as multigrid [27, 10, 54] or the Krylov-space methods. In particular, algebraic multigrid (AMG) [46, 52, 50] does not require the knowledge of the grid geometry and removes some limitations of its predecessor, geometric multigrid. AMG can often achieve the optimal $\mathcal{O}(n)$ solution cost, with good parallel scalability and low memory requirements. To ensure convergence, however, one has to properly design the smoothing and coarsening strategy. As a result, AMG often has to be fine-tuned or extended to be efficient for specific equations, or problems in question. Krylov-space methods on the other hand, such as Conjugate Gradient [30], GMRES [48], or MINRES [43], are a general group of approaches which utilize only sparse matrix-vector products. Convergence is very sensitive, however, to the conditioning of the given system, and in practice these methods have to be coupled with preconditioners. A popular class of black-box preconditioners are the ones based on incomplete factorizations such as incomplete Cholesky [35, 33, 25], or more generally, incomplete LU (ILU) [47, 19]. However, these preconditioners have rather limited applicability, not being able to target the whole spectrum of $\mathbf{A}$, and typically specialized problem-dependent preconditioners have to be developed. In fact, in the context of elliptic PDEs, the multigrid methods often prove effective preconditioners.

Another group of preconditioners are hierarchical algorithms which exploit the fact that in the context of elliptic PDEs, certain off-diagonal blocks of $\mathbf{A}$ or $\mathbf{A}^{-1}$ are numerically low-rank [13, 4, 7, 8]. Hierarchical approaches do not typically make any other assumptions about the system in question and therefore they can be more robust than multigrid methods, e.g., in the presence of non-smooth coefficients, or strong anisotropies. The theoretical framework for hierarchical approaches is provided by $\mathcal{H}$- and $\mathcal{H}^{2}$-matrices [26, 28] (developed originally for integral equations). They allow for performing algebraic operations on matrices with low-rank structure of off-diagonal blocks or well-separated blocks. In particular, the LU factorization (and applying the inverse) can be performed with linear or quasilinear complexity [24] (typically with a given accuracy). In practice, however, the constants involved in the asymptotic scalings may be somewhat large due to the recursive nature of the algorithms which also require specific data-sparse formats for storing the matrices.

Recently, new hierarchical approaches have been proposed that concentrate on efficiently factorizing the matrix using the low-rank structure of the off-diagonal blocks (we further call them hierarchical solvers). Hierarchical Interpolative Factorization [31] is a sparsified nested dissection multifrontal approach which successively reduces the sizes of the separating fronts. LoRaSp [45] is a different approach, using extended sparsification, related to the inverse fast multipole method [18] (a similar method was also described in [51]). These methods leverage the low-rank properties to factorize the matrix approximately into products of sparse block triangular factors, obviating the need for hierarchical data-sparse matrix formats. More robust extensions have been proposed since [15, 11, 20]. In particular, [11] proposed Sparsified Nested Dissection that is guaranteed to never break and can be applied to any SPD matrix.

However, hierarchical solvers still lack good convergence guarantees. For example, one can prove that—on certain elliptic equations—geometric multigrid leads to a preconditioned system with a bounded condition number. The hierarchical approaches, on the other hand, often cannot guarantee rapid convergence of iterative methods because the accuracy of the approximation $\mathbf{A}_{\ell}\approx\mathbf{A}$ is controlled in the following sense:

[TABLE]

whereas in fact, a stronger criterion is needed:

[TABLE]

In particular, this means that $\mathbf{A}_{\ell}$ needs to be particularly accurate on the eigenvectors corresponding to small eigenvalues (the near-kernel eigenspace), which is not assured by Eq. (2). However, for many elliptic PDEs, these eigenvectors are smooth and this property can be taken advantage of, to ensure Eq. (3). This would make hierarchical solvers truly competitive, by guaranteeing a bounded number of Conjugate Gradient iterations, for instance.

1.2 Contributions

We introduce a hierarchical approach to approximately factorize sparse symmetric positive definite (SPD) matrices arising from elliptic PDE discretizations. We call it Sparse Geometric Factorization, or SGF for short. Based on the factorization, we propose a family of preconditioners with varying accuracies and sparsity patterns. Depending on the chosen sparsity pattern, the preconditioners can be computed in $\mathcal{O}(n)$ or $\mathcal{O}(n\log{n})$ operations where $n$ is the number of unknowns, with $\mathcal{O}(n)$ memory requirements. The obtained cheaply-invertible operator $\mathbf{A}_{\ell}\approx\mathbf{A}$ retains the SPD property of $\mathbf{A}$, and the factorization always succeeds, with all pivots guaranteed to be invertible in exact arithmetic (in this sense, it follows [11, 15]).

SGF shares the general framework with [11, 31, 15], ensuring however that the obtained operator $\mathbf{A}_{\ell}\approx\mathbf{A}$ is also accurate on the critical near-kernel smooth eigenvectors, which addresses the limitation of hierarchical approaches mentioned above. Our choice of approximating the near-kernel subspace are the piecewise polynomial functions of the grid (e.g., piecewise constant, or piecewise linear vectors), but the space approximating the near-kernel subspace can be supplied by the user.

In the context of $\mathcal{H}$-matrix preconditioners for elliptic PDE operators, an improved accuracy on the near-kernel eigenspace can be achieved by preserving the action of $\mathbf{A}$ on piecewise constant vectors when defining the low-rank bases (see e.g. [5, 6] for details). An attempt to incorporate these ideas into hierarchical solvers was made in [55], which improved the robustness of LoRaSp [45] by preserving the action of $\mathbf{A}$ on the globally constant vector throughout computations. In our approach, sparsifying the off-diagonal blocks in the factorization is driven primarily by requiring that the action of $\mathbf{A}$ on piecewise polynomial functions be preserved.

Compared to previous work [5, 6, 55], we also investigate more options to select different types of blocks for compressions. In particular, our methods allow sparsifying matrix blocks that cannot be expected to be low-rank, still resulting in effective preconditioners. Our family of preconditioners adapts the grid partitioning methods from [11, 31], unifying them into one framework, and adding new approaches that allow for better control of the preconditioner sparsity pattern. To keep the paper focused, the partitionings are defined on cartesian grids and relatively simple stencils. However, the SGF algorithm is very general and does not assume any particular grid structure or discretization scheme, and can be applied to any SPD matrix [11].

We benchmark the preconditioners on large systems of different types. The iteration counts of preconditioned Conjugate Gradient grow very slowly with the system sizes, and the solution times scale roughly as $\mathcal{O}(n)$, also on ill-conditioned systems. We also compare our methods to equivalent approaches using rank-revealing decompositions to approximate the off-diagonal blocks, as in other hierarchical algorithms [11, 31, 15, 51, 45]. Our methods scale significantly better with growing system sizes on most tested problems. While only briefly mentioned in this paper, we believe that our preconditioners have very promising parallelization properties (inherited from [31, 45, 11]). Some possible parallelization strategies have been described in [15, 36].

1.3 Organization of the paper

An overview of the SGF algorithm is described in Section 2 along with some useful definitions. In Section 3 we describe our method to preserve piecewise polynomial vectors in the factorization. The SGF algorithm is then described in detail in Section 4. The proposed family of preconditioners, with specific realizations of the generic factorization algorithm, is described in Section 5. Some numerical analysis is provided in Section 6. Numerical experiments are presented in Section 7, followed by conclusions in Section 8.

2 Overview of Sparse Geometric Factorization (SGF)

In the PDE applications, it is often convenient to consider a partition of the set of unknowns $V=\{1,2,\ldots,n\}$ into small disjoint subsets (here called nodes) $V=B_{1}\cup B_{2}\cup\cdots\cup B_{t}$. The partition naturally induces a block representation of $\mathbf{A}$. The interaction matrix $\mathbf{A}_{B_{i}B_{j}}$ between $B_{i}$ and $B_{j}$ will typically be nonzero only if the stencils of the discretization grid involving $B_{i}$ and $B_{j}$, are close to each other, resulting in sparse $\mathbf{A}$.

2.1 Algorithm

Given the framework, we perform a sparse approximate factorization of $\mathbf{A},$ proceeding in a hierarchical fashion:

(Elimination step) Eliminate selected nodes (called interior nodes) using the block Cholesky factorization, to obtain $\mathbf{A}=\mathbf{G}\mathbf{A}_{\left(1/2\right)}\mathbf{G}^{T}$. 2. 2.

(Compression step) Sparsify selected remaining nodes to obtain $\mathbf{A}_{\left(1/2\right)}\approx\mathbf{B}\mathbf{A}_{(1)}\mathbf{B}^{T},$ where $\mathbf{A}_{(1)}$ is sparser than $\mathbf{A}_{\left(1/2\right)}$, and $\mathbf{B}$ is sparse block-diagonal. 3. 3.

Form a new coarser partition. If the new partition is composed of a single element, factorize $\mathbf{A}_{(1)}$ using exact Cholesky factorization. 4. 4.

Otherwise recurse on the new partition and the submatrix of $\mathbf{A}_{(1)}$ corresponding to the not-yet-eliminated variables.

After all $\ell=\mathcal{O}(\log{n})$ iterations have been completed (also referred to as levels), the obtained approximate operator, which is also SPD, is of the form:

[TABLE]

The factorization is illustrated in Fig. 2.1. The partitions can be based on nested dissection [21] which determines the interior nodes eliminated in Step 1 above. Elimination of an interior node introduces fill-in only between its neighboring separators. Simple domain partitioning may be used as well in which Step 1 is skipped.

The key step of the algorithm is Compression (Step 2). Without it, the factorization would become an exact block Cholesky decomposition. The role of compression is to repeatedly eliminate additional variables from each selected node, to ensure $\mathcal{O}(n)$ complexity of applying the approximate operator $\mathbf{A}_{\ell}^{-1}$. This elimination is approximate, that is, each compression introduces small error; however, the way $\mathbf{A}$ acts on piecewise polynomial functions of the grid, is preserved, thus retaining accuracy on smooth eigenvectors. Compression is described in Section 3.

2.2 Domain partition and partition graph

We denote the set of row indices of $\mathbf{A}$ (called unknowns or variables) by $V=\{1,2,\ldots,n\}$.

Definition 2.1**.**

*A partition of $V=\{1,2,\ldots,n\}$ is a finite collection $P(V)=\{B_{1},B_{2},\ldots,B_{t}\}$ of disjoint subsets $B_{i}\subseteq V$ such that $V=B_{1}\cup B_{2}\cup B_{3}\cup\cdots\cup B_{t}$. *

Definition 2.2**.**

*We say that a partition $P^{\prime}(V)$ is a coarsening of partition $P(V)$, (or a coarser partition) which we denote by $P(V)\prec P^{\prime}(V)$, if for each set $D\in P^{\prime}(V)$, there is a subcollection $\{C_{1},C_{2},\ldots,C_{s}\}\subseteq P(\Omega)$ such that $D=C_{1}\cup C_{2}\cup\cdots\cup C_{s}$. *

In the sequel, we will distinguish the variables from $B_{i}$ that have not been eliminated since the eliminated ones no longer play a role in the factorization.

Definition 2.3**.**

*Let $P_{k}(V)$ denote the partition in the level $k$ of the algorithm outlined in Section 2.1. The active subset of $B\in P_{k}(V)$, denoted by $\mathcal{B}$, is the subset of variables in $B$ that are not eliminated at the current stage of the algorithm. *

Definition 2.4**.**

*Let $P_{k}(V)=\{B_{1},B_{2},\ldots,B_{t_{k}}\}$ be the partition in the $k$-th level of the algorithm outlined in Section 2.1. The associated set of nodes is defined as $N_{k}(V)=\{(B_{1},\mathcal{B}_{1}),(B_{2},\mathcal{B}_{2}),\ldots,(B_{t_{k}},\mathcal{B}_{t_{k}})\},$ where $\mathcal{B}_{i}$ is the active subset of $B_{i}$. *

We assume that $\mathbf{A}_{\mathcal{B}_{i}\mathcal{B}_{j}}$ (the interaction matrix between two nodes) can be nonzero only if grid stencils corresponding to some elements of $B_{i}$ and $B_{j}$, are close to each other. The number of variables $|\mathcal{B}|$ will be referred to as the size of the node.

3 Key idea: polynomial compression

The key difference between our approach and similar methods performing sparse approximate factorizations ([45, 31, 11, 51, 40, 14]), is the design of compression step. The idea of compression is to drop interactions of many variables while retaining accuracy on the near-kernel eigenspace. In the case of second-order elliptic differential operators, it is known that the eigenvectors corresponding to smallest eigenvalues are smooth. We will therefore assume that an eigenvector with small associated eigenvalue can be locally approximated well by a discretized low-degree polynomial of the grid (see below). More generally, the subspace approximating the near-kernel eigenspace can be supplied by the user.

3.1 Basis of discretized polynomials

We denote by $\mathbf{\Pi}^{0}$ the $n\times 1$ matrix (a vector of length $n$) of ones. If every unknown has a specified location in $\mathbb{R}^{3}$, e.g., the location of the underlying grid vertex, let $(x_{i},y_{i},z_{i})\in\mathbb{R}^{3}$ denote the location of the $i$-th unknown, for $i\in V=\{1,2,3,\ldots,n\}$. If that be the case, we define:

[TABLE]

In other words, the columns of the matrix $\mathbf{\Pi}^{j}$ for $j=0,1,2$, are a basis of the space of discretized real polynomials of degree $j$. Once the degree has been chosen, we write $\mathbf{\Pi}:=\mathbf{\Pi}^{j}$. We denote the number of columns of $\mathbf{\Pi}$ by $\pi$. The definition of $\mathbf{\mathbf{\Pi}}$ may need to be modified when more than one variable has the same underlying location, for example in the case of a vector PDE (see Section 7.1.3).

3.2 Compression step

Given an SPD matrix $\mathbf{M}\in\mathbb{R}^{n\times n}$ and a node $(B,\mathcal{B})\in N(V)$, let $\mathcal{N}$ denote the set of variables interacting with $\mathcal{B}$ and let $\mathcal{W}$ denote the set of all the remaining variables. We have:

[TABLE]

Compressing interactions of $(B,\mathcal{B})$ will not modify the interactions of $\mathcal{W}$ and so to keep compact notation, we consider only the upper left $2\times 2$ block. Also, for simplicity, assume that there are only two nodes having nonzero interactions with $(B,\mathcal{B})$, and denote them by $(N_{1},\mathcal{N}_{1}),$ and $(N_{2},\mathcal{N}_{2})$, respectively.

First, we scale the block row and column corresponding to $(B,\mathcal{B})$:

[TABLE]

where $\mathbf{M}_{\mathcal{B}\mathcal{B}}=\mathbf{L}\mathbf{L}^{T}$ denotes the Cholesky decomposition, and $\mathbf{\widehat{M}}_{\mathcal{B}\mathcal{N}}=\mathbf{L}^{-1}\mathbf{M}_{\mathcal{B}\mathcal{N}}.$ Now, let $\mathbf{\Phi}_{B}=\mathbf{\Pi}(\mathcal{B},:),$ and likewise for $(N_{1},\mathcal{N}_{1})$ and $(N_{2},\mathcal{N}_{2})$. Then

[TABLE]

spans the space of piecewise (discretized) polynomial functions. Consider the $|\mathcal{B}|\times 3\pi$ matrix (which we call the filtered interaction matrix):

[TABLE]

and an orthogonal basis $\mathbf{Q}_{1}$ for its range, e.g., from the column-pivoted QR decomposition:

[TABLE]

The matrix $\mathbf{N}$ will typically have a small number of columns, and therefore $\mathbf{Q}_{1}$, spanning its range, is thin. Notice also that $\mathbf{Q}_{2}^{T}\mathbf{\widehat{M}}_{\mathcal{B}\mathcal{N}_{1}}\mathbf{\Phi}_{N_{1}}$, $\mathbf{Q}_{2}^{T}\mathbf{\widehat{M}}_{\mathcal{B}\mathcal{N}_{2}}\mathbf{\Phi}_{N_{2}},$ and $\mathbf{Q}_{2}^{T}\mathbf{L}^{T}\mathbf{\Phi}_{B}$ are zero matrices because $\mathbf{Q}_{2}$ spans the space orthogonal to the range of $\mathbf{N}$. With

[TABLE]

we therefore have:

[TABLE]

Based on the equation above, we can drop the blocks $\mathbf{Q}_{2}^{T}\mathbf{\widehat{M}}_{\mathcal{B}\mathcal{N}_{2}}$ and $\mathbf{Q}_{2}^{T}\mathbf{\widehat{M}}_{\mathcal{B}\mathcal{N}_{2}}$ (and their transposes) when approximating $\mathbf{M}$ as they do not affect its action on piecewise polynomial vectors. In other words, we approximate:

[TABLE]

A number of variables from $\mathcal{B}$ get eliminated from the system as they no longer interact with other variables in $\mathbf{M}_{(+)}$, which is strictly sparser than $\mathbf{M}$. The smaller the rank of $\mathbf{Q}_{1},$ the larger the number of eliminated variables. We then perform a cheap update:

[TABLE]

and we can recurse on $\mathbf{M}_{(+)}$ to sparsify interactions of the subsequent node. After interactions of $(N_{1},\mathcal{N}_{1})$ and $(N_{2},\mathcal{N}_{2})$ have also been compressed, subsets of variables from each of $\mathcal{B}$, $\mathcal{N}_{1}$, and $\mathcal{N}_{2}$ are eliminated, and we have $\mathbf{M}\approx\mathbf{U}\mathbf{M}_{(++)}\mathbf{U}^{T}$ where $\mathbf{M}_{(++)}$ is much sparser than $\mathbf{M},$ and $\mathbf{U}$ is a sparse block diagonal matrix.

Going back to the general case, assuming that there are $g$ nodes interacting with $(B,\mathcal{B}),$ compressing its interactions would be written as $\mathbf{M}\approx\mathbf{\widetilde{M}},$ where:

[TABLE]

The space unaffected by compression is formally described in the following.

Lemma 3.1**.**

For an SPD matrix $\mathbf{M}$, let $\mathbf{\widetilde{M}}$ be defined as in Eq. (9). Then for any $B\in P(V),$ we have:

[TABLE]

where $\mathbf{\Pi}_{B}$ is the $|B|\times\pi$ matrix defined by:

[TABLE]

*Moreover, $\mathbf{\widetilde{M}}$ is also SPD. *

Proof 3.2**.**

*It remains to prove the preservation of the SPD property. But $\mathbf{\widetilde{M}}$ is SPD iff $\mathbf{M}_{(+)}$ is SPD, where $\mathbf{\widetilde{M}}=\mathbf{B}\mathbf{M}_{(+)}\mathbf{B}^{T}$ as above. But $\mathbf{M}_{(+)}$, up to a permutation of variables, is composed of two non-interacting diagonal blocks of which one is the identity matrix and the other one is, again up to a permutation of variables, a principal submatrix of $\mathbf{B}^{-1}\mathbf{M}\mathbf{B}^{-T}$, which is SPD. *

Note a minor technical detail. In Lemma 3.1, we use $B$ instead of $\mathcal{B}$ (i.e., the preserved vectors correspond to all unknowns in $B$). Initially, one would have $B=\mathcal{B}$, i.e., when no variables have been eliminated, one puts $\mathbf{\Phi}_{B}=\mathbf{\Pi}(\mathcal{B},:)=\mathbf{\Pi}(B,:)$. However, later on, when variables get eliminated, the globally preserved vectors are still given by Eq. (11), as long as the matrices $\mathbf{\Phi}_{B}$ are properly updated.

Notice also that the scaling by $\mathbf{L}$ before compression Eq. (5) plays an important role in the proof of the SPD property of $\mathbf{\widetilde{M}}$. In particular, the factorization will never fail (at least in exact precision). Similar observation has been made in [11, 14, 53].

3.3 Using low-rank approximation

We would like to note here that our methods do not preclude the use of low-rank approximations in compression. Namely, provided that the matrix $\mathbf{M}_{\mathcal{B}\mathcal{N}}$ is numerically low-rank, the matrix $\mathbf{Q}_{1}$ can be ensured to also satisfy $\mathbf{\widehat{M}}_{\mathcal{B}\mathcal{N}}=\mathbf{Q}_{1}\mathbf{Y}^{T}_{1}+\mathbf{Q}_{2}\mathbf{Y}^{T}_{2}$, where $\|\mathbf{Y}_{2}^{T}\|_{2}$ is small. In that case, the error $\|\mathbf{M}-\mathbf{\widetilde{M}}\|_{2}$ is controlled more directly. In the absence of geometrical information, for instance, the constant vector, i.e., $\mathbf{\Pi}=\mathbf{\Pi}^{0}$ can be used, combined with a rank-revealing decomposition, to ensure that $\|\mathbf{Y}_{2}^{T}\|_{2}$ above is small.

4 Sparse Geometric Factorization (SGF)

Denote by $\mathbf{A}_{k}$ the approximation to $\mathbf{A}$ after $k$ levels of the algorithm outlined in Section 2. Eliminating many nodes before any compressions occur, can guarantee that (for any $k$):

[TABLE]

where $\mathbf{Y}$ is a matrix with $\Theta(n)$ orthogonal columns. However, we would particularly like Eq. (12) to also hold for $\mathbf{Y}$ whose columns span the eigenspace of $\mathbf{A}$ corresponding to $\Theta(n)$ smallest eigenvalues. Then we could expect $\mathbf{A}^{-1}_{k}$ to be a very close approximation of $\mathbf{A}^{-1}$ (see Section 6). However, we do not know the eigenvectors a priori, and we want $\mathbf{A}_{k}^{-1}$ to be as sparse as possible. This will motivate the algorithm, which we now describe in detail.

We denote by $P_{i}(V)$ the partition at level $i$ of the algorithm, and by $N_{i}(V)$ the set of the corresponding nodes. The nodes need to be divided into two subsets: the set of interiors, denoted by $I_{i}$ and the set of separators, denoted by $S_{i}$. An interior node can only interact with separator nodes; interior nodes will be eliminated by Gaussian elimination. However, the interactions of selected separators will also be compressed (which does not violate their separating properties). With each node $(B,\mathcal{B})$, we associate a matrix $\mathbf{\Phi}_{B},$ where in the first partition we put $\mathbf{\Phi}_{B}\leftarrow\mathbf{\Pi}_{B}(B,:)=\mathbf{\Pi}(B,:)$. We denote the total number of levels by $\ell$. We also put $\mathbf{A}_{(0)}:=\mathbf{A}$.

4.1 Eliminating interior nodes

We eliminate the interior nodes, if any, using the standard block Cholesky algorithm, which at the $i$-th level we denote by:

[TABLE]

4.2 Compressing selected separators

To sparsify $\mathbf{A}_{(i-\frac{1}{2})}$ while preserving some form of Eq. (12), we apply the polynomial compression (Section 3.2) for each selected node, which we denote by:

[TABLE]

We now describe this step and explain its correctness. From Section 3.2 we know that compressing interactions of a single node $(B_{0},\mathcal{B}_{0})$ with $g$ interacting nodes, leads to Eq. (12) for $k=i$, with $\mathbf{Y}=\mathbf{\Phi}$, which is our space approximating the near-kernel eigenspace:

[TABLE]

The filtered interaction matrix from Eq. (6) whose QR we need to compute, reads:

[TABLE]

This means that compression is practical only if the blocks $\mathbf{\Phi}_{B}$ have small numbers of columns. Otherwise no or very few interactions can be compressed. After compressing interactions of $(B_{0},\mathcal{B}_{0})$, we need to update (as in Eq. (8)):

[TABLE]

The update Eq. (17) does not change the block structure of $\mathbf{\Phi}$, and so compressions in the given level can continue in the same manner.

Notice, however, that when eliminating interior nodes, it seems we also need to update $\mathbf{\Phi},$ since then:

[TABLE]

The update $\mathbf{\Phi}\leftarrow\mathbf{G}_{i}^{T}\mathbf{\Phi}$ would ruin the block-diagonal-like structure of Eq. (15). Fortunately, we do not need to perform this update at all because it only changes the rows of $\mathbf{\Phi}$ that correspond to variables just eliminated by block Cholesky.

4.3 Constructing the new partition

After compressions in level $i$ have been completed, we need to define the new partition $P_{i+1}(V)\succ P_{i}(V)$ (see Def. 2.2). Every $B\in P_{i}(V)$ must have a unique “father” $D\in P_{i+1}(V)$ such that $B\subseteq D$.

4.4 Forming the new $\mathbf{\Phi}$ matrix

With $\mathbf{\Phi}$ defined as in Eq. (15) (with blocks defined by partition in the first level), when coarser partitions are formed, the numbers of columns of nonzero sub-blocks of $\mathbf{\Phi}$ corresponding to the new partition sets will also grow, and eventually the compression will become impractical. On the other hand, we can keep modifying $\mathbf{\Phi}$ so that for each set $D$ in the $k$-th partition, the matrix $\mathbf{Y}_{D}$ spanning the space of discretized polynomials over $D$ (of a pre-chosen degree equal for each level), satisfies Eq. (12). Notice that if this is in fact true for level $k$, then it is also true at the beginning of level $k+1,$ because of the nature of polynomials, and the fact that each set $D$ at level $k+1$ is a union of sets from level $k$.

However, we still need to make sure that the appropriate blocks in Eq. (15) correspond to polynomial bases of the sets at the new level. When a set $D$ is formed by merging the sets $B_{1},B_{2},\ldots,B_{c}$ in the algorithm, then its active subset $\mathcal{D}=\mathcal{B}_{1}\cup\mathcal{B}_{2}\cup\cdots\cup\mathcal{B}_{c}$, where $\mathcal{B}_{j}$ is the active subset of $\mathcal{B}_{j}$. Again in light of the nature of updates Eq. (17), this means that we can simply update:

[TABLE]

By definition, the number of columns of matrix $\mathbf{\Phi}_{D}$, for any set $D$ at any level, is bounded at all times. This is crucial for achieving $\mathcal{O}(n)$ or $\mathcal{O}(n\log{n})$ complexity of factorization, and $\mathcal{O}(n)$ complexity of applying $\mathbf{A}_{\ell}^{-1}$ (see Section 5). We have:

[TABLE]

where $\mathbf{\tilde{\Pi}}_{D}=\mathbf{B}_{i}^{T}\mathbf{G}_{i}^{T}\cdots\mathbf{B}_{2}^{T}\mathbf{G}_{2}^{T}\mathbf{B}_{1}^{T}\mathbf{G}_{1}^{T}\mathbf{\Pi}_{D}$, for $\mathbf{\Pi}_{D}$ as in Eq. (11). This means that $\mathbf{\Phi}_{D}$ is exactly the matrix needed in the compression at level $i+1$ of the algorithm, that is, it corresponds to the space of discretized polynomials over $D$.

4.5 Eliminating the last node

If the new partition is all of $V$, i.e., $P_{i}(V)=V$, we have reached the top level of the algorithm, i.e., $i=\ell$. Eliminating $(V,\mathcal{V})\in N_{\ell}(V)$ via exact Cholesky decomposition concludes the factorization.

4.6 Pseudocode of the factorization

The hierarchical factorization is summarized in the pseudocode below.

4.7 Applying the preconditioner

The approximate inverse operator reads:

[TABLE]

From Eq. (13) and Eq. (14) it should be clear that $\mathbf{A}_{\ell}^{-1}$ can be applied by a sequence of interspersed block triangular solves and multiplications by block-diagonal matrices.

5 Family of preconditioners

In this section, we propose a family of preconditioners composed of four realizations of the SGF algorithm. To illustrate the ideas, we consider a cartesian grid with a relatively simple discretization scheme. We would like to stress, however, that our methods are not limited to cartesian grids. There is no assumption about the grid structure in Section 4; the algorithm only requires a partition hierarchy and the $\mathbf{\Pi}$ matrix. The appropriate partitions can be defined geometrically given a general grid, or purely using the matrix graph. In fact, our methods are being incorporated into spaND [11] which uses METIS [34] to obtain appropriate graph partitionings. To keep focus, however, in this paper we limit ourselves to cartesian grids where partitions can be explicitly defined

We therefore consider the case in which the discretized domain is of the form $\Omega=(0,a_{1})\times(0,a_{2})\times(0,a_{3})$ and the grid vertices are located at points $(\frac{kh_{1}}{2},\frac{lh_{2}}{2},\frac{mh_{3}}{2})$ where $h_{1},h_{2},h_{3}>0$ are small and fixed, and $k,l,m$ are natural numbers. We assume that the discretization is such that the matrix entry $\mathbf{A}_{ij}$ can be nonzero only if the grid vertex corresponding to the $j$-th unknown is in the $3\times 3\times 3$ subgrid of vertices around the vertex corresponding to the $i$-th unknown.

5.1 Preconditioners using nested dissection partitions

We consider the partitions $P_{1}(V)\prec P_{2}(V)\prec\cdots\prec P_{\ell}(V)$ inspired by nested dissection [21, 49]. We choose a natural number $b>1$. The sets in $P_{t}(V)$ are formed when the domain $\Omega$ is cut by the planes: $\{x=(rd-1)h_{1}\},$ $\{x=rdh_{1}\},$ $\{y=(sd-1)h_{2}\},$ $\{y=sdh_{2}\},$ $\{z=(ud-1)h_{3}\},$ $\{z=udh_{3}\}$ for $r,s,u\in\{1,2,3,\ldots\}$, where $d=2^{t-1}\,b$. Two variables $i,j\in V$ then belong to the same partition set iff their corresponding grid vertices are in the same set formed by the cutting planes. This is illustrated in Fig. 5.1.

Notice that there are three types of sets formed by the cutting planes, which we call 0-, 1-, and 2-cells. The 0-cell is a cube and contains a single grid vertex (the 0-cells touch the corners of eight large cubes). The 1-cell is typically a longitudinal parallelpiped containing $d-1$ grid vertices (1-cells touch the edges of four large cubes). The 2-cell is typically a flat parallelpiped containing $(d-1)^{2}$ grid vertices (the 2-cells touch the faces of two cubes). Finally, a 3-cell is typically a cube containing $(d-1)^{3}$ grid vertices (these are the large cubes in Fig. 5.1). The 0- 1- and 2-cells naturally create a “buffer” between the 3-cells and the rest of the domain. In a 2D case (for example when the domain is $(0,1)^{2}$) there would be no 3-cells. In this case the 0- and 1-cells naturally separate the 2-dimensional cells from the rest of the domain. See Fig. 2.1 for an illustration with active unknowns visible.

The set of 3-cells will form interiors and will be eliminated using the block Cholesky algorithm (Section 4.1). Compression of the remaining cells then defines the preconditioner (here, we do not distinguish between a cell and a partition node):

1. Nest-All-All

All interactions of all cells are compressed 2. Nest-2-All

All interactions of 2-cells are compressed but others are not 3. Nest-2-2

Only the interactions of 2-cells with non-adjacent cells are compressed

Compressing only selected interactions in Nest-2-2 above is achieved by including in $\mathbf{Q}_{1}$ the span of the off-diagonal blocks connecting the given 2-cell with its adjacent 0- and 1-cells (so that most or all of the compressed interactions are those between 2-cells). Most other hierarchical algorithms similar to Nest-2-2—which compress only well-separated interactions—cannot guarantee the SPD property, and may break (for example LoRaSp [45]). In Nest-2-2, the compression ranks (numbers of columns of $\mathbf{Q}_{1}$) increase with each level (in contrast, for example, to [45, 51] where the ranks can be bounded) but the SPD property is preserved.

Nest-2-All is analogous to the schemes used in the context of Hierarchical Interpolative Factorization and Sparsified Nested Dissection [11, 31] (the latter defines the partitions based on the locations of unknowns, or purely based on the matrix entries). These approaches base the compression step, however, on low-rank approximations, not polynomial compression ([31] being also quite different algebraically).

Nest-All-All, on the other hand, is quite different in that it also compresses interactions of 1-cells (compressing interactions of 0-cells does not make any difference and can be skipped). The off-diagonal block corresponding to a 1-cell in the matrix should not be expected to be low-rank. However, polynomial compression does not require the low-rank property. Intuitively speaking, it replaces each interaction between two nodes by a coarser one, preserving the action of the matrix on a chosen subspace (as opposed to exploiting the low-rank structure to represent the interaction accurately).

5.2 Preconditioners using general partitions

The second type of domain partitions considered does not distinguish interior nodes that could be eliminated without introducing large fill-ins. Again, such partitions can be obtained by analyzing just the matrix entries [45], using tools such as SCOTCH [44]. As an example, with the same setup as before, $P_{t}(V)$ could be a partition formed when the domain $\Omega$ is cut into cells containing the same, or roughly the same number of grid vertices:

[TABLE]

This is shown in Fig. 5.2.

There are no interior nodes as before. Instead, at every level, the size of each node is reduced in compression, and then nodes are merged when forming the new partition for the subsequent level. This means that compression is the only step that includes Cholesky factorization (the scaling Eq. (5)). We obtain:

1. Gen-All-All

All interactions of all cells are compressed

One should not expect approaches that use low-rank approximations in compression to perform well with Gen-All-All unless the matrix satisfies the strong admissibility criterion [9]. The approaches that do use simple partitioning ([45, 51, 40, 14]) typically distinguish between neighboring and well-separated interactions, the latter assumed to be low-rank. Only the well-separated interactions are compressed. The methods of this paper can also be used in that context. The advantage of compressing all interactions of a node, however, is that we can obtain a much sparser operator, and we ensure that $\mathbf{A}_{\ell}$ is SPD.

5.3 Accuracy comparison

In terms of accuracy of factorization, one can conceptually write:

[TABLE]

meaning that Gen-All-All is the least accurate (no interiors eliminated, all interactions compressed), and Nest-2-2 is the most accurate (interiors eliminated; only interactions of 2-cells with non-adjacent cells compressed).

5.4 Complexities of factorization and applying the approximate inverse

As we describe next, Nest-All-All and Gen-All-All can be expected to have $\mathcal{O}(n)$ complexity. This is an advantage over Nest-2-All or Nest-2-2 used in [11, 31] which generally have $\mathcal{O}(n\log{n})$ complexity. The partitions used in our methods satisfy the assumptions of the proposition below, with $\delta=8$. The assumptions common to both considered cases should be satisfied by any reasonably balanced partitions. Similar (including more general) complexity analyses are described in [11, 45, 55].

Proposition 5.1**.**

Assume that the polynomial compression is used with a fixed pre-chosen polynomial degree, and that the node sizes in the first level are bounded. Also assume that any given node in the algorithm has nonzero interactions with a bounded number of other nodes, and that a set $B\in P_{k}(V)$, is a union of a bounded number of sets from $B\in P_{k-1}(V)$. For $\delta>1$, assume that there exists a constant $C_{0}>0$ such that the partition $P_{k}(V)$ is composed of at most $C_{0}n/\delta^{k-1}$ sets. Then:

For a scheme compressing all interactions (such as Nest-All-All and Gen-All-All) the complexities and memory requirements of computing the factorization, as well as applying $\mathbf{A}_{\ell}^{-1}$, are all $\mathcal{O}(n)$. 2. 2.

For any other scheme, if there exists a constant $C_{1}>0$ such that the size of any node at level $k$ is at most $C_{1}\delta^{(k-1)/3}$, then the complexity of computing the factorization is $\mathcal{O}(n\log{n})$. Memory requirements as well as the complexity of applying $\mathbf{A}_{\ell}^{-1}$ are $\mathcal{O}(n)$. This holds for Nest-2-All and Nest-2-2 with $\delta=8$.

Proof 5.2**.**

Consider first case 1. Since the polynomial compression with fixed polynomial degree is used, and every node has a nonzero interaction with $\mathcal{O}(1)$ other nodes, the ranks of the filtered interaction matrices Eq. (16) are likewise $\mathcal{O}(1)$, and therefore the sizes of nodes are also bounded. The $\mathcal{O}(n)$ factorization complexity follows from the fact that $\sum_{k=0}^{\infty}{\frac{n}{\delta^{k}}}=\mathcal{O}(n)$, and similarly for memory requirements.

*Consider now case 2. The cost associated with eliminating or compressing a node this time is $\mathcal{O}(\delta^{k})$. Hence we obtain the complexity bound from $\sum_{k=0}^{\lceil\log_{\delta}C_{0}n\rceil}{\frac{n}{\delta^{k}}\delta^{k}}=$ $\mathcal{O}(n\log{n})$. Since we now need $\mathcal{O}(\delta^{2k/3})$ memory to store a node and its interactions, we obtain the memory requirement by noticing that $\sum_{k=0}^{\lceil\log_{\delta}C_{0}n\rceil}{\frac{n}{\delta^{k}}\delta^{2k/3}}=\sum_{k=0}^{\lceil\log_{\delta}C_{0}n\rceil}{n\delta^{-k/3}}=\mathcal{O}(n)$. The same argument applies to show that applying $\mathbf{A}_{\ell}^{-1}$ is $\mathcal{O}(n)$. To see that case 2 above holds for Nest-2-All and Nest-2-2, notice that the number of grid vertices inside a 1-cell is $\mathcal{O}(2^{k}).$ *

6 Approximation error and preconditioner quality

A successful preconditioner for SPD systems must be accurate on the eigenvectors corresponding to the near-kernel (smallest) eigenvalues. The results in [55] suggest that for the (unit-length) eigenvector $\mathbf{v}_{i}$ corresponding to the given eigenvalue $\lambda_{i}>0$, the contribution of the error $\|(\mathbf{A}-\mathbf{A}_{\ell})\mathbf{v}_{i}\|_{2},$ to the condition number of the preconditioned system, is amplified by $\lambda_{i}^{-1}$. Results in [5] show that, to achieve a condition number independent of the problem size, the accuracy of $\mathbf{A}_{\ell}$ can be relatively crude provided that Eq. (12) holds for $\mathbf{Y}$ with columns approximating well the near-kernel eigenspace.

We now describe how preserving piecewise polynomial vectors in our algorithm leads to a better preconditioner. Let $\mathbf{A}_{i}$ again denote the approximation to $\mathbf{A}$ obtained after completing the $i$-th level. For $i\in\{1,2,\ldots,\ell\}$, define the error term: $\mathbf{E}_{i}:=\mathbf{A}_{i}-\mathbf{A}_{i-1},$ where by definition $\mathbf{A}_{0}:=\mathbf{A}$. We have:

[TABLE]

The error terms $\mathbf{E}_{i}$ are the results of compression (gaussian elimination does not introduce errors in exact precision). The effect of applying polynomial compression is described in the following.

Proposition 6.1**.**

Let $k\in\{1,2,\ldots,\ell\}$ and $B\in P_{k}(V)$, where $P_{k}(V)$ is the partition at level $k$ of the algorithm. Also, let $\mathbf{\Pi}$ denote the global polynomial basis matrix described in Section 3.1. Then, for any $1\leq i\leq k,$

[TABLE]

where $\mathbf{\Pi}_{B}$ is the $|B|\times\pi$ matrix defined by:

[TABLE]

In particular, the action of $\mathbf{A}$ on $\mathbf{\Pi}$ is preserved:

[TABLE]

Proof 6.2**.**

*The proposition is true for $k=1$ which follows from Lemma 3.1. Now suppose that the proposition is true for a given $1\leq k\leq\ell-1$. Then for $B\in P_{k+1}(V)$ and $1\leq i\leq k$ we have $\mathbf{E}_{i}\mathbf{\Pi}_{B}=0$. This is a consequence of the fact that $B$ is a union of sets from $P_{k}(V)$ (see Eq. (11)). To conclude the proof, we need to show that $\mathbf{E}_{k+1}\mathbf{\Pi}_{B}=0$ but this is a consequence of Eq. (19) and the remarks that follow it. *

Denoting by $\mathbf{u}_{k}$ the piecewise polynomial approximation to $\mathbf{u}\in\mathbb{R}^{n}$ at the $k$-th level of the algorithm, we can expect that $\|\mathbf{u}-\mathbf{u}_{k}\|_{2}$ will be small if $\mathbf{u}$ is a smooth function of the grid. The approximation error on $\mathbf{u}$ at the $k$-th level will also be then small because $\mathbf{E}_{k}\mathbf{u}=\mathbf{E}_{k}(\mathbf{u}-\mathbf{u}_{k})$. The following corollary, which is a direct consequence of Proposition 6.1, makes this statement precise.

Corollary 6.3**.**

For $\mathbf{u}\in\mathbb{R}^{n}$ and $k\in\{1,2,\ldots,\ell\}$, let $\mathbf{u}_{k}$ denote the projection of $\mathbf{u}$ onto $\Upsilon_{k}:=\text{span}\{\mathbf{\Pi}_{B}\}_{B\in P_{k}(\Omega)}$. We then have:

[TABLE]

Put differently, if $\mathbf{P}_{k}$ denotes the orthogonal projection matrix onto $\Upsilon_{k}^{\perp}$, then:

[TABLE]

Based on Corollary 6.3, we can write:

[TABLE]

where for the eigenvalue decomposition $\mathbf{A}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{T}$ we write $\mathbf{A}^{\frac{1}{2}}=\mathbf{V}\mathbf{\Lambda^{\frac{1}{2}}}\mathbf{V}^{T}$. We would like the matrix Section 6 to be as close to identity as possible. In our case, we also know that it is SPD (see Lemma 3.1). Notice that:

[TABLE]

The norm $\|\mathbf{E}_{k}\|_{2}$ can be made small using a low-rank approximation as described in Section 3.3. However, if $\mathbf{A}$ is ill-conditioned, then $\mathbf{A}^{-1}$ will have large eigenvalues, and the norm on the left hand side of Eq. (23) may become large even if $\|\mathbf{E}_{k}\|_{2}$ is relatively small. Ensuring that the range of $\mathbf{P}_{k}$ approximates well the range of the eigenvectors of $\mathbf{A}^{-1}$ with small associated eigenvalues, we can therefore directly target the critical eigenspace of ill-conditioned systems, and make small the norm $\|\mathbf{P}_{k}\mathbf{A}^{-1}\mathbf{P}_{k}\|_{2}$. In particular (see also [6], Lemma 2.1), if:

[TABLE]

then from Weyl’s inequality for eigenvalues we obtain for the condition number:

[TABLE]

This suggests that to obtain a bounded condition number of the preconditioned system, one needs the spaces $\Upsilon_{k}=\text{span}\{\mathbf{\Pi}_{B}\}_{B\in P_{k}(\Omega)}$ to approximate well the eigenspace corresponding to the smalles eigenvalues of $\mathbf{A}$.

7 Numerical results

We apply our methods to problems of importance in engineering. Our goals in this section are:

To test the efficiency of our methods when applied to large systems. Ideally, we would like to observe the optimal $\mathcal{O}(n)$ scalings of solution times. 2. 2.

To compare our approach to the one in which compression is based on the standard low-rank approximation of the off-diagonal blocks. 3. 3.

To benchmark the different preconditioners described in Section 5, with varying orders of polynomial compression (piecewise constant, piecewise linear, and piecewise quadratic; that is, when $j$ in the definition of $\mathbf{\Pi}^{j}$ from Section 3.1, is $j=0,1,2$, respectively).

When testing Nest-2-All, Nest-All-All, and Nest-2-2, we use $b=3$, i.e., the 3-cell at level $k$ is typically a cube encompassing $(3\cdot 2^{k-1}-1)^{3}$ grid vertices (see Section 5); we also skip the compression in the first two levels, while the nodes are still small.

When testing Gen-All-All we choose $b=3$, $b=4$, and $b=5$ when using respectively, piecewise constant, piecewise linear, and piecewise quadratic compression. In other words, the partition sets at level $k$ typically contain variables corresponding to, respectively, $(3\cdot 2^{k-1})^{3}$, $(4\cdot 2^{k-1})^{3}$, and $(5\cdot 2^{k-1})^{3}$ grid vertices. We use compression starting from the first level.

In each test, when the number of not-yet-eliminated unknowns during factorization is small enough (smaller than the size of the largest node encountered so far), we factorize the remaining system exactly using block Cholesky decomposition (i.e., form the last partition). This does not change the maximal node size.

We solve the equation $\mathbf{A}\mathbf{x}=\mathbf{b}$ using Conjugate Gradient (CG) with operator $\mathbf{A}_{\ell}^{-1}$ used as a preconditioner. Each time we choose a random right hand side $\mathbf{b}$ to ensure that every eigenvector of $\mathbf{A}$ contributes to $\mathbf{b}$.

Throughout this section, we use the following notation:

• •

$n$ is the number of unknowns in the analyzed system;

• •

$it_{C}$ is the number of iterations of Conjugate Gradient needed to converge to a relative residual of 2-norm below $10^{-10}$;

• •

$t_{F}$ is the factorization time, i.e., the CPU time taken to compute $\mathbf{A}_{\ell}^{-1}$, in seconds;

• •

$t_{S}$ is the solution time, i.e., the CPU time taken by CG to converge, in seconds;

• •

$m_{R}$ is the maximum memory usage during the computation, in GB.

Our implementation is sequential and was written in Python 3.6.1., exploiting NumPy 1.14.3 and SciPy 1.1.0 for numerical computations [42, 32]. The tests were run on CPUs with Intel(R) Xeon(R) E5-2640v4 (2.4GHz) with up to 1024 GB RAM. All pivoted QR decompositions were performed by calling LAPACK’s geqp3 [3].

7.1 Descriptions of test cases

7.1.1 3D Poisson equation

The first test case is the classical 3D (constant-coefficient) Poisson equation in a cube:

[TABLE]

with Dirichlet boundary conditions, discretized using the standard 7-point stencil method. The largest system has approximately $16\cdot 10^{6}$ unknowns.

7.1.2 Incompressible flow in the SPE10 Reservoir

The second test case is the 3D flow equation (Darcy’s law) of an incompressible single-phase fluid, in an incompressible porous medium:

[TABLE]

discretized using the finite volume method (the 2-point flux approximation), with mixed Dirichlet and Neumann boundary conditions.

To define $\Omega$ and the field of coefficients $\lambda$ (the mobility field), we use the SPE10 Reservoir [17] which is an important benchmark problem in the petroleum engineering community. The field $\lambda$ varies smoothly in some layers of the reservoir, and is highly discontinuous in other layers (see Fig. 7.2). The result is a Poisson-like equation whose corresponding discretized system is very ill-conditioned.

The smallest test case has approximately $0.4\cdot 10^{6}$ variables and is obtained by considering only the upper layers of the SPE10 Reservoir, where $\lambda$ changes relatively smoothly. The case with approximately $1.1\cdot 10^{6}$ variables is obtained from the full original SPE10 Reservoir (with grid dimensions $220\times 80\times 65$). Larger cases are obtained by periodically tiling the original reservoir in each direction to obtain cubes of the desired size (based on idea from [39]). The matrices are obtained by fixing the pressure value on one outer side of the reservoir, and imposing a constant flow on the opposite side, with no-flow conditions on other sides, but the resulting system $\mathbf{A}\mathbf{x}=\mathbf{b}$ in each case is tested with a random right hand side, as described above.

7.1.3 Linear elasticity

The third test case is a finite-element approximation to the weak form of the linear elasticity equation:

[TABLE]

where $u:\Omega\rightarrow\mathbb{R}^{3}$ is the displacement field and $\sigma$ is the stress tensor satisfying:

[TABLE]

with $\lambda$ and $\mu$ denoting the material Lamé constants.

Our test problem is a cantilever beam composed of two segments, as in Fig. 7.2. The constants corresponding to the left hand side, and the right hand side segments are dimensionless $(\lambda_{\text{LHS}},\mu_{\text{LHS}})=(1.0,1.0)$ and $(\lambda_{\text{RHS}},\mu_{\text{RHS}}=50.0,50.0),$ respectively. The boundary conditions are a fixed zero displacement on the left side of the boundary, i.e., $u=0$ there, and a vertical constant pull down force applied on the other side. The test case (including the discretization using a regular grid with tetrahedral elements) is obtained from the MFEM library [1] where the reader may find all details. We test systems of increasing sizes by uniformly refining the grid in each dimension; the largest case consists of approximately $6.5\cdot 10^{6}$ variables.

This test differs from the previous ones in that each grid vertex $v$ is now associated with a displacement vector $\mathbf{u}_{v}=(v_{x},v_{y},v_{z})$ $\in\mathbb{R}^{3}$ (which is a subvector of the global solution vector $\mathbf{u}$). As mentioned in Section 3.1, the definition of the polynomial basis has to be modified. Assuming that the variable indices are ordered so that the indices corresponding to all the $v_{x}$ above (of all grid vertices) come first, then the indices corresponding to all the $v_{y}$ and then all the $v_{z}$ (in each case retaining the same order of the underlying grid vertices), we naturally define:

[TABLE]

where $j$ is the chosen polynomial degree, and $\mathbf{\Pi}^{j}$ is defined as in Section 3.1. Conceptually, the displacement component in each direction is then independent of the displacement components in other directions.

7.2 Preserving the rigid body modes

The elasticity test Eq. (27) is a loosely constrained body. In such a case, the so called rigid body modes are known to be in the near-kernel subspace (see for instance [38]). With $\mathbf{\Pi}$ defined as above (Eq. (28)), $\mathbf{A}_{\ell}$ preserves the action of $\mathbf{A}$ on the rigid body modes (more precisely, they are preserved exactly when using at least piecewise linear compression, which follows from Eq. (21) and Eq. (28)). When developing preconditioners for the linear elasticity equation, typically a special care has to be taken to reproduce the action of $\mathbf{A}$ on the rigid body modes. In our case, the preservation of these modes is in fact a by-product of the design of our methods.

7.3 Factorization times and memory usage

The theoretical scalings predicted by Proposition 5.1 are confirmed by our tests. In Fig. 7.3 we show examples of hierarchical factorization timings. In Fig. 7.4 we show example memory usage. All plots are well within $\mathcal{O}(n)$ bounding lines.

7.4 Comparison to standard low-rank approximation

For each tested scheme (Nest-All-All, Nest-2-All, Nest-2-2, Gen-All-All), we compare our results to the ones obtained when polynomial compression is replaced by a standard low-rank approximation (used by other hierarchical solvers, e.g., [45, 11, 14, 51, 20, 31]). Namely, when computing the factorization with the given elimination and compression strategy to obtain $\mathbf{A}_{\ell},$ we record the sizes of all $\mathbf{Q}_{1}$ matrices from Eq. (7) used in compression. We then compute the factorization again with exactly the same set-up, except this time, to obtain the matrix $\mathbf{Q}$ from Eq. (7), we find a low-rank approximation to the off-diagonal blocks using the column-pivoted rank-revealing QR [12, 23]:

[TABLE]

However, the split $\mathbf{Q}=\begin{pmatrix}\mathbf{Q}_{1}&\mathbf{Q}_{2}\end{pmatrix}$ is such that the number of columns of $\mathbf{Q}_{1}$ is the same as before. As a result, the sizes of all nonzero matrix blocks when computing the factorization are the same, except inside the compression step. Denoting the obtained operator by $\mathbf{\hat{A}}_{\ell}$ we have that applying $\mathbf{A}_{\ell}^{-1}$ or $\mathbf{\hat{A}}_{\ell}^{-1}$ involves exactly the same cost. We call the obtained preconditioner the standard low-rank equivalent which is algebraically very close to the recently introduced methods [14, 11]. We note here that more accurate rank-revealing factorizations, such as SVD, could be used but are often impractical because of their computationally expensive iterative nature.

7.4.1 Improved accuracy on smooth near-kernel eigenvectors

Since in the case of the constant-coefficient Poisson equation eigenvectors and eigenvalues are known exactly, we can compute the relative backward error $\frac{1}{\lambda_{i}}\|(\mathbf{A}-\mathbf{A}_{\ell})\mathbf{v}_{i}\|_{2}$ for the unit-length eigenvectors $\mathbf{v}_{1}$ and $\mathbf{v}_{n}$ corresponding to, respectively, the largest, and the smallest, eigenvalues $\lambda_{1}$ and $\lambda_{n}$, to observe the effect of Corollary 6.3. Example results for Nest-All-All using polynomial compression, and its standard low-rank equivalent, are shown in Table 7.1. The error on $\mathbf{v}_{1}$ is small and similar in both cases; on the other hand, the error on $\mathbf{v}_{n}$ is nearly constant when using polynomial compression but for its standard low-rank equivalent, the error grows with increasing system sizes and becomes three orders of magnitude larger than when using polynomial compression.

7.4.2 Improved iteration counts

By design, applying our standard low-rank equivalent of the given preconditioner involves identical computational cost as applying the original preconditioner. Therefore, looking at the CG iteration counts is an exact way to compare the two approaches in terms of the quality of the resulting preconditioning operators. On all test cases, we record the CG iteration counts needed for convergence. The iteration counts are nearly constant or grow very slowly with increasing problem sizes when using our approaches, which is often not true for their standard low-rank equivalents. The plots of CG iteration counts are shown in Fig. 7.5, Fig. 7.6, and Fig. 7.7. Combined with the theoretical $\mathcal{O}(n)$ complexity of applying $\mathbf{A}_{\ell}^{-1},$ the results mean that we can expect our algorithms to achieve scalings close to the optimal linear scaling of the total solution time.

7.4.3 Nearly optimal scalings of solution times

In Fig. 7.8 and Fig. 7.9, we plot the CPU times needed to solve the two more challenging problems (Eq. (26), Eq. (27)). The graphs corresponding to our methods appear parallel to the $\mathcal{O}(n)$ bounding lines which is not true for many of the standard low-rank equivalents.

7.5 Choosing the appropriate preconditioner

Nest-2-2 exhibits lowest and nearly constant CG iteration counts in all tested cases. Inevitably, it has largest memory requirements as well as the cost of applying $\mathbf{A}_{\ell}^{-1}$. Nest-2-All and Nest-All-All performed almost identically. Therefore, Nest-All-All is likely preferable of the two, with $\mathcal{O}(n)$ factorization complexity guarantees (instead of $\mathcal{O}(n\log{n})$). Polynomial compression allows therefore for bounding the sizes of the nodes appearing in Nest-All-All, without losing accuracy (see the proof of Proposition 5.1). We expect this to be of importance for parallel scaling of the algorithm.

In the tests, Gen-All-All performed competitively in terms of solution times despite higher iteration counts. Notice that $\mathbf{A}_{\ell}^{-1}$ is a product of block diagonal matrices with one block per each node. This suggests that Gen-All-All also has very promising properties for massive parallelization. However, higher iteration counts may make Gen-All-All less practical on very ill-conditioned problems.

8 Conclusions and future work

The numerical results show that Sparse Geometric Factorization (SGF) gives rise to robust preconditioners, exhibiting optimal or near optimal scalings of solution times. Our methods are based on a special treatment of the near-kernel smooth eigenvectors, directly targeting the fundamental limitations of hierarchical matrix approaches. Moreover, the factorization is guaranteed to succeed in exact arithmetic, and the preconditioners are SPD operators, so Conjugate Gradient can always be used. Also, polynomial compression can allow avoiding costly rank-revealing decompositions—which have been reported as the main computational bottleneck of hierarchical solvers [11, 45, 51]. We did not find polynomial compression to be a bottleneck in our implementation.

Other strategies for designing preconditioners based on SGF than described in this paper, can be considered. While our examples were discretized on regular cartesian grids, partitionings on general grids are applicable (see e.g. [11] for details). Also, combining the polynomial compression with low-rank approximation may result in more efficient methods. In the absence of geometrical information (e.g., when only the system matrix is known), piecewise constant vectors would be then used. Moreover, polynomial compression can be applied almost verbatim in an approach that would only compress the well-separated interactions in a general domain partitioning (such approaches include [45, 51]).

Similar to [31, 11] our algorithms have very promising parallel properties. Thanks to compression, most computations are performed in the initial levels of the algorithm, contrasting direct solvers such as nested dissection multifrontal, where the exact factorization of the largest block can dominate the computations. We believe that the bounded sizes of nodes ensured by polynomial compression may be of further importance for parallel scaling of the algorithm. Efficient parallel implementations of hierarchical approaches such as ours is a topic of active research [16, 36].

Acknowledgments

This work was partly funded by the Stanford University Petroleum Research Institute’s Reservoir Simulation Industrial Affiliates Program (SUPRI-B). We are very grateful to prof. Hamdi Tchelepi for interest in our methods and useful suggestions. We would also like to thank Pavel Tomin and Sergey Klevtsov for helping us prepare the SPE10 reservoir benchmark case, and Léopold Cambier for invaluable discussions about the hierarchical factorizations and test cases.

The computing for this project was performed on the Stanford University Sherlock cluster. We would like to thank the Stanford Research Computing Center for providing computational resources and support.