Simple non-extensive sparsification of the hierarchical matrices

TL;DR

This paper introduces a non-extensive sparse factorization method for H2 matrices, enabling efficient solutions comparable to other hierarchical solvers, with potential improvements in computational efficiency.

Contribution

It presents a novel sparse factorization approach for H2 matrices that replaces direct solutions with sparse system solutions, constructed from H2 matrix parameters.

Findings

Consistent results compared to HODLR, H2Lib, and IFMM solvers

Efficient approximation of H2 matrices with sparse factors

Potential reduction in computational complexity

Abstract

In this paper, we consider the matrices approximated in H2 format. The direct solution, as well as the preconditioning, of systems with such matrices is a challenging problem. We propose a non-extensive sparse factorization of the H2 matrix that allows to substitute the direct H2 solution with the solution of the system with an equivalent sparse matrix of the same size. The sparse factorization is constructed out of parameters of the H2 matrix. In the numerical experiments, we show the consistency of this approach in comparison to the other approximate block low-rank hierarchical solvers, such as HODLR, H2Lib and IFMM.

Tables

Table 1. Table 1 : Comparison with H2Lib.

N	3072	12288	49152	196608
H2Lib approx., sec	6.72	24.26	104.97	487.24
H2Lib factor., sec	13.56	164.00	1712.17	10970.43
Sp. approx., sec	4.8	26.34	110.91	399.43
Sp. sp., sec	2.16	10.12	63.47	323.23
Sp. factor., sec	0.19	1.30	7.83	56.89

Table 2. Table 2 : Comparison with IFMM method.

	$d={10}^{-3}$	$d={10}^{-2}$
IFMM solution time, sec	118	301
Sparsification solution time, sec	96	218

Equations

$$A=US{V}^{\top },$$

$$Ax=b,$$

$$Sy=U^{\top }b,$$

$$A=C_k+F_k,k\in 0,\dots ,L$$

$$B_k=2^{k}B,$$

$$F_{ij}\approx {\tilde{U}}_i{\tilde{F}}_{ij}{\tilde{V}}_j^{\top },\forall i,j\in 1,\dots ,M$$

$${\tilde{F}}_{ij}=\left[\begin{array}{cc}{\dot{F}}_{ij}& 0\\ 0& 0\end{array}\right],$$

$$U_0^{\top }=\left[\begin{array}{ccc}{\tilde{U}}_1^{\top }& 0& 0\\ 0& \ddots & 0\\ 0& 0& {\tilde{U}}_M^{\top }\end{array}\right].$$

$$V_0=\left[\begin{array}{ccc}{\tilde{V}}_1& 0& 0\\ 0& \ddots & 0\\ 0& 0& {\tilde{V}}_M\end{array}\right].$$

$$A_1=U_0^{\top }A{V}_0.$$

$$A_1=U_0^{\top }(C_0+F_0)V_0=U_0^{\top }{C}_0{V}_0+{\tilde{F}}_1,$$

$${\tilde{A}}_1=P_{r1}{A}_1{P}_{c1}$$

$${\tilde{A}}_1=\left[\begin{array}{cc}{A}_{{\mathbf{n}}_1{\mathbf{n}}_1}& A_{{\mathbf{n}}_1{\mathbf{b}}_1}\\ A_{{\mathbf{b}}_1{\mathbf{n}}_1}& A_{{\mathbf{b}}_1{\mathbf{b}}_1}\end{array}\right],$$

$${\hat{F}}_1=(P_{r1}{\tilde{F}}_1{P}_{c1}).$$

$${\hat{C}}_1=(P_{r1}{U}_0^{\top }{C}_0{V}_0{P}_{c1}).$$

$$C_1={\hat{C}}_1+{\hat{F}}_{\mathbf{ml}1}.$$

$$A_{{\mathbf{b}}_1{\mathbf{b}}_1}=({\hat{C}}_1{)}_{{\mathbf{b}}_1{\mathbf{b}}_1}+{\hat{F}}_1=({\hat{C}}_1{)}_{{\mathbf{b}}_1{\mathbf{b}}_1}+{\hat{F}}_{\mathbf{ml}1}+F_1=(C_1{)}_{{\mathbf{b}}_1{\mathbf{b}}_1}+F_1.$$

$${\tilde{A}}_1={\hat{C}}_1+{\hat{F}}_1={\hat{C}}_1+{\hat{F}}_{\mathbf{ml}1}+F_1=C_1+F_1.$$

$${\tilde{F}}_2=U_{{\mathbf{b}}_1}^{\top }{F}_1{V}_{{\mathbf{b}}_1},$$

$$U_1=\left[\begin{array}{cc}{I}_{(N-Mr)\times (N-Mr)}& 0\\ 0& U_{{\mathbf{b}}_1}\end{array}\right],V_1=\left[\begin{array}{cc}{I}_{(N-Mr)\times (N-Mr)}& 0\\ 0& V_{{\mathbf{b}}_1}\end{array}\right].$$

$$A_2=U_1^{\top }{\tilde{A}}_1{V}_1.$$

$$A_2=U_1^{\top }(C_1+F_1)V_1=U_1^{\top }{C}_1{V}_1+{\hat{F}}_2.$$

