FAST HIERARCHICAL SOLVERS FOR SPARSE MATRICES USING
EXTENDED SPARSIFICATION AND LOW-RANK
APPROXIMATION∗

arXiv:1510.07363v3 [math.NA] 15 Dec 2016

HADI POURANSARI† , PIETER COULIER† ‡ , AND ERIC DARVE† §
Abstract. Inversion of sparse matrices with standard direct solve schemes is robust, but computationally expensive. Iterative solvers, on the other hand, demonstrate better scalability, but
need to be used with an appropriate preconditioner (e.g., ILU, AMG, Gauss-Seidel, etc.) for proper
convergence. The choice of an effective preconditioner is highly problem dependent. We propose a
novel fully algebraic sparse matrix solve algorithm, which has linear complexity with the problem
size. Our scheme is based on the Gauss elimination. For a given matrix, we approximate the LU
factorization with a tunable accuracy determined a priori. This method can be used as a stand-alone
direct solver with linear complexity and tunable accuracy, or it can be used as a black-box preconditioner in conjunction with iterative methods such as GMRES. The proposed solver is based on the
low-rank approximation of fill-ins generated during the elimination. Similar to H-matrices, fill-ins
corresponding to blocks that are well-separated in the adjacency graph are represented via a hierarchical structure. The linear complexity of the algorithm is guaranteed if the blocks corresponding to
well-separated clusters of variables are numerically low-rank.
Key words. sparse, hierarchical, low-rank, elimination, compression
AMS subject classifications. 65F05, 65F08, 65F50, 65N55, 68Q25

1. Introduction. In the realm of scientific computing, solving a sparse linear
system,
(1)

Ax = b,

is known to be one of the challenging parts of many calculations, and is often the
main bottleneck. Such a system of equations may be the result of the discretization
of some partial differential equation (PDE), or more generally, can represent the local
interactions of units in a network.
Solving a system of equations of size n using a naive implementation of Gauss
elimination has O(n3 ) time complexity. The best proved time complexity to solve a
general linear system is O(nω ), where ω < 2.376 [12, 17, 40]. In the case of sparse
matrices, the time and memory complexity can be reduced when a proper elimination
order is employed. Finding the optimal ordering (that results in the minimum number
of new non-zeros in the LU factorization) is known to be an NP-complete problem
[57]. For matrices resulting from the discretization of some PDE in physical space,
nested dissection [23, 43] is known as an efficient elimination strategy. [1] discusses
the complexity of nested dissection based on the sparsity pattern of the matrix. For
a three-dimensional problem, the time and memory complexities are expected to be
O(n2 ) and O(n4/3 ), respectively, when nested dissection is employed. As the size of
the problem grows, such complexities make direct solvers prohibitive.
∗ Funding from the “Army High Performance Computing Research Center” (AHPCRC), sponsored
by the U.S. Army Research Laboratory under contract No. W911NF-07-2-0027, at Stanford, supported in part this research. The second author is a post-doctoral fellow of the Research Foundation
Flanders (FWO) and a Francqui Foundation fellow of the Belgian American Educational Foundation
(BAEF). The financial support is gratefully acknowledged.
† Stanford University, Department of Mechanical Engineering, Stanford, CA 94305, USA
([hadip,pcoulier,darve]@stanford.edu).
‡ KU Leuven, Department of Civil Engineering, Kasteelpark Arenberg 40, 3001 Leuven, Belgium.
§ Stanford University, Institute for Computational and Mathematical Engineering, Stanford, CA
94305, USA.

1

2

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE

Iterative methods, such as conjugate gradient [36], minimum residual [47], and
general minimum residual [53], are generally more time and memory efficient. In addition, iterative solvers such as those based on Krylov subspace can be accelerated
using fast linear algebra techniques. The fast multipole method (FMM) [20, 22, 28,
45, 49, 58], for example, can accelerate matrix-vector multiplication —from quadratic
complexity to linear—, which is the bulk calculation in iterative solvers based on
Krylov subspace. However, in practice, iterative methods need to be used in conjunction with preconditioners to limit the number of iterations. The choice of an efficient
preconditioner is highly problem dependent. There are many ongoing efforts to develop preconditioners that are optimized for particular applications. Hence, there is a
need for general purpose preconditioners. Hierarchical matrices enable us to develop
such preconditioners.
FMM matrices are a subclass of a larger category of matrices called hierarchical matrices (H-matrices) [7, 9, 32]. H-matrices have a hierarchical low-rank structure. For instance, in a hierarchically off-diagonal low-rank (HODLR) matrix [5],
off-diagonal blocks can be represented through a hierarchy of low-rank interactions.
If the bases used in the hierarchy are nested (i.e., the low-rank basis at each level is
constructed using the low-rank basis of the child level) the method is called hierarchically semi-separable (HSS) [3, 15, 56]. In a more general case of hierarchical matrices,
more complex low-rank structures can be considered. A full-rank dense matrix with
many low-rank structures is in fact data-sparse [31, 33]. A data-sparse matrix can be
represented via an extended sparse matrix, which has extra O(n) rows/columns, but
with only few non-zero entries [2]. The hierarchical structure of such matrices can be
used for efficient calculation and storage.
Recently, hierarchical interpolative factorization [38, 39] was proposed, which can
be used to directly solve systems obtained from differential and integral equations
based on elliptic operators. The fast factorization is obtained by skeletonization
of fronts in the multifrontal scheme. Using low-rank structure of the off-diagonal
blocks to develop fast direct solvers for linear systems arising from integral equations
has been widely studied [18, 26, 27, 41]. [25] proposed a direct solver for elliptic
PDEs with variable coefficients on two-dimensional domains by exploiting internal
low-rank structures in the matrices. [44] used hierarchical low-rank structures of the
off-diagonal blocks to introduce a preconditioner for sparse matrices based on a multifrontal variant of sparse LU factorization. [46] introduced a black-box linear solver
using tensor-train format. [42] used a recursive low-rank approximation algorithm
based on the Sherman-Morrison formula to obtain a preconditioner for symmetric
sparse matrices.
Sparse matrices can be considered as a very special case of hierarchical matrices,
where instead of low-rank blocks they initially have zero blocks. However, during
the elimination process in a direct solve scheme, many of the zero blocks get filled.
For a large category of matrices, including those obtained from the discretization of
PDEs, most of the new fill-ins are numerically low-rank. This is justified when the
Green’s function associated to the PDE is smooth (non-oscillatory). In this paper,
we will use the H-matrix structure to compress the fill-ins. A similar process can be
applied in the elimination of an extended sparse matrix resulting from an originally
dense matrix [4, 19]. This reduces the complexity of the direct solver to linear. The
linear complexity of the method is guaranteed if the blocks corresponding to the
interaction of well-separated nodes are numerically low-rank. We define the wellseparated condition in section 3.
The proposed algorithm can be considered as an extension to the block incomplete

FAST HIERARCHICAL SPARSE SOLVERS

3

LU (ILU) [52] preconditioners. In a block ILU factorization, most of the new fill-ins
(i.e., blocks that are created during the elimination process which are originally zero)
are ignored, and therefore, the block sparsity of the matrix is preserved, while the
accuracy is not. In the proposed algorithm, instead, we use low-rank approximations
to compress new fill-ins. Using a tree structure, new fill-ins at the fine level are
compressed and pushed to the parent (coarse) level. The elimination and compression
processes are done in a bottom-to-top traversal.
In addition, the proposed algorithm has formal similarities with algebraic multigrid (AMG) methods [10, 11, 51, 54]. However, the two methods differ in the way they
build the coarse system, and use restriction and prolongation operators. In AMG, the
original system is solved at different levels (from fine to coarse). Here, the compressed
fill-ins —corresponding to the Schur complements— of each level are solved at the
coarser level above. Note that the proposed algorithm is purely algebraic, similar
to AMG. If the matrix comes from discretization of a PDE on a physical grid, the
grid information can be exploited to improve the performance of the solver, similar
to geometric multi-grid.
The algorithm presented in this paper computes a hierarchical representation of
the LU factorization of a sparse matrix using low-rank approximations. We introduce intermediate operations to compress new fill-ins. The compressed fill-ins are
represented using a set of extra variables. This technique is known as extended sparsification [14]. The accuracy of the factorization phase (i.e., Gauss elimination and
compression), , can be determined
a priori. The time and memory complexity of the

factorization are O n log2 1/ and O (n log 1/), respectively, as will be clarified in
subsection 6.1.
The method presented in this paper is similar to the fast hierarchical methods
developed by Hackbusch et al. [9, 31, 32, 33] in a sense that both methods use a
tree decomposition to identify and represent low-rank blocks. The key difference,
however, is that in the Hackbusch’s algorithm the LU factorization is computed using
a depth-first tree traversal order, whereas here we use a breadth-first (level by level
from leaf to root) traversal. The connection and differences of the proposed algorithm
and Hackbusch’s fast H-algebra is discussed in our companion paper [19], in which a
similar method is used for dense matrix factorization.
Our solver can be used as a stand-alone direct solver with tunable accuracy. The
factorization part is completely separate from the solve part and is generally more
expensive. This makes the algorithm appealing when multiple right hand sides are
available (e.g., using the proposed solver as a black-box preconditioner in an iterative
method). We have implemented the algorithm in C++ (the code can be downloaded
from bitbucket.org/hadip/lorasp), and benchmarked it as both a stand-alone solver
(see subsection 6.1), and a preconditioner in conjunction with the generalized minimum residual (GMRES) iterative solver [53] (see subsection 6.2).
Furthermore, the proposed algorithm has interesting parallelization properties.
On one hand, all calculations are block matrix computations which can be highly accelerated using BLAS3 operations [6]. On the other hand, since the sparsity pattern
at every level is preserved, the data dependency is very local, which is an interesting property to reduce the amount of communications. In addition, the amount of
calculation scales with the third power of the size of blocks, while the communications scales with the second power of block sizes. This helps with the concurrency of
the parallel implementation. Moreover, the order of elimination does not change the
complexity of the presented algorithm. This is in particular an appealing property
for parallel implementation. The parallel implementation of the proposed method is

4

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE

not further discussed in this paper.
The remainder of this paper is organized as follows. In section 2 we briefly
introduce a graph representation of sparse matrices, and an interpretation of the
Gauss elimination using the adjacency graph. In section 3 some concepts related to
the hierarchical representation of matrices are defined. The algorithm is explained in
section 4 in detail, and the linear complexity analysis is provided in section 5. We
present numerical results obtained from various benchmarks in section 6. There are
many avenues for optimization and extension of the algorithm. We discuss some of
these opportunities in section 7.
2. Sparse linear systems. In this section we briefly introduce the graphical
framework that is required in the rest of the paper. We assume a sparse linear system
of size n as in (1) is given.
2.1. Adjacency graph. In many algorithms, including the method proposed
in this paper, it is necessary (or more efficient) to operate on sub-blocks of the matrix rather than single elements. The blocks of the matrix can be identified using a
partitioning as defined below.
Definition 1. (partitioning) A partitioning P is defined as a surjective map
{1, . . . , n} → {1, . . . , nP }. P groups rows/columns of A into nP clusters, CjP :=
{k ∈ {1, . . . , n}|P(k) = j} for 1 ≤ j ≤ nP .
We denote an entry of matrix A located in row k and column t by A[k,t] . For
1 ≤ i, j ≤ nP , we use Ai,j to represent a sub-matrix formed by concatenating all
entries A[k,t] such that k ∈ CiP and t ∈ CjP . Additionally, for a vector x of size n,
we use xi to represent a sub-vector formed by concatenating all x[k] entries such that
k ∈ CiP .
It is often fruitful to represent sparse matrices using graphs. An adjacency graph,
as defined below, represents a sparse matrix with partitioning.
Definition 2. (adjacency graph) A sparse matrix A with a partitioning P can
be represented by its adjacency graph G(V, E), where V = {v1 , . . . , vnP }. Each vi ∈ V
for 1 ≤ i ≤ nP represents a cluster CiP of rows and columns of A. A vertex vi is
connected to a vertex vj by a directed edge evi →vj ∈ E if and only if the block Aj,i in
A is non-zero.1
In Figure 1, an example of the adjacency graph is illustrated. In the rest of the
paper we use vertex and node interchangeably for the elements of V in the adjacency
graph.
The linear system in (1) can also be represented using the adjacency graph of
A, G(V, E). For a node vi ∈ V , Var(vi ) = xi denotes the vector of variables corresponding to cluster CiP . Similarly, RHS(vi ) = bi denotes the vector of right hand
sides corresponding to cluster CiP . Also for an edge evi →vj ∈ E, Mat(evi →vj ) = Aj,i
denotes the sub-matrix corresponding to cluster CiP of columns and cluster CjP of
rows. For the example shown in Figure 1, the following two notations represent the
same set of equations corresponding to the node v2 and its incoming edges.
(2a)

A2,1 x1 + A2,2 x2 + A2,3 x3 + A2,4 x4 = b2

1 The adjacency graph of a matrix A with a partitioning P is essentially the quotient graph of the
adjacency graph of matrix A with identity partitioning, where the equivalence relation is induced by
partitioning P.

5

FAST HIERARCHICAL SPARSE SOLVERS

(2b)

Mat(ev1 →v2 ) · Var(v1 ) + Mat(ev2 →v2 ) · Var(v2 )+

Mat(ev3 →v2 ) · Var(v3 ) + Mat(ev4 →v2 ) · Var(v4 ) = RHS(v2 )
c1

c2

c3

c4

c1

x1

c2

x2

c3

x3

c4

x4

v3

b1

=

v1

b2
v2

b3
b4

v4

Fig. 1. Example of the adjacency graph (right) of a matrix (left). Vertices’ colors are the same
as their corresponding cluster of rows/columns in the matrix. Edges’ colors are also in correspondence with the sub-blocks in the matrix.

2.2. Elimination. The block Gauss elimination process, or block LU factorization, can also be explained using the graph representation of the matrix. At step i
of the elimination process, a set of unknowns, Var(vi ), is eliminated from the system
of equations. This corresponds to eliminating vertex vi from the adjacency graph
G(V, E). The self-edge from vi to itself corresponds to the pivot diagonal sub-block
in the matrix. After eliminating vi , for every pair of outgoing edge evi →vj to a vertex
vj and incoming edge evk →vi from a vertex vk , a new edge from vk to vj is created,
corresponding to the Schur complement of the eliminated edges, that is
(3)

−Mat(evi →vj ) · Mat(evi →vi )−1 · Mat(evk →vi ) = −Aj,i A−1
i,i Ai,k

Note that if the edge between vk and vj exists before elimination, the Schur complement adds to the existing sub-block.
The process described above reveals the fact that during the elimination process
many new edges are introduced in the graph. This corresponds to generating new
non-zero blocks in the matrix during the LU factorization. The generation of many
dense blocks is what makes the direct factorization of sparse matrices a prohibitive
process. Essentially, a matrix A can be sparse, while L and U in the LU factorization
of A are dense. In the next section, we explain how we can preserve the sparsity of the
matrix during the elimination process by compressing the well-separated interactions.
This process is known as extended sparsification.
2.3. Key idea. An important observation in the elimination process is the fact
that fill-ins (i.e., new edges created during the elimination process) that correspond to
well-separated vertices are often numerically low-rank. For a linear system obtained
from a discretized PDE, well-separated vertices refers to points that are physically far
enough from each other. For a general sparse matrix two vertices are well-separated
if their distance in the adjacency graph is large enough. It is formally defined in
Definition 6. We replace such fill-ins with a sequence of low-rank matrices.
For example, consider the following symmetric linear system that is partitioned
into 3 blocks

   
S B C
x1
b1
 x2  = b2 
B | P
(4)
x3
b3
C|
Q

6

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE

-I
x2

P

B
T

B

x1

S
(a)

C
C

T

x3

x2

Q

P'

DT
D
(b)

U

T z2

x3

x2 U

Q'

P'

y2

KT

-I K

y3

-I

-I
z3
T

V

V
x3

Q'
(c)

Fig. 2. An example of compression of a well-separated interaction. (a) The adjacency graph of
the original linear system described in (4). (b) The resulting graph after eliminating x1 node. (c)
The adjacency graph of the extended system. All edges are labeled with their corresponding block in
the matrix.

In Figure 2a the adjacency graph of the system of equations in (4) is shown. Now,
consider eliminating x1 . The resulting system is as follows
 0
   0 
P
D
x2
b2
(5)
=
,
D| Q0
x3
b03
where D = −B | S −1 C, P 0 = P −B | S −1 B, Q0 = Q−C | S −1 C, b02 = b2 −B | S −1 b1 , and
b03 = b3 − C | S −1 b1 . The adjacency graph of the system of equations in (5) is depicted
in Figure 2b. Nodes 2 and 3 can be considered well-separated (see Definition 6). They
get connected due to the elimination of node 1. We assume their interaction is lowrank, and can be written as
(6)

D ' U KV | ,

where U and V are tall matrices. We can use any low-rank approximation in (6), for
example the singular value decomposition (SVD). We combine (5) and (6) to define
a new set of equations, in which direct interaction of the nodes 2 and 3 is replaced by
a sequence of low-rank interactions
 0
   0 
P
U
x2
b2
0

 x3  b03 
Q
V
 |
   
U
 z 2   0 
−I

  =  
(7)
|

   
−I 
V

 z 3   0 


−I
K y 2   0 
|
−I K
y3
0
In the above equation, we introduced extra variables (y 2 , y 3 , z 2 , and z 3 ) to represent the far-field interactions. This technique is known as extended sparsification
[14]. Note that the system of equations (7) is equivalent to the system of equations
(5) up to the accuracy of (6). In fact, if we do a Gaussian elimination on the matrix
in (7) using 2 × 2 blocks starting from the right, we get back to the original matrix
in (5).
In an analogy with the fast multipole method, y 2 and y 3 can be considered as
the multipole coefficients, and z 2 and z 3 as the local coefficients. In Figure 2c the
adjacency graph of the extended linear system is depicted. The red and black nodes

7

FAST HIERARCHICAL SPARSE SOLVERS

correspond to the extended part of the matrix in (7) (extra variables introduced
above). As a result of extended sparsification, nodes 2 and 3 are disconnected in
Figure 2c. Therefore, in the general case x2 could be eliminated without introducing
new fill-ins between x3 and other neighbors of x2 .2 As shown In Theorem 9, removing
well-separated interactions of nodes before elimination preserves the sparsity.
In the general case, we start by a standard block Gauss elimination. When we
create fill-ins, similar to the above example, we apply the extended sparsification
(compression). This results in the creation of new nodes in the graph belonging to a
coarser level.
As we proceed with the elimination process at level i, edges between the auxiliary
nodes at level i − 1 are created. After the elimination process at level i is completed,
we proceed with the elimination at level i − 1, and continue up to the root (level 0).
Essentially, in this process we form a hierarchical approximation of L and U matrices
(which are generally dense) in the LU factorization of A, without directly computing
them.
In addition, similar to the agglomeration process in multi-grid methods, we consider one red-node and one black-node for every pair of clusters (i.e., the well-separated
interactions of each pair of clusters are compressed together). Therefore, the number
of red-nodes at the parent level is half of the number of red-nodes at the current level.
We formally define the hierarchy of nodes in section 3, and the details of the algorithm
in section 4.
3. Hierarchical representation. To form a hierarchical tree, we recursively
partition the rows/columns indices of the matrix. This is formally defined as a sequence of nested partitionings.
Definition 3. (nested partitionings) A sequence of partitionings P0 , . . . , Pl with
P
nPi = 2i is called nested if every cluster CjPi is the union of two clusters Cj 0 i+1 and
P

Cj 00i+1 for 0 ≤ i < l. Cluster CjPi is then called the parent of two child clusters
P

P

P

P

Cj 0 i+1 and Cj 00i+1 . We call Cj 0 i+1 and Cj 00i+1 sibling clusters. The level of a cluster
is defined as the index i of its defining partitioning Pi . The clusters associated to the
finest partitioning, Pl , have no children and are called leaf clusters, while the cluster
associated to P0 has no parent and is called the root cluster. All other clusters have
exactly two children and one parent.3
A visual example of nested partitionings with six levels (i.e., P0 , . . . , P5 ) is illustrated
in Appendix B. Now, we define the hierarchical tree4 (denoted by H-tree) of a sparse
matrix A given a sequence of nested partitionings.
Definition 4. (hierarchical tree) Given a sparse matrix A, and a nested sequence
of partitionings P0 , . . . , Pl , the hierarchical tree (H-tree) is defined as a directed graph
with red and black vertices and two types of edges (parent-child and interaction edges).
The vertices of H-tree are corresponding to the clusters associated with the nested
partitionings. For every cluster CjPi with 0 ≤ i < l and 1 ≤ j ≤ 2i there is a
[i+1]

corresponding black-node bj
red-nodes

[i+1]
r0 j

and

[i+1]
r1 j

in the vertex set of H-tree. Also, there is a pair of
P

P

corresponding to the sibling clusters Cj 0 i+1 and Cj 00i+1

2 Note that in this example, we intentionally chose order of elimination x , x , x , to introduce
1
2
3
a new fill-in, and demonstrate the low-rank compression process. Using x2 , x1 , x3 , as the order of
elimination would result in no fill-ins for this example. Finding such an order of elimination, however,
is a hard problem and may not be possible for practical matrices.
3 The binary subdivision is not necessary. We can generalize this to quad-tree, octree, etc.
4 In fact, it is not a tree. As we see later, there are edges between nodes at each level.

8

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE

(children of the cluster CjPi ) in the vertex set of H-tree. We call such a pair of red[i+1]

nodes a super-node, and denote it by sj
[i+1]

bj

[i+1]

is connected to r0 j

the parent of super-node

[i+1]

and r1 j

[i+1]
sj .