$$\begin{array}{c}{A}_1=U_0^{\top }A{V}_0=U_0^{\top }{C}_0{V}_0+{\hat{F}}_1\\ A_2=U_1^{\top }{U}_0^{\top }A{V}_0{V}_1=U_1^{\top }{C}_1{V}_1+{\hat{F}}_2\\ ⋮\\ A_L=\left(\prod _{k=0}^L{U}_k^{\top }\right)A\left(\prod _{k=L}^0{V}_k\right)=U_{L-1}^{\top }{C}_{L-1}{V}_{L-1}+{\hat{F}}_L=S,\end{array}$$

$$A=\left(\prod _{k=L}^0{U}_k\right)S\left(\prod _{k=0}^L{V}_k^{\top }\right).$$

$$U=\prod _{k=0}^L{U}_k,V=\prod _{k=0}^L{V}_k,$$

$$A=US{V}^{\top },$$

$$\text{bsp}:{\mathbb{R}}^{M_{1}B\times M_{2}B}\to {\mathbb{B}}^{M_1\times M_2},$$

$$\{\begin{array}{c}{R}_{ij}=1,\text{if}{A}_{ij}\in {\mathbb{R}}^{B\times B}\text{is nonzero block},\\ R_{ij}=0,\text{if}{A}_{ij}\in {\mathbb{R}}^{B\times B}\text{is zero block. }\end{array}$$

$$A=U^{\top }SV$$

$$U=\left(\prod _{j=0}^K{U}_j{P}_j^{\top }\right),V=\left(\prod _{j=0}^K{V}_j{P}_j^{\top }\right),$$

Full text

1 Introduction

Problems arising in the discretization of boundary integral equations (and a number of other problems with approximately separable kernels) lead to matrices that can be well-approximated by hierarchical block low-rank ($\mathcal{H}$[15, 17], mosaic skeleton[27]) matrices. These are the matrices hierarchically divided into blocks, some of which has low-rank. The development of the $\mathcal{H}$ matrices is the $\mathcal{H}^{2}$[16, 6] matrices, which are the hierarchical block low-rank matrices with nested bases. The nested basis property leads to the additional improvement in terms of storage and complexity of different operations such as matrix-vector products. Approximate solution and preconditioning of systems with $\mathcal{H}^{2}$ matrices is a rapidly developed area[9, 3, 11, 19], however, construction of the accurate, time and memory efficient factorization that leads to approximate solution is still a challenging problem. In this paper, we propose a new representation of $\mathcal{H}^{2}$ matrices. Namely, we show that $\mathcal{H}^{2}$ factorization of matrix $A\in\mathbb{R}^{N\times N}$ is equivalent to the factorization

[TABLE]

where $S\in\mathbb{R}^{N\times N}$ is a sparse matrix. Note that the size of matrix $S$ matches the size of matrix $A$. $U\in\mathbb{R}^{N\times N}$ and $V\in\mathbb{R}^{N\times N}$ are orthogonal matrices that are products of block-diagonal and permutation matrices. Once the factorization (1) is built, we can substitute a solution of the system

[TABLE]

by a solution of the system with the sparse matrix:

[TABLE]

where $x=Vy$. The system (2) can be easily solved using standard sparse tools. In this paper we propose:

• •

Sparse non-extensive111The term non-extensive means that the sizes of the factors $S$, $U$ and $V$ are equal to the size of the matrix $A$ (in opposition to extensive[2, 11, 26] sparse factorizations of $\mathcal{H}^{2}$ matrix). factorization for $\mathcal{H}^{2}$ matrix that leads to the solver and the preconditioner.

• •

The algorithm that allows to construct factors $U$, $S$ and $V$ from parameters of $\mathcal{H}^{2}$

matrix.

• •

Numerical comparison of the proposed method with HODLR[1], IFMM[11] and H2Lib[5] packages.

The main idea of the sparsification algorithm is the compression of the low-rank blocks (the very close idea of the compression of the fill-in blocks during the block Cholesky factorization of a sparse matrix is presented in works[25, 29]). The main difference between the presented sparse factorization and the other methods of sparsification[2, 11, 26] is a size of factors. In the other methods, the hierarchical matrix is factorized into the sparse matrices of larger sizes. The characteristic inflating coefficient is $k=5$ (it depends on the number of levels in the cluster tree[2, 26]). The matrix extension is a major drawback since it increases the complexity of matrix computations. We propose the factorization that takes $\mathcal{H}^{2}$ matrix and returns the sparse factors of the same size. Another drawback of the extended sparse factorizations is that resulting sparse matrix may lose a positive definiteness of the original $\mathcal{H}^{2}$ matrix. Proposed sparsification preserves symmetry and positive definiteness of the $\mathcal{H}^{2}$ matrix.

2 Compression algorithm

Consider the dense matrix $A\in\mathbb{R}^{N\times N}$ that can be approximated in $\mathcal{H}^{2}$ format (has corresponding low-rank blocks). Formal definition of the $\mathcal{H}^{2}$ matrix is presented in Section 4.1, here we give basic facts that are used in the current section. Matrix $A$ is a block matrix with following properties. It consists of two non-intersecting “close” and “far” matrices:

[TABLE]

block size of zero level is $B$, block size of $k$-th level $B_{k}$ is

[TABLE]

$C_{k}\in\mathbb{R}^{N\times N}$ is block-sparse matrix of full-rank close blocks, and $F_{k}\in\mathbb{R}^{N\times N}$ is block matrix of far blocks, see Figure 1(b). The matrix $F_{k}$ has low-rank block rows and block columns. Moreover, the nested basis property holds: basis rows for block rows and columns on level $l$ are a subset of basis rows at level $(l-1)$. This is used in the multilevel computations, which are described in Section 2.2.

2.1 Compression at zero level

First consider the compression procedure at zero block level ($l=0$). Assume that the number of block rows and columns at zero level is $M$. Nonzero block $F_{ij}\in\mathbb{R}^{B\times B}$ of far matrix $F_{0}$ has low rank:

[TABLE]

where $\widetilde{F}_{ij}\in\mathbb{R}^{N\times N}$ is the compressed far block with the following structure:

[TABLE]

where $\dot{F}_{ij}\in\mathbb{R}^{r\times r}$. Matrices $\widetilde{U}_{i}\in\mathbb{R}^{N\times N}$ and $\widetilde{V}_{j}\in\mathbb{R}^{N\times N}$ are orthogonal. The blocks in $i$-th row have the same left orthogonal compression factor $\widetilde{U}_{i}$ and all blocks in $j$-th column have the same right factor $\widetilde{V}_{j}^{\top}$.

The goal of the compression procedure is to sparsify the matrix $A$ by obtaining the compressed blocks $\widetilde{F}_{ij}$ instead of original blocks $F_{ij}$. One can achieve this by finding $\widetilde{U}_{i}$ and $\widetilde{V}_{j}$ compression matrices and applying matrix $\widetilde{U}_{i}^{\top}$ to $i$-th row and $\widetilde{V}_{j}$ to $j$-th column. We introduce the block-diagonal orthogonal compression matrix

[TABLE]

Similarly, for block columns we obtain the block-diagonal orthogonal compression matrix

[TABLE]

Applying matrices $U_{0}$ and $V_{0}$ to the matrix $A$ we obtain the matrix $A_{1}$ with compressed far matrix:

[TABLE]

The process is illustrated in Figure 2.

Finally, we obtain:

[TABLE]

where $\widetilde{F}_{1}$ is a compressed far matrix which consists of blocks $\widetilde{F}_{ij}$. Note that the matrix $C_{0}$ is available as one of the parameters of the $\mathcal{H}^{2}$ format, it is a so-called “close matrix”.

2.2 Compression at the first level ($l=1$)

For each block row in $A_{1}$ we denote the rows with zero far blocks by “non-basis”, and the other rows by “first level basis”. Assume that each block row (column) has $r$ basis rows (columns) and $(B-r)$ non-basis. Introduce the permutation $P_{r1}$ that puts non-basis block rows before the basis ones preserving the row order and permutation $P_{c1}$ that does the same for columns, see Figure 4(a). For the permuted matrix

[TABLE]

we obtain

[TABLE]

where $A_{\mathbf{n_{1}n_{1}}}\in\mathbb{R}^{M(B-r)\times M(B-r)}$ is a submatrix on the intersection of non-basis rows and non-basis columns, $A_{\mathbf{b_{1}n_{1}}}\in\mathbb{R}^{Mr\times M(B-r)}$ is on the intersection of basis rows and non-basis columns and so on, see Figure 4(a). Denote the permuted far matrix:

[TABLE]

Note that permutations $P_{r1}$ and $P_{c1}$ concentrate all nonzero blocks of compressed far zone $\widetilde{F}_{1}$ inside of the submatrix $A_{\mathbf{b_{1}b_{1}}}$. Denote permuted close matrix:

[TABLE]

Consider the submatrix $A_{\mathbf{b_{1}b_{1}}}\in\mathbb{R}^{Mr\times Mr}$, note that this matrix has exactly the same close and far block structure as the matrix $A$, but the block size in $A_{\mathbf{b_{1}b_{1}}}$ is $r$. Now we join block rows and columns of the matrix $A_{\mathbf{b_{1}b_{1}}}$ by groups of $J$ blocks (e.g. $J=2$ in Figure 3(a) ). Assume that $Jr=B$.

We will call the grouped blocks “big blocks”. Among these blocks, the big block that consists only of far sub-blocks will be called far, the big block that contains at least one close small block will be referred to as close. Denote blocks of the far matrix $\widehat{F}_{1}$ that become close after grouping by $\widehat{F}_{\mathbf{ml}1}$, see Figure 3. We also introduce a new close matrix with big blocks by

[TABLE]

Denote far matrix with big blocks by $F_{1}$. Consider this joining for the block $A_{\mathbf{b_{1}b_{1}}}$:

[TABLE]

Similarly, for the matrix $\widetilde{A}_{1}$:

[TABLE]

It can be shown that block rows and columns of the matrix $F_{1}$ have low-rank by the properties of the $\mathcal{H}^{2}$ matrix $A$. Similarly to (3) compute orthogonal block-diagonal matrices $\\ U_{\mathbf{b_{1}}},V_{\mathbf{b_{1}}}\in\mathbb{R}^{Mr\times Mr}$ that compress matrix $F_{1}$.

Multiplication of matrix $F_{1}$ by matrices $U_{\mathbf{b_{1}}}$ and $V_{\mathbf{b_{1}}}$ leads to compression:

[TABLE]

where the matrix $\widetilde{F}_{2}$ consists of compressed blocks.

Now we introduce extended matrices $U_{\mathbf{b_{1}}}$ and $V_{\mathbf{b_{1}}}$ that can be applied to matrix $\widehat{A}_{1}$:

[TABLE]

Applying matrices $U_{1}$ and $V_{1}$ to matrix $\widetilde{A}_{1}$ we obtain the matrix with compressed first level:

[TABLE]

The process of the first level compression is shown in Figure 4(b).

For the first level we obtain

[TABLE]

2.3 Compression at all levels

We apply permutation and repeat this procedure $L$ times and obtain:

[TABLE]

thus

[TABLE]

If we denote

[TABLE]

then the final result of the algorithm is a sparse approximate factorization

[TABLE]

where $S$ is a sparse matrix of the same size as matrix $A$, $U$ and $V$ are orthogonal matrices that are products of permutation and block-diagonal orthogonal matrices.

Remark 2.1**.**

If matrix $A$ is approximated into $\mathcal{H}^{2}$ format, then matrices $S$, $U$ and $V$ can be constructed from parameters of the $\mathcal{H}^{2}$ representation, see details in Section 4.

Remark 2.2**.**

Sparsity of the matrix $S$ is proven in Section 3.

Proposition 1**.**

If the matrix $A$ is symmetric and positive definite, then the factors $U$ and $V$ are equal and the matrix $S$ is symmetric and positive definite.

Proof.

If the matrix $A$ is symmetric and positive definite, then from the steps of compression algorithm, compression matrices $U_{i}$ and $V_{i}$, $i\in 0,\dots,L$ are equal, thus, by equations (10), $U=V$. Since $S=U^{\top}AV$, $U=V$ and $A$ is symmetric and positive definite, then $S$ is symmetric and positive definite matrix. ∎

Remark 2.3**.**

The proposed sparsification algorithm is applicable to the special cases of $\mathcal{H}^{2}$ matrices such as HSS (Hierarchically Semiseparable), HOLDR (Hierarchical Off-Diagonal Low-Rank).

2.4 Pseudo code of the compression algorithm

3 Sparsity of the matrix $S$

First, define the block sparsity pattern of a block sparse matrix. For a matrix $A$ with $M_{1}$ block columns, $M_{2}$ block rows and block size $B$ define $\textbf{bsp}(A)^{M_{1}\times M_{2}}_{B\times B}$ (block sparsity pattern) as a function

[TABLE]

where $\mathbb{B}=\left\{0,1\right\}$. The function takes block matrix $A\in\mathbb{R}^{M_{1}B\times M_{2}B}$ as input and returns as output the matrix $R=\textbf{bsp}(A)^{M_{1}\times M_{2}}_{B\times B}\in\mathbb{B}^{M_{1}\times M_{2}}$ such that

[TABLE]

By $\#\textbf{bsp}(A)^{M_{1}\times M_{2}}_{B\times B}$ define the number of nonzero blocks of matrix $A$ and the number of ones in matrix $R$.

Proposition 2**.**

If the matrix $A$ has $\mathcal{H}^{2}$ structure, the compression Algorithm 1 has $L$ levels, the block size on each level is $B$, the matrix $A$ has zero level close matrix $C$, and the far blocks are compressed with rank $r=B/2$, then the compression algorithm for the matrix $A$ leads to the factorization:

[TABLE]

where $U$ and $V$ are orthogonal matrices equal to the multiplication of block-diagonal compression and permutation matrices

[TABLE]

$S$* is a sparse matrix that has*

[TABLE]

nonzero blocks222The symbol # before the matrix means the number of nonzero $(r\times r)$ blocks in this matrix. of size $(r\times r)$.

Proof.

Consider the matrix $S$ from (11). Let matrix $S_{ij}$ correspond to $i$-th level non-basis hyper row and $j$-th level non-basis hyper column. Thanks to basis-non-basis row and column permutations $P_{ri}$ and $P_{ci}$ matrix $S$ is separated into blocks $S_{ij}$, where $i,j\in{0,\dots,L}$.

The number of nonzero blocks in $S$ is equal to sum of nonzero blocks in $S_{ij}$:

[TABLE]

Let us compute the number of nonzero blocks in block $S_{ij}$. Since on each level we join block rows by groups of $J$ blocks, we obtain:

[TABLE]

where $i,j\in 0,\dots,L$.

Note that

[TABLE]

Thus

[TABLE]

[TABLE]

Obtain

[TABLE]

Thus, the number of nonzero blocks in matrix $S$ is less than the number of nonzero blocks in close matrix $C$ multiplied by constant ($4L+6(\frac{1}{2^{L}}-1)$), if matrix $C$ is block-sparse, and all proposition conditions are met, then matrix $S$ is also sparse. ∎

4 Building sparse factorization from $\mathcal{H}^{2}$ coefficients

4.1 Definition of $\mathcal{H}^{2}$ matrix

In this section we consider the matrix $A$ approximated in $\mathcal{H}^{2}$ format. There exists a number of efficient ways to build this approximation[22, 6]. In this paper, we do not consider the process of building the $\mathcal{H}^{2}$ matrix and assume that it is given. Let us explain in details how to construct factors in the decomposition (11) from parameters of $\mathcal{H}^{2}$ matrix $A$. First, we give the definition of $\mathcal{H}^{2}$ matrix, the more detailed definition can be found in[6].

Definition 3** (Row and column cluster trees).**

Cluster trees of rows and columns $\mathcal{T}_{r}$ and $\mathcal{T}_{c}$ define the hierarchical division of block rows and columns. At each level of the row cluster tree $\mathcal{T}_{r}$, each node corresponds to a block row of the matrix $A$, child nodes correspond to the subrows of this row. Same for the column cluster tree $\mathcal{T}_{c}$.

Definition 4** (Block cluster tree).**

Let $\mathcal{T}_{rc}$ be a tree. $\mathcal{T}_{rc}$ is a block cluster tree for $\mathcal{T}_{r}$ and $\mathcal{T}_{c}$ if it satisfies the following conditions:

• •

$\text{root}(\mathcal{T}_{rc}$) = ($\text{root}(\mathcal{T}_{r}$),$\text{root}(\mathcal{T}_{c}$)).

• •

Each node $b\in\mathcal{T}_{rc}$ has the form $b=(t,s)$ for $t\in\mathcal{T}_{r}$ and $s\in\mathcal{T}_{c}$.

• •

Let $b=(t,s)\in\mathcal{T}_{rc}$. If $\text{sons}(b)\neq\varnothing$, then

[TABLE]

Definition 5** (Admissibility condition).**

Let $\mathcal{T}_{r}$ and $\mathcal{T}_{r}$ be a row and column cluster trees. A predicate

[TABLE]

is an admissibility condition for $\mathcal{T}_{r}$ and $\mathcal{T}_{r}$ if

[TABLE]

and

[TABLE]

If $\mathcal{A}(t,s)$, the pair (t,s) is called admissible.

Definition 6** (Admissibility block cluster tree).**

Let $\mathcal{T}_{rc}$ be a block cluster tree for $\mathcal{T}_{r}$ and $\mathcal{T}_{c}$, let $\mathcal{A}$ be an admissibility condition. If for each $(t,s)\in\mathcal{T}_{rc}$ either $\text{sons}(t)\neq\varnothing\neq\text{sons}(s)$ or $\mathcal{A}(t,s)=$ True holds, the block cluster tree $\mathcal{T}_{rc}$ called $\mathcal{A}$-Admissible.

Definition 7** (Farfield and nearfield).**

Let $\mathcal{T}_{rc}$ be a block cluster tree for $\mathcal{T}_{r}$ and $\mathcal{T}_{c}$, let $\mathcal{A}$ be an admissibility condition. The index set

[TABLE]

is called set of farfield blocks. The index set

[TABLE]

is called set of nearfield blocks.

Definition 8** (Cut-off matrices).**

Let $\mathcal{T}_{r}$ be a row cluster tree. For all tree nodes $t\in\mathcal{T}_{r}$ the cut-off matrix $\chi_{t}\in\mathbb{R}^{N\times N}$ corresponding to $t$ is defined by

[TABLE]

Definition 9** (Cluster basis).**

Let $K=(K_{t})_{t\in\mathcal{T}_{r}}$ be a family of finite index sets (rank distribution for $\mathcal{T}_{r}$). Let $R=(R_{t})_{t\in\mathcal{T}_{r}}$ be a family of matrices satisfying $R_{t}\in\mathbb{R}_{\hat{t}}^{N\times K_{t}}$ for all $t\in\mathcal{T}_{r}$. Then $R$ is called row cluster basis and the matrices $R_{t}$ are called row cluster basis matrices. Analogically for column cluster basis $E$ and column cluster basis matrices $E_{s}$, $s\in\mathcal{T}_{c}$.

Definition 10** (Close matrix).**

The matrix $C\in\mathbb{R}^{N\times N}$ is close if

[TABLE]

Definition 11** ($\mathcal{H}^{2}$ matrix).**

The matrix $A\in\mathbb{R}^{N\times N}$ is approximated in $\mathcal{H}^{2}$ format if $\mathcal{T}_{c}$ and $\mathcal{T}_{r}$ are block cluster trees of columns and rows of matrix $A$, if there exist a row cluster basis $R$ (row transition matrices), the column cluster basis $E$ (column transition matrices), a family $D=(D_{b})_{b\in\mathcal{I}_{N\times N}^{+}}$ of matrices satisfying $D_{b}\in\mathbb{R}^{K_{t}\times L_{s}}$ for all $b=(s,t)\in\mathcal{I}_{N\times N}^{+}$ (interaction list), and close matrix $C$, and if

[TABLE]

4.2 Construction of matrices $U$ and $V$ from the coefficients of $\mathcal{H}^{2}$ matrix

Let us first construct orthogonal matrices $U$ and $V$ from the factorization (11). According to equation (10):

[TABLE]

and

[TABLE]

Note that matrices $U_{k}$ and $V_{k},k\in 0,\dots,L$, are very close in their meaning to cluster basis matrices $R_{t}$ and $E_{s}$ both are level compression matrices. The difference between these matrices is that the diagonal blocks of matrices $U_{k}$ and $V_{k}$ are square orthogonal blocks, and diagonal blocks of matrices $R_{t}$ and $E_{s}$ are rectangular non-orthogonal blocks.

Thus we can take matrices $R_{t}$ and $E_{s}$, orthogonalize blocks, complete each block to square orthogonal block and obtain the matrices $U_{k}$ and $V_{k}$. The algorithm that orthogonalizes blocks of matrices $R_{t}$ and $E_{s}$ is known as the $\mathcal{H}^{2}$ compression algorithm, it can be found in[6, 7]. Compression of blocks can be done by QR decomposition of blocks with the square $Q$ factor. Permutations $P_{rk}$ and $P_{ck}$ can be constructed from the cluster trees.

4.3 Construction of the matrix $S$ from the coefficients of $\mathcal{H}^{2}$ matrix

According to equations (9), equation (5) and equation (6):

[TABLE]

[TABLE]

Construction of matrices $U_{i}$ and $V_{i}$ is shown in the previous subsection, matrix $C_{0}$ is stored in the $\mathcal{H}^{2}$ matrix explicitly as close matrix $C$, matrices $\widehat{F}_{\mathrm{ml}i}$ are exactly matrices $D_{i}$ from interaction list, matrix $\widehat{F}_{L}$ is matrix $D_{L}$. Thus, matrix $S$ can be easily computed from the coefficients of the $\mathcal{H}^{2}$ matrix.

5 Numerical experiments

Sparsification algorithm is implemented in the Python programming language. For the $\mathcal{H}^{2}$ matrix implementation we use the h2tools[22] library. All computations are performed on MacBook Air with a 1.3GHz Intel Core i5 processor and 4 GB 1600 MHz DDR3 RAM.

First, we numerically show that the matrix $S$ in factorization (11) is indeed sparse. Then we give the timing and storage requirements of the sparse factorization. Thereafter we consider the sparse factorization of the $\mathcal{H}^{2}$ matrix combined with a sparse direct solver as a direct solver for the system with the $\mathcal{H}^{2}$ matrix and compare this approach with HODLR direct solver and $\mathcal{H}^{2}$-LU solver from $\mathcal{H}$2Lib library. Finally, we consider sparse factorization with the matrix $S$ factorized by the ILUt method as a preconditioner for GMRES solver.

We want to note that the sparsification algorithm does not worsen the accuracy of the $\mathcal{H}^{2}$ approximation. It follows from the algorithm and it is confirmed in the experiments. Therefore, in the experiments below, we omit the accuracy of the sparse factorization and show only the $\mathcal{H}^{2}$ approximation accuracy to avoid redundancy.

5.1 Sparsity of the factor $S$

In Section 3 we have studied the sparsity of the factor $S$ from the factorization (11) analytically. Here we present the numerical illustration that matrix $S$ is indeed sparse. Tests are performed for two $\mathcal{H}^{2}$ matrices.

Example 5.1**.**

$\mathcal{H}^{2}$* approximation of the matrix*

[TABLE]

where $r_{i}\in\mathbb{R}^{2}$ or $r_{i}\in\mathbb{R}^{3}$ is the position of the $i$-th element. Elements are randomly distributed in identity square $\Omega=[0,1]\times[0,1]$ in $\mathbb{R}^{2}$ case and in identity cube $\Omega=[0,1]\times[0,1]\times[0,1]$ in $\mathbb{R}^{3}$ case. The $\mathcal{H}^{2}$ approximation accuracy is $\epsilon=10^{-6}$. In Figure 5 the factor $S$ of the sparse factorization of the $\mathcal{H}^{2}$ matrix with the core (12) is marked by “inv”.

Example 5.2**.**

$\mathcal{H}^{2}$* approximation of the matrix*

[TABLE]

where $r_{i}\in\mathbb{R}^{2}$ or $r_{i}\in\mathbb{R}^{3}$ is the position of the $i$-th element. Elements are randomly distributed in identity square $\Omega=[0,1]\times[0,1]$ in $\mathbb{R}^{2}$ case and in identity cube $\Omega=[0,1]\times[0,1]\times[0,1]$ in $\mathbb{R}^{3}$ case. The $\mathcal{H}^{2}$ approximation accuracy is $\epsilon=10^{-6}$, In Figure 5 the factor $S$ of the sparse factorization of the $\mathcal{H}^{2}$ matrix with the core (13) is marked by “exp”.

We consider the sparse factorization (11) of the $\mathcal{H}^{2}$ matrices from Examples 5.1 and 5.2 in 2D and 3D. In Figure 5 we show the number of nonzero elements per row in the factor $S$. Since for all considered matrices the number of nonzero elements per row is constant (and does not grow with the matrix size), we conclude that the matrix $S$ in factorization (11) is indeed sparse.

5.2 The storage requirements and timing of the sparsification algorithm

In Figure 6 we show the time of the sparse factorization of $\mathcal{H}^{2}$ matrices with the cores (12) and (13) (approximation accuracy is $\varepsilon=10^{-6}$) denoted by ”inv” and ”exp”.

Sparsification time grows almost linearly. In Figure 7 we show the memory requirements of the matrix (12) in 2D and 3D, also we show the memory requirements of its $\mathcal{H}^{2}$ approximation ($\varepsilon=10^{-6}$), of the sparse factorization (sum of $U,$ $S,$ and $V$) and of the factor $S$ separately.

In Figure 8 we show the memory requirements of the matrix (13) in 2D and 3D, also we show the memory requirements of its $\mathcal{H}^{2}$ approximation ($\varepsilon=10^{-6}$), of the sparse factorization (sum of $U,$ $S,$ and $V$) and of the factor $S$ separately.

For all examples, the memory requirements of both $\mathcal{H}^{2}$ and the sparse factorization grows almost linearly, unlike the memory requirements of the original matrix which scales quadratically.

5.3 Comparison to HODLR

Paper[3] considers the HODLR approximation of the dense matrix and its factorization as an efficient way to compute the determinant of the dense matrix. We propose the $\mathcal{H}^{2}$ approximation of the dense matrix, its sparsification and factorization of the sparse matrix as an alternative. The triangular factorization of the sparse matrix is computed by CHOLMOD[12] package. Tests are performed for 3D data, for matrix

[TABLE]

where $r_{i}\in\mathbb{R}^{3}$ is the position of the $i$-th element. Both HODLR and $\mathcal{H}^{2}$ approximation accuracy is $\varepsilon=10^{-6}$. The HODLR factorization accuracy is $\vartheta=10^{-5}$, the accuracy of the triangular factorization of the sparse matrix is $\vartheta=10^{-10}$ (which is redundant, but the used package has no accuracy options). In Figure 9 we show the time comparison for this two approaches.

Total solution time comparison in 2D and 3D is presented in Figure 10.

5.4 Comparison to H2Lib liberary

In this subsection we compare the approach proposed in this paper (sparse non-extensive factorization of the $\mathcal{H}^{2}$ matrix, triangular factorization of the sparse matrix and then the solution of the system) with the approach based on $\mathcal{H}^{2}$-LU factorization, proposed in the work[9] and implemented in H2lib package[5]. The $\mathcal{H}^{2}$-LU factorization takes the $\mathcal{H}^{2}$ matrix and returns the triangular factors $L$ and $U$ in $\mathcal{H}^{2}$ format. The triangular factorization of the sparse matrix is computed by CHOLMOD[12] package. Tests are performed on the following problem.

Example 5.3**.**

Consider the Dirichlet boundary value problem for Laplace’s equation

[TABLE]

where $\Omega=[-1,1]^{3}$ is a cube. The standard technic: using the single layer potential we obtain boundary integral formulation of the equation (14).

[TABLE]

where $G(x)=\frac{1}{4\pi}\frac{1}{|x|}$ is the fundamental solution for the Laplace operator. Then we discretize the integral equation (15) on the triangular grid on $\Gamma$ using the Galerkin method. Obtained dense matrix is approximated in $\mathcal{H}^{2}$ format with accuracy $\varepsilon=10^{-6}$.

The accuracy of the $\mathcal{H}^{2}$-LU factorization is $\vartheta=10^{-6}$, the accuracy of the triangular factorization of the sparse matrix is $\vartheta=10^{-10}$.

In Table 1 we show the time comparison of the solution of the system with matrix from Example 5.3, using $\mathcal{H}^{2}$-LU and sparsification approaches. For the $\mathcal{H}^{2}$-LU we show the approximation in $\mathcal{H}^{2}$ format and factorization time, for the sparsification we show the approximation in $\mathcal{H}^{2}$ format, sparsification and sparse factorization time. In both cases, time of the solution of the system with factorized matrix is negligible, so we do not show it.

In Figure 11 we show the comparison of the total time, required for the system solution.

The sparsification approach has not only better timing, but also better asymptotics.

5.5 Sparsification method as a preconditioner

$\mathcal{H}^{2}$ matrix is an efficient tool to multiply a matrix by a vector. This allows to apply iterative solvers like GMRES to solution of the systems with $\mathcal{H}^{2}$ matrix. But the preconditioning is still a challenging problem due to complexity of the factorization of the $\mathcal{H}^{2}$ matrix. We propose to use the approximately factored sparsification of the $\mathcal{H}^{2}$ matrix as a preconditioner to iterative method. For tests we use randomly distributed 3D data with following interaction matrix.

Example 5.4**.**

[TABLE]

where $r_{i}\in\mathbb{R}^{3}$ is the position of the $i$-th element.

This example is used for testing of IFMM (Inverse Fast Multipole Method) method as a preconditioning in[11], so we have chosen this example for the convenient comparison.

This matrix is useful for the iterative tests since condition number of this matrix significantly depends on the parameter $d$: the larger $d$ is, the larger condition number is.

Example with well-conditioned matrix. First, consider the matrix from Example 5.4 with $d=10^{-3}$, condition number $\mathrm{cond}(A)=10$. We solve this system using GMRES iterative solver for the matrix approximated in $\mathcal{H}^{2}$ format with accuracy $\epsilon=10^{-9}$ and as a preconditioner we use $\mathcal{H}^{2}$ approximation of the matrix $A$ with accuracy $\epsilon=10^{-3}$ sparsified and factorized with ILUt decomposition. Figure 12 shows convergence of the GMRES method with different ILUt threshold parameters. The required residual of the GMRES method $r=10^{-10}$.

The standard trade off: the more time on the preconditioner building we spend, the faster iterations converge. Figure 13 illustrates the total time required for the system solution (including the sparsification construction).

We show the total time required for $\mathcal{H}^{2}$ matrix sparsification, time required for building the ILUt preconditioner with dropping tolerance $\tau=2\times 10^{-2}$ and iterations timing in Figure 14(a).

Example with ill-conditioned matrix. Consider the matrix from Example 5.4 with $d=10^{-2}$, condition number $\mathrm{cond}(A)=10^{4}$. As in the previous paragraph, we solve this system using GMRES iterative solver for the matrix approximated in $\mathcal{H}^{2}$ format with accuracy $\epsilon=10^{-9}$ and as a preconditioner we used $\mathcal{H}^{2}$ approximation of the matrix $A$ with accuracy $\epsilon=10^{-3}$ sparsified and factorized with ILUt decomposition. In Figure 15 convergence till the tolerance $10^{-10}$ for different ILUt parameters is shown.

In Figure 16 the total solution time for different ILUt parameters is shown.

As we can see, the optimal ILUt parameter for this problem is $\tau=10^{-2}$. Figure 17(b) presents the total time required for solution of the system with optimal ILUt parameter.

We compared the results of our solver (GMRES preconditioned by sparsification method) with results of IFMM solver presented in[11] for the same problem, same parameters and similar hardware. For test we used matrix from Example 5.4 with 3D data, matrix size $N=10^{5}$, iterative method GMRES till the residual $r=10^{-10}$, parameters $d=\{10^{-3},10^{-2}\}$ for well- and ill-conditioned systems. We present the best total solution time achieved in experiments for both methods in Table 2.

Note that sparsification method is implemented in Python programing language and can be significantly improved by applying Cython or switching to another programing language (as C++ or Fortran).

6 Related work

Hierarchical low-rank matrix formats such as $\mathcal{H}$[15, 17, 16, 8], HODLR[1, 3] (Hierarchical Off-Diagonal Low-Rank), HSS[21, 10, 23] (Hierarchically

Semiseparable), $\mathcal{H}^{2}$ [16, 6] matrices and etc., that are matrix analogies of the fast multipole method[14, 13], have two significant features: they do store information in data-sparse formats and they provide the fast matrix by vector product. Fast ($\mathcal{O}(N)$, where $N$ is size of the matrix) matrix by vector product allows to apply iterative solvers. Data-sparse representation allows to store matrix in $\mathcal{O}(N)$ cells of memory, but storage scheme is usually complicated.

If the hierarchical matrix is ill-conditioned, then pure iterative solver fails and it is required to apply either approximate direct solver or preconditioner (that is also approximate direct solver, probably with lower accuracy). Due to complex storage schemes of hierarchical matrices, construction of the approximate direct solver is a challenging problem. There exists two general approaches to approximate direct solution of $\mathcal{H}^{2}$ matrix: factorization of hierarchical matrix[4, 21, 23], and sparsification of the hierarchical matrix followed by factorization of the sparse matrix[2, 26].

The factorization approach is more popular for hierarchical matrices with strong low-rank structure, also known as hierarchical matrices with weak-admissibility criteria[18]

($\mathcal{H}$[15, 17], HODLR[1, 3], HSS[21, 10, 23] matrices). For the $\mathcal{H}$ matrix, the algorithm $\mathcal{H}$-LU[4] with almost linear complexity was proposed. This algorithm has been successfully applied to many problems. The major drawback of the $\mathcal{H}$-LU algorithm is that factorization time and memory required for $L$ and $U$ factors can be quite large. Approximate direct solvers based on factorization of HSS and HODLR matrices are also well studied and found many[23, 28, 20, 24, 1, 3, 21] successful applications.

One of the approaches to the solution of the systems with $\mathcal{H}^{2}$ matrices is the sparse factorization (sparsification). The sparsification approach is usually applied to hierarchical matrices with weak-admissibility criteria ($\mathcal{H}^{2}$ matrices). Sparsification algorithms transform $\mathcal{H}^{2}$ matrix into the sparse matrix and then factorize the sparse matrix. Algorithm, proposed in this paper is the sparsification algorithm. The main difference between presented work and the other sparse factorizations[2, 11, 26] is a size of sparse factors. The main benefit of the non-extensive sparsification is that it preserves the size of the factorized $\mathcal{H}^{2}$ matrix, while the other sparsification algorithms return extended factors.

7 Conclusions

We have proposed a new approach to the solution of the systems with $\mathcal{H}^{2}$ matrices that is based on sufficient non-extensive sparsification of the $\mathcal{H}^{2}$ matrix. Proposed sparsification is suitable for any $\mathcal{H}^{2}$ matrices (including non-symmetric) and preserves such important properties of the matrix as its size, symmetry (if exists) and positive definite (if exists).