, that is corresponding to the cluster CjPi .
[i+1]

by parent-child edges. Hence, bj

is also

Additionally, the red-node corresponding to CjPi is

[i+1]
bj

connected to
by a parent-child edge. We denote the parent of a node v (which
can be a red, black, or super node) by P(v).
We also consider one special red-node associated to the root cluster. The level
of each vertex of H-tree is denoted by a superscript index. Additionally, the depth of
H-tree is defined as l.
There is an interaction edge between two red-nodes with level l (leaf red-nodes),
if the corresponding vertices in the adjacency graph of A with partitioning Pl are
[l]
[l]
connected. Therefore, the subgraph of H-tree induced by r0 j and r1 j for 1 ≤ j ≤ 2l−1
is the adjacency graph of A with partitioning Pl .
There is no interaction edge between non-leaf vertices (shown transparent in Figure 3) of an H-tree before applying the elimination. Non-leaf nodes are reserved to be
used for the extended sparsification similar to the red and black nodes in Figure 2c.
Similar to the vertices of an adjacency graph, each node of an H-tree also corresponds to a set of variables and equations. Leaf-nodes of the H-tree correspond
to the variables and equations in (1) partitioned using Pl . Non-leaf nodes of the Htree, however, correspond to the auxiliary variables and equations. In subsection 4.3,
when we explain the extended sparsification in the general case, we will introduce the
variables and equations corresponding to the non-leaf nodes of the H-tree.
An example of an H-tree is depicted in Figure 3. A parent-child edge between
two nodes is shown by a dashed line, whereas interaction edges are shown by solid
lines. In the rest of the paper, we use edge and interaction edge interchangeably, while
parent-child edges are explicitly mentioned.

Fig. 3. An example of a hierarchical tree. Dashed lines show parent-child edges, and solid lines
represent interaction edges. Non-leaf nodes (shown transparent) have no interaction initially, and
are reserved to represent the well-separated interactions at the level below them.

Definition 5. (adjacent clusters) Consider a sparse matrix A and a nested sequence of partitionings P0 , P1 , P2 , . . . , Pl . Two leaf clusters CpPl and CqPl are called
adjacent (or neighbors) iff Ap,q or Aq,p are non-zero blocks (i.e., nodes vp and vq in
the adjacency graph of A with partitioning Pl are connected). Two clusters C and C 0

FAST HIERARCHICAL SPARSE SOLVERS

9

(not necessarily with the same level) are adjacent iff a leaf descendant5 cluster of C
is adjacent to a leaf descendant cluster of C 0 .
Definition 6. (well-separated nodes in the H-tree) Nodes u and v in the vertex
set of an H-tree (which can be red, black, or super-nodes at different levels) are wellseparated if their corresponding clusters are not adjacent. An interaction edge that
connects two well-separated nodes is called a well-separated edge. Well-separated edges
do not exist initially in the H-tree, and are created as a result of elimination.
Low-rank interaction is typically observed in physical systems when two clusters
are sufficiently far apart from each other (e.g., well-separated clusters in fast multipole method). We can replace this condition by an equivalent distance requirement
in the graph. However, for simplicity we weakened the requirement and declare that
two clusters are well-separated only if they are not adjacent (Definition 6). If the
partitioning of the domain is adequate, this simple definition is sufficient to approximate the more accurate separation requirement based on distance. In many cases,
well-separated interactions are low-rank, but there is no guarantee. It depends on
various details not studied in this work, such as the shape of the clusters.
4. Algorithm. In the previous section we defined the H-tree, which is used
to represent graphically the extended system, similar to the example provided in
Figure 2c. In this section, we explain the details of the algorithm to compute a
hierarchical representation of an approximate LU factorization of matrix A. Note that
while A is sparse, the L and U matrices are typically dense. However, we assume that
they can be represented using a hierarchical tree (through the extended sparsification).
We could find the matrices L and U through an elimination process, and compute
their hierarchical representation using an extended sparsification. However, in order
to achieve linear complexity, we perform elimination and the extended sparsification
(compression) together.
The algorithm presented in this paper takes advantage of similar technique as in
the inverse fast multipole method (IFMM) [19]. The IFMM can be used to compute
the hierarchical representation of LU factorization of a dense matrix which is given
in hierarchical form (i.e., FMM matrix). Essentially, in the IFMM, we are given an
H-tree that has interaction edges at all levels. This H-tree represents a dense matrix.
In the sparse case, however, the H-tree initially has interaction edges only at
the leaf level. For both of the sparse and dense cases, we start with an H-tree that
represents matrix A, and end up with an H-tree representing an approximate LU
factorization of A. In Algorithm 1 the overall factorization scheme is introduced.
Various sub-algorithms are explained afterwards. In addition, a step by step example
of the elimination process on the H-tree and the corresponding extended matrix is
presented in Appendix A. Similar to the standard LU factorization, after the elimination process we are able to efficiently compute A−1 b through forward and backward
substitutions.
4.1. Initializing the H-tree. The Initialize(l) function in Algorithm 1 consists of computing l nested partitionings and form the H-tree with depth l, as defined
in Definition 4. An example of an H-tree and the corresponding matrix is depicted in
Figure 20. Leaf nodes and interaction edges of the H-tree initially represent the given
linear system of equations (1) with partitioning Pl . Through the elimination process,
we extend the system of equations, and use non-leaf nodes of the H-tree to represent
the new variables and equations.
5 Cluster

C1 is a descendant of cluster C2 iff C1 ⊆ C2 .

10

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE

Algorithm 1: Factorization using H-tree.
Input: sparse matrix A
Initialize(l)
// form an initial H-tree with l levels
for i ← l to 1 do
// iterate over levels from leaf to root
for j ← 1 to 2i−1 do
// iterate over nodes at each level
[i]
[i]
[i]
sj ← MergeRedNodes(r0 j , r1 j )
// form super-node
for j ← 1 to 2i−1 do
[i]
Compress(sj ) // compress well-separated interactions of super-node
[i]

Eliminate(sj )
[i]

Eliminate(bj )

// Gauss elimination for super-node
// block Gauss elimination for black-node

Output: LU factorization of A stored hierarchically in the H-tree
4.2. Forming super-nodes. The outer-loop in Algorithm 1 is over different
levels from the bottom to the top of the tree. At each level, we start by merging
red-siblings into super-nodes. In Algorithm 1 this process is denoted by the function MergeRedNodes(). This process results in a coarser representation of the linear
system. We substitute interactions (i.e., edges) between any two pairs of red-nodes
[i]
[i]
[i]
[i]
[i]
[i]
(r0 j , r1 j and r0 j 0 , r1 j 0 ) with an interaction between super-nodes (sj and sj 0 ) as
follows.
!
Mat(er0 [i] →r0 [i] ) Mat(er1 [i] →r0 [i] )
j
j
j0
j0
(8)
Mat(es[i] →s[i] ) =
Mat(er0 [i] →r1 [i] ) Mat(er1 [i] →r1 [i] )
j
j0
j

j0

j

j0

Similarly, the variable and right hand side vectors corresponding to a super-node
is formed by concatenating the variables and right hand sides of its constituting rednodes.


[i]
[i]
[i]
(9)
Var(sj ) = concatenate Var(r0 j ), Var(r1 j )

(10)



[i]
[i]
[i]
RHS(sj ) = concatenate RHS(r0 j ), RHS(r1 j )

The process of merging red-nodes to form the super-nodes is illustrated in Figure 4.
The merging process is also depicted in Figure 21 in the example provided in Appendix A.
4.3. Compressing well-separated edges. The next sub-algorithm to consider
is the compression. In Algorithm 1 this process is denoted by the function Compress().
During the compression, well-separated interactions of a super-node are pushed to the
parent level nodes which represent a set of auxiliary variables. This is essentially the
extended sparsification method which we discussed in the example in subsection 2.3
(i.e., replacing the well-separated edges in Figure 2b by the sequence of edges between
red and black nodes at the parent level as shown in Figure 2c).
[i]
[i]
Assume we are at level i, and about to apply Compress(sj ). Also, assume sj is of
size m (i.e., corresponds to m variables and m equations), and interacts with (i.e., has

11

FAST HIERARCHICAL SPARSE SOLVERS

Fig. 4. An example of the merge process corresponding to the MergeRedNodes() function in
Algorithm 1. The leaf red-nodes of the H-tree in Figure 3 (top) are merged to super-nodes (bottom).

an edge to) t well-separated nodes p1 , p2 , . . . , pt with sizes m1 , m2 , . . . , mt , respectively.
As it will be clear later, the pk nodes for k = 1, . . . , t can either be a red-node (at the
parent level) or a super-node (at the same level). Assume blocks A1 , A2 , . . . , At are
[i]
associated to the outgoing well-separated edges from sj to p1 , p2 , . . . , pt , respectively,
i.e., Ak = Mat(es[i] →p ), where Ak is an mk by m matrix. Similarly, B1| , B2| , . . . , Bt|
j

k

[i]

are associated to the incoming well-separated edges to sj , i.e., Bk| = Mat(ep

[i]
k →sj

),

where Bk is an mk by m matrix. This is depicted schematically in Figure 5 (left).
Similar to the example in subsection 2.3, we assume well-separated edges can be
approximated using a low-rank factorization. We compress all well-separated interactions of a super-node together. We then introduce auxiliary variables, and replace the
well-separated edges by new edges between the auxiliary variables (i.e., going from
the left configuration to the right configuration in Figure 5).

[i]

Fig. 5. Schematic of the compression process for a super-node sj corresponding to the
Compress() function in Algorithm 1. Well-separated interactions are replaced with low-rank interactions with the red parent.

In order to compress all t well-separated edges together, we form a temporary
matrix by vertically concatenating A1 , . . . , At as well as B1 , . . . , Bt . Use a low-rank
approximation method (e.g., SVD) and write:
   
A1
R1
 ..   .. 
 .   . 
   
 At   R t  |
  '  V ,
(11)
B1  Q1 
   
 .   . 
 ..   .. 
Bt
Qt
where Rk and Qk are mk by r matrices for k = 1, 2, . . . , t, and V is an m by r matrix.

12

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE

r is the rank in the above low-rank approximation. From (11) we can write
Ak ' Rk V | and Bk| = V Q|k

(12)

for k = 1, 2, . . . , t

[i]

Var(sj ) contributes to the equation corresponding to a node pk with a term like
[i]

Ak · Var(sj ) for k = 1, 2, . . . , t. We can rewrite this term as


[i]
[i]
[i]
(13)
Ak · Var(sj ) ' (Rk V | ) · Var(sj ) = Rk V | · Var(sj )
Similarly, the contribution of Var(p1 ), . . . , Var(pt ) in the equation corresponding to
[i]
the node sj can be written as
(14)

t
X
k=1

Bk| · Var(pk ) '

t
X
k=1

(V Q|k ) · Var(pk ) = V

t
X
k=1

Q|k · Var(pk )

Now, we apply the extended sparsification, and introduce new variables and equations
as follows
[i]

[i]

[i]

[i]

(15)

Var(P(bj )) := V | · Var(sj ) ⇒ V | · Var(sj ) − Var(P(bj )) = 0

(16)

Var(bj ) :=

[i]

t
X
k=1

Q|k · Var(pk ) ⇒

t
X
k=1

[i]

Q|k · Var(pk ) − Var(bj ) = 0
[i]

[i]

In the above, we defined two new vector of variables, Var(P(bj )) and Var(bj ), each
of size r (similar to the vectors y 2 and z 2 in the example of subsection 2.3). We assign
[i]
[i]
equations (15) and (16) to the black-node bj and the red-node P(bj ), respectively.
Therefore,
(17)

[i]

[i]

RHS(bj ) = RHS(P(bj )) = 0

Now, we apply the following edge updates:
[i]
• Remove edges from sj to pk and vice versa for k = 1, 2, . . . , t.
[i]

• Add edges from P(bj ) to pk with blocks Rk for k = 1, 2, . . . , t.
[i]

• Add edges from pk to P(bj ) with blocks Q|k for k = 1, 2, . . . , t.
[i]

[i]

• Add an edge from sj to bj and vice versa with blocks V | and V , respectively.
[i]

[i]

• Add an edge from bj to P(bj ) and vice versa with blocks −Ir (minus identity
matrix of size r).
Therefore, in the compression process all well-separated edges connected to the super[i]
[i]
node sj are substituted with edges to/from P(bj ), as shown in Figure 5. Note that
each step of the compression process described here is in some sense “half” of the
extended sparsification step as described in subsection 2.3. For example, in Figure 22
[3]
[3]
there is a well-separated interaction between s2 and s4 . After the compression
[2]
[3]
step, this interaction is substituted with a well-separated edge between r1 1 and s4
[3]
(Figure 23). Ultimately, when we are about to eliminate s4 , we further compress
[2]
[2]
this edge and make connection between r1 1 and r1 2 (Figure 26), which is now the
same as the example in Figure 2.

13

FAST HIERARCHICAL SPARSE SOLVERS
[i]

As a result of the compression process on a super-node sj , we defined 2r new
[i]

[i]

auxiliary variables, corresponding to two nodes bj and P(bj ), each of size r. Note
that if the matrix is symmetric, Ak = Bk for k = 1, . . . , t. Therefore, we would not
need to concatenate Bk ’s in (11), and only half of the above calculations are required.
4.4. Elimination. After compressing all well-separated edges, we apply the
[i]
standard elimination. As a result of the compression process, the super-node sj
is only connected to its original neighbor nodes. This is a key property of the algorithm that preserves the sparsity of the matrix, which results in a slightly larger
system of equations (a constant number times the original size of the matrix). The
elimination process for a node is explained in subsection 2.2. We first eliminate the
[i]
[i]
super-node sj , and then eliminate its black-parent bj . In the example provided
in Appendix A, Figures 22, 24, 25 and 27 illustrate the graph and matrix after the
elimination process.
4.5. Solve. After the factorization part is completed, we can solve for multiple
right hand sides. The solve process consists of two steps: a forward and a backward
traversal of all nodes. This is identical to the standard forward and backward substitutions in the LU factorization. In the forward traversal we visit all nodes in the
order they have been eliminated, and in the backward traversal we visit nodes in the
exact reverse order. The solve process is introduced in Algorithm 2.
Note that in the factorization part we introduced auxiliary variables and equations
(i.e., all variables and equations associated to non-leaf nodes). We denote this extended
system of equations by Ae xe = be , which is equivalent to the system (1) up to the
accuracy of (11) (i.e., eliminating the auxiliary variables from the extended system
recovers the original system of equations). In the solve part, we have to solve for all
variables (original and auxiliary variables, i.e., xe ) even though we are just interested
in the original variables, x. The number of extra variables is limited (see Theorem 10),
and is of the same order as the number of original variables.
Algorithm 2: Solve for a given right hand side using H-tree.
Input: An H-tree representing an approximate LU factorization of A and
a right hand side vector b
SetRHS()
// set leaf RHSs to b and non-leaf RHSs to 0
/* forward substitution

for i ← l to 1 do
for j ← 1 to 2i−1 do
[i]
SolveL(sj )
[i]
SolveL(bj )

for i ← 1 to l do
for j ← 2i−1 to 1 do
[i]
SolveU(bj )

Output: x ' A−1 b

// iterate over nodes at each level
// update RHS of the super-node
// update RHS of the black-node

/* backward substitution

[i]
SolveU(sj )
[i]
SplitVar(sj )

*/
// iterate over levels from leaf to root

*/
// iterate over levels from root to leaf
// iterate nodes with reverse order
// solve for variables of the black-node
// solve for variables of the super-node

// split solution between the constituting red-nodes

14

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE

The solve algorithm begins with SetRHS(), that is to set the right hand side of
all nodes in the H-tree. As explained in Definition 4, leaf red-nodes of the H-tree
correspond to the original equation (1); therefore, the right hand side of each leaf
red-node is a sub-vector of b determined by the leaf-partitioning of the H-tree, Pl .
Based on (17), the right hand side of every non-leaf red-node and black-node is 0.
The right hand sides of super-nodes are computed by (10).
After setting the right hand side vectors, we apply functions SolveL() and
SolveU() to all super-nodes and black-nodes. The function SolveL() (see Algorithm 3) is applied through a forward traversal, and updates the right hand side
vectors. The function SolveU() (see Algorithm 4) is applied through a backward
traversal, and solve for the variables of each node. After solving for variables of a
super-node, we split the solution between its two constituting red-nodes (denoted by
function SplitVar()) according to (9). When Algorithm 2 is completed, the solution
vector x is formed by concatenating variable vectors of all leaf red-nodes.
Algorithm 3: Forward traversal (update the right hand sides)
Function SolveL(node p)
f ← Mat(ep→p )−1 · RHS(p)
for ep→q ∈ OutGoingEdges(p) do
if Order(q)>Order(p) then
// if q is eliminated after p
RHS(q) ← RHS(q) − Mat(ep→q ) · f

Algorithm 4: Backward traversal (solve for variables)
Function SolveU(node p)
Var(p)← RHS(p)
for eq→p ∈ InComingEdges(p) do
if Order(q)>Order(p) then
// if q is eliminated after p
Var(p) ← Var(p) − Mat(eq→p ) · Var(q)
Var(p) ← Mat(ep→p )−1 · Var(p)
In Algorithms 3 and 4 OutGoingEdges(p) and InComingEdges(p), respectively,
denote the set of all outgoing and incoming edges of a node p in the H-tree. The
function Order(p) returns the order of elimination of a node p, i.e., if in Algorithm 1
a node q is eliminated after a node p, then Order(q)>Order(p).
5. Linear complexity. In this section we show that the block sparsity of the
extended matrix is preserved through the elimination process. Therefore, the factorization algorithm has provable linear complexity provided that the block sizes (and
thus the rank of the low-rank approximations) are bounded.
In Definition 6 we defined well-separated nodes in an H-tree. In this section, we
generalize this concept, and define distance of nodes in a given H-tree.
Definition 7. (distance of nodes in H-tree) Consider nodes u and v in the vertex
set of an H-tree of a sparse matrix A with sequence of nested partitionings P0 , . . . , Pl .
P
Using Definition 4, u and v correspond to clusters CjPi and Cj 0 i0 for some 0 ≤ i, i0 ≤ l,

FAST HIERARCHICAL SPARSE SOLVERS

15

P

0

1 ≤ j ≤ 2i , and 1 ≤ j 0 ≤ 2i . Assume i ≤ i0 , and Cj 0 i0 is a descendant of cluster
CkPi . The distance between nodes u and v is defined as the length of (i.e., number of
edges) the minimum path between nodes vj and vk in the adjacency graph of A with
partitioning Pi .
Corollary 8. Nodes u and v in the vertex set of an H-tree are well-separated
iff their distance is greater than 1.
Note that other criteria are possible to define well-separated nodes (see section 7).
Theorem 9. (preservation of sparsity) In Algorithm 1, we never create an edge
between two nodes with distance greater than 2.
Proof. Elimination can results in connecting nodes at large distances. However, in
Algorithm 1, before applying elimination we remove all edges to nodes with distance
larger than 1 (i.e., the well-separated edges). Therefore, after eliminating a node we
create edges between nodes with distance at most 2.
Theorem 9 shows that for each node we need to process at most κ1 + κ2 edges,
where κ1 and κ2 are, respectively, the maximum number of super-nodes at distance 1
and 2 from a super-node. Note that κ1 and κ2 depend on the original matrix sparsity
pattern, and are independent of the size of the matrix. To establish linear complexity
of the factorization, we need to bound the size of the nodes (i.e., number of variables
associated to them) in the H-tree.
For matrices arising from the discretization of a PDE, well separated edges correspond to the interaction of points that are physically far from each other. Therefore, if
the Green’s function of the associated PDE is smooth enough, one can expect a wellseparated interaction to be numerically low-rank. We provide numerical evidence in
section 6 to support this argument. Note that the low-rank property of well-separated
nodes depends on the quality of partitioning and the definition of well-separation.
These are topics for followup studies.
For general sparse matrices we can guarantee the linear complexity through
bounding the rank growth. This is explained in the next theorem.
Theorem 10. (linear complexity condition) Consider di to be the maximum size
of super-nodes at level i of an H-tree with depth l resulted from Algorithm 1 such that

(18)

di ≤ αl−i dl

and

dl = O(n/2l ) = O(1),

√
where 0 < α < 3 2 is a constant number. Also assume κ1 and κ2 , the maximum
number of super-nodes at distance 1 and 2 from a super-node, are O(1) quantities.
Under these conditions the cost of the algorithm is linear with respect to the problem
size.

Proof. For a given super-node
at level i the compression cost is O κ2 d3i , and

the elimination cost is O κ21 d3i . Note that the required memory scales with d2i for
each node. Ignoring the constant factors κ1 and κ2 , the order of the total cost of
factorization is as follows:

(19)

factorization cost = O

l
X
i=1

!
2i−1 d3i

16

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE

Plug (18) in (19):
factorization cost = O
(20)

=O

l
X

!
2i−1 α3(l−i) d3l

i=1

2l−1 d3l

l−1  3 i
X
α
i0 =0

0

!

2

change of variable: i0 = l − i

= O(2l d3l )
Note that for α <

√
3

2 we have

Pl−1  α3 i0
i0 =0

2

= O(1). Furthermore, dl is the number

of variables in super-nodes at the leaf level which is O(n/2l ). Therefore:

(21)
factorization cost = O nd2l , factorization memory = O (ndl )

6. Numerical results. We have implemented the algorithm described in section 4 in C++. The code (we call it LoRaSp6 ) can be downloaded from bitbucket.
org/hadip/lorasp. We use Eigen [29] as the backend for linear algebra calculations,
and SCOTCH [48] for graph partitioning. We present results for various benchmarks,
where LoRaSp is used as a direct solver, or as a preconditioner in conjunction with
an iterative solver.
6.1. LoRaSp as a stand-alone solver. In this section we employ LoRaSp as
a stand-alone solver. The accuracy of the solver depends on the accuracy of the lowrank approximations during the compression step as explained in subsection 4.3. Any
low-rank approximation method can be used for the compression. Here, we use SVD.
For every well-separated interaction, we first compute the SVD, and then truncate the
singular values at some point. There are many possible criteria to truncate singular
values. We discuss some possible criteria. Figure 6 shows the decay of singular values
for blocks corresponding to the interaction between randomly chosen well-separated
nodes at different levels of an H-tree. The tree corresponds to a matrix obtained from
the second-order uniform discretization of the Poisson equation:
∇. (∇T ) = f

(22)

T =0

in D,

on ∂D

The domain D is a three-dimensional unit cube. The matrix size is 32,768, and the
depth of the corresponding H-tree is 11. Evidently, singular values have exponential
decay at different levels of the tree. The zero (up to machine precision) singular values
are not shown in the plot.
To demonstrate the linear complexity of the method, we considered a sequence
of problems with a growing number of variables. Consider the following sequence of
uniform discretization of the domain D in (22): 32 × 32 × 16, 32 × 32 × 32, 64 × 32 × 32,
64 × 64 × 32, 64 × 64 × 64, 128 × 64 × 64, and 128 × 128 × 64. The matrix size is
increased by a factor of 2 in the consecutive problems. Hence, to keep the size of
the leaf super-nodes constant among all problems, we consider H-trees with depth
10, 11, 12, 13, 14, 15, 16 for this sequence of problems, respectively. In general, the
depth of H-tree should scale linearly with log2 n, where n is the size of matrix.
6 Low

Rank Sparse solver.

17

FAST HIERARCHICAL SPARSE SOLVERS
0

10

Le ve l
Le ve l
Le ve l
Le ve l

−1

σ k /σ 0

10

11
10
9
8

−2

10

−3

10

−4

10

0

10

20
k

30

40

Fig. 6. Decay of singular values of random blocks corresponding to well-separated interactions
(to be compressed) at different levels of an H-tree with 11 levels resulted from Algorithm 1.

Well-separated edges corresponding to a block B, as shown in (11), with singularvalues σ0 , σ1 , ..., are compressed by keeping only the singular-values that satisfy:
σk
≥
(23)
σ0
Smaller values of  lead a to more accurate approximation of each block, and consequently a more accurate approximation of the final solution. For a given linear system
Ax = b, the precision of any solution x̃ is quantified by the relative error and relative
residual defined as follows
(24)

error =

kx − x̃k2
,
kxk2

residual =

kAx̃ − bk2
kbk2

Figure 7b shows that smaller values of  (i.e., more accurate low-rank approximations)
result in a more accurate estimation of the solution to the linear system in the cost of
larger factorization and solve times, as shown in Figure 7a. For a constant , the time
spent for the factorization and solve parts are asymptotically linear with the problem
size. Note that for smaller values of  the linear scaling is achieved for larger values
of n. In addition, the error and residual of the estimated solution for a fixed  barely
change with the problem size (see Figure 7b).
As it is clear from Figure 7, we can obtain more accurate solutions by decreasing
the parameter  in (23). To show the convergence of the solver, we picked a fixed
problem size, and measured the accuracy of the estimated solution as  decreases.
In addition, for comparison purposes, we consider a 2D variation of (22) which is
discretized using a finite volume approach with Voronoi tessellation. The points are
2
drawn from a random uniform distribution in the [0, 1] interval. The discretization
results in a matrix A = DB, where D is a diagonal matrix with inverse of the Voronoi
cells on the diagonal, and B is a symmetric matrix. We apply the factorization directly
to B. Note that the average number of non-zeros per row for a matrix corresponding
to a 2D Voronoi discretization is 7, which is the same as for a uniform second order
3D discretization.7
In Figure 8 the convergence of the solution for a 3D Poisson problem with n =
1.3 × 105 (corresponding to a 64 × 64 × 32 grid) and a 2D Poisson problem with the
7 We

can show this by double counting the angles in a 2D Voronoi tessellation, once through
points, and once through triangles of the corresponding Delaunay triangulations.

18

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE

10

4

10

0

ǫ = 0.01
10 3

time (sec)

10

2

ǫ = 0.04
ǫ = 0.16

10 -1

ǫ = 0.64
linear line

10 1

10 -2

10 0
10 -3
10

-1

10 -2 4
10

10 5

10 6

10 -4 4
10

10 5

n

(a)

10 6

n

(b)

Fig. 7. Performance of the solver for different problem sizes using different levels of precision
(shown by different colors and symbols) in low-rank approximations. (a) Time spent on factorization
(solid line) and solve (dashed line) parts, (b) error (solid line) and residual (dashed line) of the
solution.

same size (corresponding to Voronoi tessellation) are shown. The error and residual
decrease proportional to the precision in the low-rank approximations, . Furthermore,
it is clear that the residual is smaller than the error. The ratio of the error to residual
is generally an increasing function with respect to the condition number of the matrix.
For the same number of points, the condition number of a 3D discretization is lower
than that of a 2D discretization. Therefore, the error and residual are closer in the
3D case in comparison to the 2D case as illustrated in Figure 8b.
Figure 8a demonstrates the factorization time as a function of the low-rank approximation precision, . Three dimensionality leads to a higher number of wellseparated interactions; hence, after every elimination more new fill-ins are introduced
in the 3D case. This results in a higher factorization time compared to the 2D case.
Figure 9 shows the breakdown of the time spent on different parts of the algorithm,
namely, low-rank approximation of well-separated blocks (here, SVD), general matrix
multiplication (gemm), and computing the inverse of pivot blocks. Clearly, for the 3D
case most of the time is spent on SVD, which is known to be an expensive algorithm
for low-rank approximation. This can be improved significantly if faster low-rank
approximation methods are employed. We discuss some of the alternative methods
in section 7.
In Figure 10, the average ranks of the well-separated interactions at each level are
depicted as a function of the low-rank approximation precision. Note that for both the
3D and 2D cases an H-tree with depth l = 13 is used, i.e., there are 16 variables per
leaf red-nodes on average. The rank for the 3D case increases dramatically compared
to the 2D case when a more accurate solution is desired. If dL is the maximum
rank among all levels (i.e., maximum size of a red-node), similar to the analysis of
2
Theorem 10, the factorization
 complexity is O(ndL ). From Figures 6 and 10 we can
1
observe that dL = O log  . Therefore, we have:
(25)



1
factorization cost = O n log2
,




1
factorization memory = O n log


19

FAST HIERARCHICAL SPARSE SOLVERS

10 4

10 0

factorizatiom time (sec)

2D
3D
10

10 -1

3

10 -2
10 -3

10 2
10 -4
10 -5

10 1

10 -6
0

10 0
10

10 -1

10 -2

10 -3

10

10 -4

error 2D
residual 2D
error 3D
residual 3D

slope = 1

-7

10 0

10 -1

10 -2

ǫ

ǫ

(a)

(b)

10 -3

10 -4

Fig. 8. (a) Factorization time as a function of the low-rank approximation precision () for
the 3D and 2D Poisson problems of size 1.3 × 105 ; (b) error and residual of the solution.

4.5
4.0
3.5

1200

svd
gemm
inversion

1000
800

2.5

time (sec)

time (sec)

3.0
2.0
1.5

600
400

1.0

200

0.5
0.0

svd
gemm
inversion

4096 2048 1024 512 256 128 464 32 16
10 ²

(a)

8

4

2

1

0

4096 2048 1024 512 256 128 464 32 16
10 ²

8

4

2

1

(b)

Fig. 9. Breakdown of the time spent on different parts as a function of the low-rank approximation precision for the (a) 2D and (b) 3D cases.

At the end of this section, we present another benchmark in which LoRaSp is
used as a stand-alone solver to obtain solutions with floating point single precision
accuracy. We consider (22) discretized on uniform 2D grids of different sizes, and use
the proposed solver with  = 10−4 for low-rank approximation as introduced in (23).
For different problem sizes we pick depth of the H-tree such that the average size of
a super-node at leaf level is 64. The factorization and solve times as functions of the
matrix size are depicted in Figure 11a demonstrating linear complexity of the solver.
The relative residual of the solution, as defined in (24) in all cases is less than 10−6
(see Figure 11b).
6.2. LoRaSp as a preconditioner. In subsection 6.1 we showed that for a
fixed low-rank approximation precision, and therefore solution accuracy, the total
cost of the algorithm grows linearly with the problem size. However, as suggested
by Figures 8 and 9 obtaining a high accuracy solution may be expensive for some
problems. One standard remedy in that case is to use the low-accuracy solver as a
high-accuracy preconditioner in conjunction with an iterative solver. This is particularly very appealing here, since the factorization part is completely separated form the

20

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE
20

10
5

200
150
rank

15
rank

250

level = 3
level = 4
level = 5
level = 6
level = 7
level = 8
level = 9
level = 10
level = 11
level = 12
level = 13

100

level = 3
level = 4
level = 5
level = 6
level = 7
level = 8
level = 9
level = 10
level = 11
level = 12
level = 13

50

00
10

10-1

10-2

²

10-3

00
10

10-4

(a)

10-1

10-2

²

10-3

10-4

(b)

Fig. 10. Average rank of well-separated interactions per level as a function of the low-rank
approximation precision for the (a) 2D and (b) 3D cases.

10 2

relative residual

time (sec)

10 1

10 0

factorization time
solve time
linear line

10 0

10 -1

10 -2 4
10

10 5

10 6

10 -2

10 -4

10 -6

10 -8 4
10

10 5

n

n

(a)

(b)

10 6

Fig. 11. LoRaSp factorization and solve times as stand-alone solver for discretized Poisson
equation on grids with different sizes (a), and the relative residual of the solution (b). Low-rank
precision is  = 10−4 , resulting in relative residual less than 10−6 and relative error ∼ 10−4 .

solve part. Therefore, we can factorize the matrix only once, and apply the (cheap)
solve part at every iteration. Here, we use the GMRES method [53] as the iterative
solver in conjunction with the proposed algorithm as a preconditioner for all benchmarks. Note that for matrices with specific properties such as a symmetric positive
definite (SPD) matrix, one can use a more optimized iterative method (e.g., conjugate
gradient in the case of SPD matrix).
6.2.1. Poisson equation (structured grid). As the first benchmark, we consider the sequence of 3D Poisson problems introduced in subsection 6.1. We use
 = 10−1 to factorize the matrix, and find an approximation Ã of A−1 . Factorization
and solve times are shown in Figure 7a. We solve the system of equation ÃAx = Ãb
through GMRES afterwards. Since the solve part of the proposed algorithm is much
cheaper compared to the factorization part, each GMRES iteration is also relatively
cheap. In Figure 12 the sparsity pattern of the original and preconditioned matrices are shown. The preconditioned matrix, ÃA, approaches the identity matrix as 
decreases.

21

FAST HIERARCHICAL SPARSE SOLVERS

7

6

5

4

3

2

1

0

(a)

(b)

Fig. 12. (a) Sparsity pattern of a matrix A obtained from a discretization of (22). (b) Sparsity
pattern of the preconditioned matrix ÃA using LoRaSp with  = 10−1 . A non-zero entry a is colored
by − log |a| (i.e., larger values are darker).

In Figure 13, the GMRES residual as a function of the iteration number is plotted
for different problem sizes, when LoRaSp with  = 10−1 is used as a preconditioner.
At every iteration of the preconditioned GMRES method we have an approximation
x̃ of the solution. The residual in this case is defined as follows:
(26)

preconditioned GMRES residual =

kÃb − ÃAx̃k2
kÃbk2

The number of iterations that GMRES needs to converge slightly grows with size of
the problem. This is due to growth of the condition number of the problem. Similar
to multi-grid methods, we can use specific knowledge about the underlying PDE
and discretization to obtain problem-size independent convergence. This is the topic
of our future work, and is not discussed in this paper. In Figure 14 the condition
number, κ(A), is plotted as a function of the matrix size. The condition numbers are
approximated using the 1-norm [34, 37]. Note that for a matrix of size n corresponding
to the second order finite difference discretization of the Poisson equation, we expect
the condition number to grow as n2/3 . The n2/3 trend is also depicted in Figure 14.
6.2.2. Variable coefficient Poisson equation (structured grid). As our
next benchmark we consider the variable coefficient Poisson equation with periodic
boundary conditions discretized on a three-dimensional uniform grid:
(27)

∇ · (φ∇T ) = f

In the above equation, the scalar fields φ and f are given, and we solve for T , similar
to (22). We consider three cases for the coefficient field, φ:
• case 1: At each point of the domain, φ is drawn from a uniform distribution,
unif(0, 1), independently.
• case 2: At each point of the domain, ρ is drawn from a uniform distribution,
unif(0, 1), independently. φ is then defined as φ = ρ1 .
• case 3: At each point of the domain, φ is drawn from a uniform distribution,
unif(−1, 1), independently.
For the two first cases, the corresponding matrices are symmetric negative definite.
Case 2 shows up in the numerical simulation of a variable-density flow in the low-Mach

22

residual

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE
10

1

10

2

10

3

10

4

10

5

10

6

10

7

10

8

10

9

10

10

10

11

10

12

10

13

10

14

10

15

0

n = 16k
n = 32k
n = 64k
n = 128k
n = 256k
n = 512k
n = 1024k

5

10

15
20
iteration number

25

30

35

Fig. 13. GMRES residual as a function of the iteration number for different problem sizes.
LoRaSp is used as a preconditioner with  = 10−1 .

4

10

3

κ(A)

10

2

10

1-nor m es timation of κ(A)
n 2 / 3 line
1

10 3
10

4

10

5

n

10

6

10

Fig. 14. Condition number versus matrix size for a sequence of matrices obtained from the
discretization of (22).

number limit [16, 30, 50], where in that case T and ρ are hydrodynamic pressure and
density of the flow, respectively. The third case, however, results in an indefinite
matrix. In Figure 15a, the norm of the eigenvalues of the matrices corresponding to
a 163 grid for all cases are shown. In the third case, nearly half of the eigenvalues are
positive, and half of them are negative, corresponding to the left and right sides of
the red curve in Figure 15a.
Figure 15b shows the 1-norm approximation of the condition number of the matrices for all cases. Evidently, a larger grid results in a higher condition number. Also,
as expected, the condition number in case 2 is higher than case 1, and the condition
number in case 3 is higher than case 2.
We used LoRaSp as a preconditioner in conjunction with GMRES. A summary of
the results for the two first cases is provided in Table 1. The convergence criterion of

23

FAST HIERARCHICAL SPARSE SOLVERS

107

106

10

case 1
case 2
case 3

4

κ(A)

|λi |

106
102

105

100

case 1
case 2

10

case 3

-2

0

1000

2000

3000

104

4000

i

size = 4k

(a)

size = 32k

(b)

Fig. 15. Properties of the matrices corresponding to the discretization of (27): (a) Absolute
value of the eigenvalues for the 163 grid. (b) The 1-norm approximation of the condition number
for 163 and 323 grids.

GMRES, with residual defined in (26), is set to 10−14 . For cases 1 and 2, we used LoRaSp with low-rank precision  = 10−1 as defined in (23). Similar to subsection 6.2.1,
we increase the depth of the H-tree with log n, where n is the size of matrix. Similar
to the results presented in subsection 6.2.1, the factorization time has an almost linear complexity with the size of the matrix. As depicted in Figure 15b, the condition
number grows rapidly from case 1 to case 2, and with the size of the matrix. This
explains a slight growth of the number of iterations, and relative error.
case
1
1
1
2
2
2

(4k)
(32k)
(256k)
(4k)
(32k)
(256k)

fact. time

GMRES time

tot. time

# iters

rel. error

0.44
4.62
37.08
0.40
3.72
33.52

0.06
1.07
16.21
0.07
1.31
19.44

0.50
5.69
53.29
0.47
5.03
52.95

10
15
15
12
20
30

1.7e-11
2.6e-10
8.6e-10
6.6e-11
9.2e-10
1.2e-9

Table 1
GMRES performance using LoRaSp as a preconditioner to solve a variable coefficient Poisson
equation. Matrix sizes, corresponding to 163 , 323 , and 643 grids are written in parentheses. Times
are reported in seconds.

Case 3 corresponds to an indefinite matrix, with a large condition number. This
is typically a more difficult problem, compared to the first two cases. Since case 3 is
inherently a harder problem compared to cases 1 and 2, we choose a higher low-rank
precision  = 10−3 . In Figure 16a, the averaged rank of interactions for each level is
plotted. We also plot the compression ratio for each level, which is defined as follows:
(28)

compression ratio (l) =

hinteraction rankil
,
hsize of super-nodesil

where h il denotes averaging in level l of the H-tree. Note that it is clear from Figure 16a that even though the rank is increasing, the compression ratio is approximately
0.6–0.7 for all levels. Hence, the algorithm takes advantage of low-rank structures appropriately.

24

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE

In Figure 16b, the preconditioned GMRES residual defined in (26) is plotted as
a function of the iteration number. GMRES for case 3 does not converge when no
preconditioner or a diagonal preconditioner is used. We also applied ILU [52] as a
preconditioner for GMRES. We tried various (including very large) values for the fill
parameter in ILU. No convergence was obtained when ILU is used as a preconditioner.

103

rank
compression ratio

residual

102

101

100

10-1 11

10

9

8

7

level

6

5

4

3

100
10−1
10−2
10−3
10−4
10−5
10−6
10−7
10−8
10−9
10−10
10−11
10−12
10−13
10−14
10−15
0

no precond. (12%)
diagonal (154%)
ILU (318%)
H-tree (1.9e-5%)
100

200
300
iteration number

(a)

400

500

(b)

Fig. 16. Results corresponding to the case 3 for a 323 grid. (a) Averaged rank and compression
ratio in the H-tree. (b) GMRES residual as a function of the iteration number using various preconditioners. The relative error at the end of the iterations is shown in parentheses for each case.
Although the residuals may be small, the relative error for some of the methods is very large.

6.2.3. Elasticity equation (unstructured grid). Our next benchmark is obtained from an unstructured grid (see [24]) to solve the three-dimensional elasticity
equation:
(29)

(λ + µ)∇(∇.u) + µ∇2 u + F = 0

The matrix is symmetric with size n = 334, 956, and total of 10, 977, 198 non-zero
entries. This problem is significantly more difficult in comparison to the previous
benchmarks. We used various preconditioning strategies to evaluate the performance
of our proposed solver. In Table 2 a summary of the results is provided. We consider 10−12 to be the convergence criterion for the GMRES residual, and consider a
maximum of 500 iterations.
When no preconditioner is used, GMRES does not converge, and we obtain a
solution with 2% relative error after 500 iterations. Employing a diagonal preconditioner (i.e., ignore all non-diagonal entries of the matrix, and approximate A−1 with
the inverse of its diagonal part) also does not lead to convergence. A 1.4% relative
error after 500 GMRES iterations is obtained, which is an improvement in comparison
to the case without a preconditioner.
Next, we tried ILU as a preconditioner for GMRES. We used dual-threshold ILU
with drop tolerance fixed equal to the GMRES convergence precision, and varying
fill values (1, 2, etc.). ILU with fill value of 1 does not converge after 500 iterations;
however, for fill values greater than 1, convergence is obtained before 500 iterations.
Increasing the fill value leads to a higher factorization time.
Finally, we used the proposed algorithm as a preconditioner. For the low-rank
approximation, we used a variation of (23). Consider we want to find a low-rank
approximation of a block B with singular-value decomposition B = U SV | . We keep

25

FAST HIERARCHICAL SPARSE SOLVERS

the first k singular-values (and therefore, singular vectors) such that, k is the smallest
integer that:
kB − Uk Sk Vk| kF
<
kall levels belowkF

(30)

The subscript k in Uk and Vk means keeping only the first k columns, and in Sk means
keeping the first k singular-values. k kF refers to the Frobenius norm. Therefore,
kall levels belowkF refers to square root of the sum of Frobenius norm squared of
all blocks at the current level as well as the levels below. Both the above criteria
and the one in (23) work properly, and lead to the same conclusions; however, the
above method is slightly more efficient. For this benchmark, we used a sequence of
decreasing values for  in (30): 1 = 1024 × 10−7 , 2 = 256 × 10−7 , 3 = 64 × 10−7 , 4 =
16 × 10−7 , 5 = 4 × 10−7 , 6 = 1 × 10−7 .
precond.

fact. time

GMRES time

tot. time

# iters

rel. error

none
diagonal
ILU 1

0
0.02
98.6

115.9
116.3
156.7

115.9
116.3
255.3

500
500
500

2.1e-2
1.4e-2
5.3e-06

ILU 2
ILU 3
ILU 4
H-tree 1
H-tree 2
H-tree 3
H-tree 4
H-tree 5
H-tree 6

211.1
313.0
399.9
90.6
93.5
114.6
166.0
379.8
961.1

132.0
102.2
82.3
298.9
299.6
295.5
130.6
74.9
43.8

343.1
415.2
482.2
389.5
393.1
410.1
296.6
454.7
1004.9

408
309
245
467
466
367
144
46
16

4.9e-10
5.7e-10
3.2e-10
2.3e-10
2.2e-10
2.4e-10
3.5e-10
1.3e-10
2.1e-10

Table 2
GMRES performance using various preconditioners for a matrix of size 330k with more than 10
million non-zeros obtained from a 3D unstructured discretization of the elasticity equation. Times
are reported in seconds.

Figure 17 illustrates the variation of residual versus the number of iterations using various preconditioners. Clearly, diagonal preconditioner accelerates convergence
compared to the case with no preconditioner. ILU preconditioners bring about convergence faster than diagonal. The H-tree based preconditioners lead to a significant
acceleration in convergence. Decreasing  (and similarly increasing fill value in the
ILU) results in a shorter iteration time at the cost of a more expensive factorization
as listed in Table 2. In practice, one should pick an intermediate value for  (and
similarly for the fill value when ILU is used) to get an optimal total runtime (i.e.,
factorization + GMRES iterations).
In Figure 18 the breakdown of the total time is plotted for the cases that convergence is achieved. The non-monotonic functionality of the total time as a function of
 is clear form this figure. For instance, for this set of  values, 4 = 16 × 10−7 has the
optimal time, which is more efficient in comparison to the best time of the ILU. Note
that the current implementation of the algorithm is completely sequential. There are
various optimizations to enhance the performance of the solver. We discuss some of
the possible optimizations in section 7. Typically,  meets its optimal value when
factorization and iteration times are almost equal. This is also evident in Figure 18.

26

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE
100
no precond.
diagonal
ILU 1
ILU 2
ILU 3
ILU 4

10−1
10−2
10−3
10−4

H-tree 1
H-tree 2
H-tree 3
H-tree 4
H-tree 5
H-tree 6

residual

10−5
10−6
10−7
10−8
10−9
10−10
10−11
10−12
10−13
0

100

200
300
iteration number

400

500

Fig. 17. GMRES residual as a function of iteration number using various preconditioners
for a matrix of size 330k with more than 10 million non-zeros obtained from a 3D unstructured
discretization of the elasticity equation.
1200
factorization
GMRES
1000

time (sec)

800

600

400

200

0

ILU 2

ILU 3

ILU 4

H-tree 3

H-tree 4

H-tree 5

H-tree 6

Fig. 18. Total solve time using various preconditioners for a matrix of size 330k with more
than 10 million non-zeros obtained from a 3D unstructured discretization of the elasticity equation.

6.2.4. Advection-diffusion problem (non-symmetric linear system). In
our last benchmark, we consider the advection-diffusion problem in a cubic domain of
size L3 with Dirichlet boundary conditions. This problem is governed by the following
PDE.
(31)

∂T
+ U · ∇T = ν∇2 T + s
∂t

In the above equation U = (Ux , Uy , Uz ) is a constant velocity vector, ν is the diffusion
coefficient, and s is a source term. Using an implicit numerical time integration
method, and assuming Ux = Uy = Uz = U , the above equation transforms to the
following non-dimensional form
(32)

σT + R∇T − ∇2 T = g,

27

FAST HIERARCHICAL SPARSE SOLVERS
2

L
where σ = ν∆t
and R = LU
ν assuming ∆t is the time step.
Discretizing (32) using a central finite differencing method results in symmetric and skew-symmetric matrices for the diffusion and advection terms, respectively.
Therefore, the full system of equations is represented by a non-symmetric matrix.
We verified the accuracy of the H-tree solver (similar to Figure 8) for a case with
σ = 0 and R = 1 on a 323 grid. Figure 19a shows the error and residual as a function
of the low-rank precision parameter  as defined in (23). Similar to the symmetric
case, residual and error are proportional to .
Furthermore, we compare the convergence of the GMRES solver for the advectiondiffusion problem when H-tree and ILU are used as preconditioner. We consider a
sequence of problems with σ = 1 and varying R on a 323 grid. In Figure 19b the
numbers of GMRES iterations required to converge to a solution with residual less
than 10−10 are shown for different cases. For the H-tree preconditioner a low-rank
precision of  = 10−1 is used, while for the ILU drop tolerance is equal to the GMRES
convergence precision and fill parameter is set to 3 (the minimum fill parameter such
that all cases converge with less than 500 iterations). A 1-norm approximation of the
condition number of the system is also illustrated in Figure 19b. For the problem
studied here, the condition number (and therefore, the number of GMRES iterations)
is a non-monotonic function of R. H-tree preconditioner exhibits a stable number
of iterations and accuracy for all cases. The ILU preconditioner, however, gives rise
to a larger number of iterations and higher variation with R. Furthermore, for the
case with R = 1024 the ILU preconditioner results in a solution with final residual
(defined in (24)) of order 10−2 , while the preconditioned residual (defined in (26))
is less than 10−10 . Such discrepancy means the preconditioner (in this case ILU) is
ill-conditioned.

10-1

error
residual
linear line

10-2
10-3
10-4

10 4
condition number
H-tree iterations
ILU iterations
10 3

10 2

10-5
10 1

10-6
10-7100

10-1

10-2
²

(a)

10-3

10-4

10 0 0
10

10 1

10 2

10 3

R

(b)

Fig. 19. (a) Error and residual of the solution of (32) as a function of low-rank precision 
using H-tree solver. (b) Number of GMRES iterations using ILU and H-tree as preconditioner. For
R = 1024, the ILU preconditioner results in a final residual 10−2 , while the pre-conditioned residual
is less than 10−10 . This case is highlighted by a red square in the figure. The 1-norm approximation
of the matrix condition number is also plotted.

7. Conclusion and future works. We proposed a new algorithm to solve
sparse linear systems with numerically low-rank structures in linear time. The algorithm is based on the LU factorization of the sparse matrix, where the matrices L

28

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE

and U are computed and stored using a hierarchical low-rank structure. The accuracy
of the factorization is determined a priori. For a precision
 tolerance , the complexities
of the factorization cost and memory are O n log2 1/ and O (n log 1/), respectively.
The proposed algorithm is fully algebraic and preserves the sparsity of the original matrix during the elimination. Therefore, it can be considered as an extension
to the ILU method. In the ILU factorization, new fill-ins are ignored. In the proposed algorithm, however, new fill-ins are compressed using low-rank approximations.
Compressed fill-ins form a new set of equations —at a coarser level— which are factorized through elimination. Furthermore, the multilevel process of the factorization
is similar to AMG, where the original system is solved at different levels (grid size).
We provided various benchmarks to illustrate the performance of the proposed
algorithm. We used matrices obtained from the discretization of the Poisson and elasticity equations on structured and unstructured grids, respectively. A non-symmetric
benchmark corresponding to the advection-diffusion problem is also presented. The
proposed factorization method is used both as a stand-alone solver with tunable accuracy, and as a preconditioner in conjunction with an iterative method (e.g., GMRES).
There are various aspects of the method which are general, and can be modified
to optimize the solver for particular matrices without loosing the properties demonstrated in this paper. Here is a list of some aspects that can be modified in the
algorithm:
1. Partitioning of the sparse matrix graph is generic. Higher quality partitioning
in general results in higher accuracy solution. If the matrix is associated to a
physical grid, the physical coordinates of the solution points can be used to
improve the quality of partitioning.
2. Here, we used a binary tree corresponding to the recursive bi-partitioning of
the graph. Many other options are possible, e.g., using octree if the matrix
comes from a three-dimensional problem, or an adaptive tree with arbitrary
number of children per node.
3. We used SVD for low-rank representation of the well-separated nodes. Other
low-rank approximation methods could be used as well, e.g., randomized SVD
[35], randomized block algorithm [55], adaptive cross approximation [8], rankrevealing QR/LU [13], etc.
4. Having a low-rank representation of a well-separated block, there are different
measures to define the error. For instance, we used the ratio of the singular
values to the largest singular value as a measure of accuracy in the low-rank
approximation in (23). One can use different criteria, e.g., absolute singular
values, ratio of singular values to largest singular value of the whole level or
the full matrix, Frobenius norm of the low-rank block, etc.
5. We defined two super-nodes as well-separated if their distance is greater than
1 (see Corollary 8). One can change this definition, and make it stronger.
For example, define two super-nodes as well-separated if their distance is at
least 3. This is similar to the fill value parameter in the ILU.
6. At each level, we can use any ordering to eliminate super-nodes and blacknodes. Here we used a generic ordering. There are various other orderings
which reduce the calculation cost, including the minimum degree, minimum
deficiency, nested dissection, etc. [21]. The complexity of the algorithm
remains linear irrespective of the ordering. This is particularly an alluring
property for a parallel implementation of the algorithm.

29

FAST HIERARCHICAL SPARSE SOLVERS

Appendix A. In this section, we present an example of the factorization process
(see Algorithm 1) for one level. Each figure represents one step of the algorithm. The
H-tree is shown on the right, and the corresponding extended matrix is shown on the
left.
level 3

r 01
r 11
r 02
r 12
r 03
r 13
r 04
r 14

level 2

r 01 r 11 r 02 r 12 r 03 r 13 r 04 r 14 b 2 b 4 r 11 r 12

b2
b4
r 11
r 12

Fig. 20. Original matrix (left) and the corresponding adjacency graph (right).

level 3

s1

s2

s3

level 2

s4

b 2 b 4 r 11 r 12

s1
s2
s3
s4
b2
b4
r 11
r 12
[3]

Fig. 21. Creating super-nodes at level 3 (see subsection 4.2), where the super-node si
of red-nodes

[3]
r0 i

and

[3]
r1 i .

consists

30

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE

level 3

s1

s2

s3

level 2

s4

b 2 b 4 r 11 r 12

s1
s2
s3
s4
b2
b4
r 11
r 12
[3]

Fig. 22. Eliminating s1 (see subsection 4.4). Green edges (and their corresponding blocks
in the matrix) represent a numerically low-rank interaction between two well-separated nodes to be
compressed.

level 3

s1

s2

s3

level 2

s4

b 2 b 4 r 11 r 12

s1
s2
s3
s4
b2
b4
r 11
r 12
[3]

[3]

Fig. 23. Compressing the well-separated edge between s2 and s4 (see subsection 4.3). Blue
edges (and their corresponding blocks in the matrix) correspond to minus identity.

level 3

s1

s2

s3

level 2

s4

b 2 b 4 r 11 r 12

s1
s2
s3
s4
b2
b4
r 11
r 12
[3]

[3]

Fig. 24. Eliminating s2 and b2 (see subsection 4.4).

FAST HIERARCHICAL SPARSE SOLVERS

level 3

s1

s2

s3

level 2

s4

b 2 b 4 r 11 r 12

s1
s2
s3
s4
b2
b4
r 11
r 12
[3]

Fig. 25. Eliminating s3 (see subsection 4.4).

level 3

s1

s2

s3

level 2

s4

b 2 b 4 r 11 r 12

s1
s2
s3
s4
b2
b4
r 11
r 12
[2]

[3]

Fig. 26. Compressing the well-separated edge between r1 1 and s4 (see subsection 4.3).

level 3

s1

s2

s3

level 2

s4

b 2 b 4 r 11 r 12

s1
s2
s3
s4
b2
b4
r 11
r 12
[3]

[3]

Fig. 27. Eliminating s4 and b4 (see subsection 4.4).

31

32

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE

Appendix B. In this appendix we provide a graphical example of a nested partitioning using the SCOTCH library for a sparse matrix corresponding to discretization
of (22) on a 2D Voronoi grid.

Fig. 28. An example of 6 levels of a nested partitioning. Clusters are distinguished by different
colors. The edges between different clusters are intentionally omitted in this figure for visualization
purpose.

FAST HIERARCHICAL SPARSE SOLVERS

33

REFERENCES
[1] Noga Alon and Raphael Yuster, Solving linear systems through nested dissection, in Foundations of Computer Science (FOCS), 2010 51st Annual IEEE Symposium on, IEEE, 2010,
pp. 225–234.
[2] Sivaram Ambikasaran, Fast Algorithms for Dense Numerical Linear Algebra and Applications, PhD thesis, Stanford University, 2013.
[3] Sivaram Ambikasaran and Eric Darve, An O(n log n) fast direct solver for partial hierarchically semi-separable matrices, Journal of Scientific Computing, 57 (2013), pp. 477–501.
[4] Sivaram Ambikasaran and Eric Darve, The inverse fast multipole method, arXiv preprint
arXiv:1407.1572, (2014).
[5] Amirhossein Aminfar, Sivaram Ambikasaran, and Eric Darve, A fast block low-rank dense
solver with applications to finite-element matrices, Journal of Computational Physics, 304
(2016), pp. 170–188.
[6] Edward Anderson, Zhaojun Bai, Christian Bischof, Susan Blackford, James Demmel,
Jack Dongarra, Jeremy Du Croz, Anne Greenbaum, S Hammerling, Alan McKenney, et al., LAPACK Users’ guide, vol. 9, Siam, 1999.
[7] Mario Bebendorf, Hierarchical matrices, Springer, 2008.
[8] Mario Bebendorf and Sergej Rjasanow, Adaptive low-rank approximation of collocation
matrices, Computing, 70 (2003), pp. 1–24.
[9] Steffen Börm, Lars Grasedyck, and Wolfgang Hackbusch, Introduction to hierarchical matrices with applications, Engineering Analysis with Boundary Elements, 27 (2003),
pp. 405–422.
[10] Achi Brandt, Algebraic multigrid theory: The symmetric case, Applied mathematics and
computation, 19 (1986), pp. 23–56.
[11] A Brandt, S McCormick, and J Huge, Algebraic multigrid (AMG) for sparse matrix equations, Sparsity and its Applications, (1985), p. 257.
[12] James R Bunch and John E Hopcroft, Triangular factorization and inversion by fast matrix
multiplication, Mathematics of Computation, 28 (1974), pp. 231–236.
[13] Tony F Chan, Rank revealing QR factorizations, Linear algebra and its applications, 88 (1987),
pp. 67–82.
[14] Shiv Chandrasekaran, Patrick Dewilde, Ming Gu, William Lyons, and Timothy Pals,
A fast solver for hss representations via sparse matrices, SIAM Journal on Matrix Analysis
and Applications, 29 (2006), pp. 67–81.
[15] Shiv Chandrasekaran, Ming Gu, and Timothy Pals, A fast ULV decomposition solver
for hierarchically semiseparable representations, SIAM Journal on Matrix Analysis and
Applications, 28 (2006), pp. 603–622.
[16] Y-H Choi and Charles L Merkle, The application of preconditioning in viscous flows, Journal of Computational Physics, 105 (1993), pp. 207–223.
[17] Don Coppersmith and Shmuel Winograd, Matrix multiplication via arithmetic progressions,
in Proceedings of the nineteenth annual ACM symposium on Theory of computing, ACM,
1987, pp. 1–6.
[18] Eduardo Corona, Per-Gunnar Martinsson, and Denis Zorin, An O(n) direct solver for
integral equations on the plane, Applied and Computational Harmonic Analysis, 38 (2015),
pp. 284–317.
[19] Pieter Coulier, Hadi Pouransari, and Eric Darve, The inverse fast multipole method:
using a fast approximate direct solver as a preconditioner for dense linear systems, arXiv
preprint arXiv:1508.01835, (2015).
[20] Eric Darve, The fast multipole method: numerical implementation, Journal of Computational
Physics, 160 (2000), pp. 195–240.
[21] Timothy A Davis, Direct methods for sparse linear systems, vol. 2, Siam, 2006.
[22] William Fong and Eric Darve, The black-box fast multipole method, Journal of Computational Physics, 228 (2009), pp. 8712–8725.
[23] Alan George, Nested dissection of a regular finite element mesh, SIAM Journal on Numerical
Analysis, 10 (1973), pp. 345–363.
[24] Debraj Ghosh, Philip Avery, and Charbel Farhat, A FETI-preconditioned conjugate gradient method for large-scale stochastic finite element problems, International journal for
numerical methods in engineering, 80 (2009), pp. 914–931.
[25] Adrianna Gillman and Per-Gunnar Martinsson, A direct solver with O(n) complexity for
variable coefficient elliptic pdes discretized via a high-order composite spectral collocation
method, SIAM Journal on Scientific Computing, 36 (2014), pp. A2023–A2046.
[26] Adrianna Gillman, Patrick M Young, and Per-Gunnar Martinsson, A direct solver with

34

[27]

[28]
[29]
[30]
[31]
[32]
[33]

[34]
[35]

[36]
[37]

[38]
[39]
[40]
[41]

[42]
[43]
[44]
[45]
[46]
[47]
[48]

[49]
[50]

[51]
[52]
[53]

HADI POURANSARI, PIETER COULIER, AND ERIC DARVE
O(n) complexity for integral equations on one-dimensional domains, Frontiers of Mathematics in China, 7 (2012), pp. 217–247.
Leslie Greengard, Denis Gueyffier, Per-Gunnar Martinsson, and Vladimir Rokhlin,
Fast direct solvers for integral equations in complex three-dimensional domains, Acta Numerica, 18 (2009), pp. 243–275.
Leslie Greengard and Vladimir Rokhlin, A fast algorithm for particle simulations, Journal
of computational physics, 73 (1987), pp. 325–348.
Gaël Guennebaud, Benoı̂t Jacob, et al., Eigen v3. http://eigen.tuxfamily.org, 2010.
Hervé Guillard and Cécile Viozat, On the behaviour of upwind schemes in the low mach
number limit, Computers & fluids, 28 (1999), pp. 63–86.
Wolfgang Hackbusch, A sparse matrix arithmetic based on H-matrices. part I: Introduction
to h-matrices, Computing, 62 (1999), pp. 89–108.
Wolfgang Hackbusch and Steffen Börm, H2-matrix approximation of integral operators
by interpolation, Applied Numerical Mathematics, 43 (2002), pp. 129–143.
Wolfgang Hackbusch and Boris N Khoromskij, A sparse H-matrix arithmetic: general
complexity estimates, Journal of Computational and Applied Mathematics, 125 (2000),
pp. 479–501.
William W Hager, Condition estimates, SIAM Journal on Scientific and Statistical Computing, 5 (1984), pp. 311–316.
Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp, Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions,
SIAM review, 53 (2011), pp. 217–288.
Magnus Rudolph Hestenes and Eduard Stiefel, Methods of conjugate gradients for solving
linear systems, (1952).
Nicholas J Higham and Françoise Tisseur, A block algorithm for matrix 1-norm estimation, with an application to 1-norm pseudospectra, SIAM Journal on Matrix Analysis and
Applications, 21 (2000), pp. 1185–1201.
Kenneth L Ho and Lexing Ying, Hierarchical interpolative factorization for elliptic operators: differential equations, Communications on Pure and Applied Mathematics, (2015).
Ho, Kenneth L and Ying, Lexing, Hierarchical interpolative factorization for elliptic operators: integral equations, Communications on Pure and Applied Mathematics, (2015).
Oscar H Ibarra, Shlomo Moran, and Roger Hui, A generalization of the fast LUP matrix
decomposition algorithm and applications, Journal of Algorithms, 3 (1982), pp. 45–56.
Wai Yip Kong, James Bremer, and Vladimir Rokhlin, An adaptive fast direct solver for
boundary integral equations in two dimensions, Applied and Computational Harmonic
Analysis, 31 (2011), pp. 346–369.
Ruipeng Li and Yousef Saad, Divide and conquer low-rank preconditioners for symmetric
matrices, SIAM Journal on Scientific Computing, 35 (2013), pp. A2069–A2095.
Richard J Lipton, Donald J Rose, and Robert Endre Tarjan, Generalized nested dissection, SIAM journal on numerical analysis, 16 (1979), pp. 346–358.
Artem Napov and Xiaoye S Li, An algebraic multifrontal preconditioner that exploits the
low-rank property, Numerical Linear Algebra with Applications, 23 (2016), pp. 61–82.
Naoshi Nishimura, Fast multipole accelerated boundary integral equation methods, Applied
Mechanics Reviews, 55 (2002), pp. 299–324.
Ivan V Oseledets and SV Dolgov, Solution of linear systems and matrix inversion in the
TT-format, SIAM Journal on Scientific Computing, 34 (2012), pp. A2718–A2739.
Christopher C Paige and Michael A Saunders, Solution of sparse indefinite systems of
linear equations, SIAM journal on numerical analysis, 12 (1975), pp. 617–629.
François Pellegrini and Jean Roman, Scotch: A software package for static mapping by
dual recursive bipartitioning of process and architecture graphs, in High-Performance Computing and Networking, Springer, 1996, pp. 493–498.
Hadi Pouransari and Eric Darve, Optimizing the adaptive fast multipole method for fractal
sets, SIAM Journal on Scientific Computing, 37 (2015), pp. A1040–A1066.
Hadi Pouransari, Milad Mortazavi, and Ali Mani, Parallel variable-density particle-laden
turbulence simulation, Annual Research Briefs, Center for Turbulence Research, (2015),
pp. 43–54.
JW Ruge and Klaus Stüben, Algebraic multigrid, Multigrid methods, 3 (1987), pp. 73–130.
Yousef Saad, ILUT: A dual threshold incomplete LU factorization, Numerical linear algebra
with applications, 1 (1994), pp. 387–402.
Yousef Saad and Martin H Schultz, GMRES: A generalized minimal residual algorithm for
solving nonsymmetric linear systems, SIAM Journal on scientific and statistical computing,
7 (1986), pp. 856–869.

FAST HIERARCHICAL SPARSE SOLVERS

35

[54] Klaus Stüben, A review of algebraic multigrid, Journal of Computational and Applied Mathematics, 128 (2001), pp. 281–309.
[55] Sergey Voronin and Per-Gunnar Martinsson, A randomized blocked algorithm
for efficiently computing rank-revealing factorizations of matrices, arXiv preprint
arXiv:1503.07157, (2015).
[56] Jianlin Xia, Shivkumar Chandrasekaran, Ming Gu, and Xiaoye S Li, Fast algorithms
for hierarchically semiseparable matrices, Numerical Linear Algebra with Applications, 17
(2010), pp. 953–976.
[57] Mihalis Yannakakis, Computing the minimum fill-in is NP-complete, SIAM Journal on Algebraic Discrete Methods, 2 (1981), pp. 77–79.
[58] Lexing Ying, George Biros, and Denis Zorin, A kernel-independent adaptive fast multipole
algorithm in two and three dimensions, Journal of Computational Physics, 196 (2004),
pp. 591–626.

