Accuracy Controlled Structure-Preserving ${\cal H}^2$-Matrix-Matrix Product in Linear Complexity with Change of Cluster Bases

TL;DR

This paper introduces an accuracy-controlled ${ m H}^2$-matrix-matrix product algorithm that dynamically adjusts cluster bases during computation, maintaining linear complexity and preserving the matrix structure, suitable for large-scale electromagnetic problems.

Contribution

It presents a novel ${ m H}^2$-matrix-matrix product algorithm with explicit accuracy control via dynamic cluster basis changes, ensuring linear complexity and structure preservation.

Findings

Achieves accurate ${ m H}^2$-matrix products with linear complexity

Demonstrates efficiency in large-scale electromagnetic simulations

Maintains matrix structure during computation

Abstract

${\cal H}^2$-matrix constitutes a general mathematical framework for efficient computation of both partial-differential-equation and integral-equation-based operators. Existing linear-complexity ${\cal H}^2$ matrix-matrix product (MMP) algorithm lacks explicit accuracy control, while controlling accuracy without compromising linear complexity is challenging. In this paper, we develop an accuracy controlled ${\cal H}^2$ matrix-matrix product algorithm by instantaneously changing the cluster bases during the matrix product computation based on prescribed accuracy. Meanwhile, we retain the computational complexity of the overall algorithm to be linear. Different from the existing ${\cal H}^2$ matrix-matrix product algorithm where formatted multiplications are performed using the original cluster bases, in the proposed algorithm, all additions and multiplications are either exact or computed based on prescribed accuracy. Furthermore, the original ${\cal H}^2$-matrix structure is preserved in the matrix product. While achieving optimal complexity for constant-rank matrices, the computational complexity of the proposed algorithm is also minimized for variable-rank ${\cal H}^2$-matrices. The proposed work serves as a fundamental arithmetic in the development of fast solvers for large-scale electromagnetic analysis. Applications to both large-scale capacitance extraction and electromagnetic scattering problems involving millions of unknowns on a single core have demonstrated the accuracy and efficiency of the proposed algorithm.

Tables

Table 1. TABLE I: ℋ 2 superscript ℋ 2 {\cal H}^{2} MMP error at different ϵ t ​ r ​ u ​ n ​ c subscript italic-ϵ 𝑡 𝑟 𝑢 𝑛 𝑐 \epsilon_{trunc} for large-scale capacitance extraction matrices as a function of N 𝑁 N .

$N$	4,480	17,152	67,072	265,216	1,054,720
${ϵ}_{tr}:$ 1E-2	4.35E-2	5.72E-2	5.73E-2	5.80E-2	5.97E-2
${ϵ}_{tr}:$ 1E-4	3.71E-3	3.72E-3	3.86E-3	3.80E-3	3.67E-3
${ϵ}_{tr}:$ 1E-6	2.82E-4	3.28E-4	3.86E-4	4.50E-4	5.66E-4

Table 2. TABLE II: ℋ 2 superscript ℋ 2 {\cal H}^{2} MMP error for 2-D slab. for different ϵ t ​ r ​ u ​ n ​ c subscript italic-ϵ 𝑡 𝑟 𝑢 𝑛 𝑐 \epsilon_{trunc} as a function of N 𝑁 N .

$N$	22560	89920	359040	1098720
1E-2	8.54E-3	1.06E-2	1.49E-2	1.07E-2
1E-3	2.52E-3	3.17E-3	4.23E-3	3.79E-3
1E-4	7.86E-4	9.76E-4	1.38E-3	1.23E-3
1E-5	2.91E-4	3.37E-4	4.22E-4	4.11E-4
1E-6	8.04E-5	9.85E-5	1.27E-4	1.36E-4

Table 3. TABLE III: C s ​ p subscript 𝐶 𝑠 𝑝 C_{sp} as a function of N 𝑁 N for the dieletric cube array

$N$	3024	24192	193536	1548288
$C_{sp}$	16	42	95	126

Table 4. TABLE IV: ℋ 2 superscript ℋ 2 {\cal H}^{2} MMP error at different ϵ t ​ r ​ u ​ n ​ c subscript italic-ϵ 𝑡 𝑟 𝑢 𝑛 𝑐 \epsilon_{trunc} for 3-D cube array.

$N$	3024	24192	193536	1548288
Existing [1]	9.02E-2	1.01E-1	1.77E-1	2.74E-1
1E-2	1.91E-2	2.38E-2	3.82E-2	6.58E-2
1E-3	5.51E-3	7.23E-3	1.06E-2	2.16E-2
1E-4	1.48E-3	2.46E-3	3.69E-3	8.09E-3

Equations

$$max\{diam({\Omega }_t),diam({\Omega }_s)\}\le \eta dist({\Omega }_t,{\Omega }_s),$$

$${\text{R}}_{t,s}=({\text{V}}_t{)}_{\#t\times k}({\text{S}}_{t,s}{)}_{k\times k}({\text{V}}_s{)}_{\#s\times k}^T$$

$$({\text{V}}_t{)}_{\#t\times k}=\left[\begin{array}{cc}({\text{V}}_{t_1}{)}_{\#t_1\times k_1}& 0\\ 0& ({\text{V}}_{t_2}{)}_{\#t_2\times k_2}\end{array}\right]\left[\begin{array}{c}({\text{T}}_{t_1}{)}_{k_1\times k}\\ ({\text{T}}_{t_2}{)}_{k_2\times k}\end{array}\right]$$

$$\text{case-4: }{\text{R}}_{i,j}^{\text{A}}\times {\text{R}}_{j,k}^{\text{B}}={\text{V}}_{i_r}^{\text{A}}{\text{S}}_{i,j}^{\text{A}}({\text{V}}_{j_c}^{\text{A}}{)}^T\times {\text{V}}_{j_r}^{\text{B}}{\text{S}}_{j,k}^{\text{B}}({\text{V}}_{k_c}^{\text{B}}{)}^T,$$

$${\text{R}}_{i,j}^{\text{A}}\times {\text{R}}_{j,k}^{\text{B}}={\text{V}}_{i_r}^{\text{A}}{\text{S}}_{i,k}^{\text{C}}({\text{V}}_{k_c}^{\text{B}}{)}^T,$$

$$\text{case-1: }{\text{F}}_{i,j}^{\text{A}}\times {\text{F}}_{j,k}^{\text{B}};$$

$$\text{case-2: }{\text{F}}_{i,j}^{\text{A}}\times {\text{R}}_{j,k}^{\text{B}}=({\text{F}}_{i,j}^{\text{A}}{\text{V}}_{j_r}^{\text{B}})\times {\text{S}}_{j,k}^{\text{B}}\times ({\text{V}}_{k_c}^{\text{B}}{)}^T;$$

$$\text{case-3: }{\text{R}}_{i,j}^{\text{A}}\times {\text{F}}_{j,k}^{\text{B}}={\text{V}}_{i_r}^{\text{A}}\times {\text{S}}_{i,j}^{\text{A}}\times \left(({\text{V}}_{j_c}^{\text{A}}{)}^T{\text{F}}_{j,k}^{\text{B}}\right).$$

$${\text{G}}_{i_{r1}}^{\text{C},L}=\sum _{j=1}^{O(C_{sp})}\sum _{k=1}^{O(C_{sp})}({\text{F}}_{i,j}^{\text{A}}{\text{F}}_{j,k}^{\text{B}})({\text{F}}_{i,j}^{\text{A}}{\text{F}}_{j,k}^{\text{B}}{)}^H,$$

$${\text{G}}_{i_{r2}}^{\text{C},L}=\sum _{j=1}^{O(C_{sp})}({\text{F}}_{i,j}^{\text{A}}{\text{V}}_{j_r}^{\text{B}})({\text{F}}_{i,j}^{\text{A}}{\text{V}}_{j_r}^{\text{B}}{)}^H,$$

$${\text{G}}_{i_{r3}}^{\text{C},L}={\text{V}}_{i_r}^{\text{A}}({\text{V}}_{i_r}^{\text{A}}{)}^H.$$

$${\text{G}}_{i_{r3}}^{\text{C},L}=\widehat{{\text{G}}_{i_{r1}}^{\text{C},L}}+\widehat{{\text{G}}_{i_{r2}}^{\text{C},L}}+\widehat{{\text{G}}_{i_{r3}}^{\text{C},L}}.$$

$${\text{G}}_{k_{c1}}^{\text{C},L}=\sum _{i=1}^{O(C_{sp})}\sum _{j=1}^{O(C_{sp})}({\text{F}}_{i,j}^{\text{A}}{\text{F}}_{j,k}^{\text{B}}{)}^T({\text{F}}_{i,j}^{\text{A}}{\text{F}}_{j,k}^{\text{B}}{)}^{\ast }.$$

$${\text{G}}_{k_{c2}}^{\text{C},L}=\sum _{j=1}^{O(C_{sp})}{\left(({\text{V}}_{j_c}^{\text{A}}{)}^T{\text{F}}_{j,k}^{\text{B}}\right)}^T{\left(({\text{V}}_{j_c}^{\text{A}}{)}^T{\text{F}}_{j,k}^{\text{B}}\right)}^{\ast }.$$

$${\text{G}}_{k_{c3}}^{\text{C},L}={\text{V}}_{k_c}^{\text{B}}({\text{V}}_{k_c}^{\text{B}}{)}^H.$$

$${\text{G}}_{k_c}^{\text{C},L}=\widehat{{\text{G}}_{k_{c1}}^{\text{C},L}}+\widehat{{\text{G}}_{k_{c2}}^{\text{C},L}}+\widehat{{\text{G}}_{k_{c3}}^{\text{C},L}}.$$

$${\text{S}}_{i,k}^{\text{C}}=\{\begin{array}{ll}({\text{V}}_{i_r}^{\text{C}}{)}^H{\text{F}}_{i,j}^{\text{A}}{\text{F}}_{j,k}^{\text{B}}({\text{V}}_{k_c}^{\text{C}}{)}^{\ast }& \text{case-1}\\ ({\text{V}}_{i_r}^{\text{C}}{)}^H{\text{F}}_{i,j}^{\text{A}}{\text{V}}_{j_r}^{\text{B}}{\text{S}}_{j,k}^{\text{B}}({\text{V}}_{k_c}^{\text{B}}{)}^T({\text{V}}_{k_c}^{\text{C}}{)}^{\ast }& \text{case-2}\\ ({\text{V}}_{i_r}^{\text{C}}{)}^H{\text{V}}_{i_r}^{\text{A}}{\text{S}}_{i,j}^{\text{A}}({\text{V}}_{j_c}^{\text{A}}{)}^T{\text{F}}_{j,k}^{\text{B}}({\text{V}}_{k_c}^{\text{C}}{)}^{\ast }& \text{case-3}\\ ({\text{V}}_{i_r}^{\text{C}}{)}^H{\text{V}}_{i_r}^{\text{A}}{\text{S}}_{i,j}^{\text{A}}{\text{B}}_j{\text{S}}_{j,k}^{\text{B}}({\text{V}}_{k_c}^{\text{B}}{)}^T({\text{V}}_{k_c}^{\text{C}}{)}^{\ast }& \text{case-4}.\end{array}$$

$${\text{B}}_j=({\text{V}}_{j_c}^{\text{A}}{)}^T\times {\text{V}}_{j_r}^{\text{B}}.$$

$$\begin{array}[]{l}\textbf{P}_{i}^{\textbf{A}}=(\textbf{V}_{i_{r}}^{\textbf{C}})^{H}\textbf{V}_{i_{r}}^{\textbf{A}};\\ \textbf{P}_{k}^{\textbf{B}}=(\textbf{V}_{k_{c}}^{\textbf{B}})^{T}(\textbf{V}_{k_{c}}^{\textbf{C}})^{*}\end{array}$$

$$\begin{array}[]{l}(\textbf{F}_{i,j}^{\textbf{A}}\textbf{F}_{j,k}^{\textbf{B}})_{coll.}=(\textbf{V}_{i_{r}}^{\textbf{C}})^{H}(\textbf{F}_{i,j}^{\textbf{A}}\textbf{F}_{j,k}^{\textbf{B}})\left(\textbf{V}_{k_{c}}^{\textbf{C}}\right)^{*}\\ (\textbf{F}_{i,j}^{\textbf{A}})_{coll.}=(\textbf{V}_{i_{r}}^{\textbf{C}})^{H}\textbf{F}_{i,j}^{\textbf{A}}\textbf{V}_{j_{r}}^{\textbf{B}}\\ (\textbf{F}_{j,k}^{\textbf{B}})_{coll.}=(\textbf{V}_{j_{c}}^{\textbf{A}})^{T}\textbf{F}_{j,k}^{\textbf{B}}\left(\textbf{V}_{k_{c}}^{\textbf{C}}\right)^{*},\end{array}$$

$${\text{S}}_{i,k}^{\text{C}}=\{\begin{array}{ll}({\text{V}}_{i_r}^{\text{C}}{)}^H{\text{F}}_{i,j}^{\text{A}}{\text{F}}_{j,k}^{\text{B}}({\text{V}}_{k_c}^{\text{C}}{)}^{\ast }& \text{case-1}\\ ({\text{F}}_{i,j}^{\text{A}}{)}_{coll.}{\text{S}}_{j,k}^{\text{B}}{\text{P}}_k^{\text{B}}& \text{case-2}\\ {\text{P}}_i^{\text{A}}{\text{S}}_{i,j}^{\text{A}}({\text{F}}_{j,k}^{\text{B}}{)}_{coll.}& \text{case-3}\\ {\text{P}}_i^{\text{A}}{\text{S}}_{i,j}^{\text{A}}{\text{B}}_j{\text{S}}_{j,k}^{\text{B}}{\text{P}}_k^{\text{B}}& \text{case-4}.\end{array}$$

$${\text{R}}_{i,j}^{\text{A}}\times {\text{R}}_{j,k}^{\text{B}}=$$

$${\text{V}}_{i_r^{ch}}^{\text{A}}{\text{T}}_{i_r}^{\text{A}}{\text{S}}_{i,j}^{\text{A}}({\text{T}}_{j_c}^{\text{A}}{)}^T({\text{V}}_{j_c^{ch}}^{\text{A}}{)}^T\times {\text{V}}_{j_r^{ch}}^{\text{B}}{\text{T}}_{j_r}^{\text{B}}{\text{S}}_{j,k}^{\text{B}}({\text{T}}_{k_c}^{\text{B}}{)}^T({\text{V}}_{k_c^{ch}}^{\text{B}}{)}^T,$$

$${\text{S}}_{i,k}^{\text{C}}={\text{S}}_{i,j}^{\text{A}}({\text{B}}_j){\text{S}}_{j,k}^{\text{B}}.$$

$${\text{R}}_{i,j}^{\text{A}}\times {\text{R}}_{j,k}^{\text{B}}={\text{V}}_{i_r^{ch}}^{\text{A}}{\text{T}}_{i_r}^{\text{A}}{\text{S}}_{i,k}^{\text{C}}({\text{T}}_{k_c}^{\text{B}}{)}^T({\text{V}}_{k_c^{ch}}^{\text{B}}{)}^T.$$

$$\text{case-1: }{\text{NL}}_{i,j}^{\text{A}}\times {\text{NL}}_{j,k}^{\text{B}};$$

$$\text{case-2: }{\text{NL}}_{i,j}^{\text{A}}\times {\text{R}}_{j,k}^{\text{B}}=({\text{NL}}_{i,j}^{\text{A}}{\text{V}}_{j_r}^{\text{B}})\times {\text{S}}_{j,k}^{\text{B}}\times ({\text{V}}_{k_c}^{\text{B}}{)}^T;$$

$$\text{case-3: }{\text{R}}_{i,j}^{\text{A}}\times {\text{NL}}_{j,k}^{\text{B}}={\text{V}}_{i_r}^{\text{A}}\times {\text{S}}_{i,j}^{\text{A}}\times \left(({\text{V}}_{j_c}^{\text{A}}{)}^T{\text{NL}}_{j,k}^{\text{B}}\right).$$

$${\text{NL}}_{i,j}^{\text{A}}\times {\text{NL}}_{j,k}^{\text{B}}=$$

$$\left[\begin{array}{cc}{\text{V}}_{i{1}_r}^{\text{C}}& \\ & {\text{V}}_{i{2}_r}^{\text{C}}\end{array}\right]({\text{NL}}_{i,j}^{\text{A}}{\text{NL}}_{j,k}^{\text{B}}{)}_{coll.}\left[\begin{array}{cc}({\text{V}}_{k{1}_c}^{\text{C}}{)}^T& \\ & ({\text{V}}_{k{2}_c}^{\text{C}}{)}^T\end{array}\right],$$

$$({\text{NL}}_{i,j}^{\text{A}}{\text{NL}}_{j,k}^{\text{B}}{)}_{coll.}=$$

Full text

I Introduction

The ${\cal H}^{2}$-matrix [1, 2] constitutes a general mathematical framework for compact representation and efficient computation of large dense systems. Both partial differential equation (PDE) and integral equation (IE) operators in electromagnetics can be represented as ${\cal H}^{2}$-matrices with controlled accuracy [3, 4, 5].

The development of ${\cal H}^{2}$-matrix arithmetic such as addition, multiplication, and inverse are of critical importance to the development of fast solvers in electromagnetics [7]. Under the ${\cal H}^{2}$-matrix framework, it has been shown that an ${\cal H}^{2}$-matrix-based addition, matrix-vector product (MVP), and matrix-matrix product (MMP) all can be performed in linear complexity for constant-rank ${\cal H}^{2}$ [1]. However, the accuracy of existing ${\cal H}^{2}$-MMP algorithm like [1] is not controlled. This is because given two ${\cal H}^{2}$-matrices $\textbf{A}_{{\cal H}^{2}}$ and $\textbf{B}_{{\cal H}^{2}}$, the matrix structure and cluster bases of their product $\textbf{C}=\textbf{A}_{{\cal H}^{2}}\times\textbf{B}_{{\cal H}^{2}}$ are pre-assumed, and a formatted multiplication is performed, whose accuracy is not controlled. For example, the row cluster bases of $\textbf{A}_{{\cal H}^{2}}$ and the column cluster bases of $\textbf{B}_{{\cal H}^{2}}$ are assumed to be those of C. This treatment lacks accuracy control since the original cluster basis may not be able to represent the new matrix content generated during the MMP. For instance, when multiplying a full-matrix block F by a low rank block $\textbf{V}_{t}\textbf{S}{\textbf{V}_{s}}^{T}$, treating the result as a low-rank block is correct. However, it is inaccurate to use the original row cluster basis $\textbf{V}_{t}$ as the product’s row cluster basis, since the latter has been changed to $\textbf{F}\textbf{V}_{t}$. Therefore, the algorithm in [1] can be accurate if the cluster bases of the original matrices can also be used to accurately represent the matrix product. However, this is unknown in general applications, and hence the accuracy of existing linear-complexity MMP algorithm is not controlled. One can find many cases where a formatted multiplication would fail.

The posteriori multiplication in [2] is more accurate than the formatted multiplication in [1]. But it is only suitable for special ${\cal H}^{2}$ matrices. Besides, this posteriori multiplication requires much more computational time and memory than the formatted one. It needs to first present the product in ${\cal H}$-matrix and then convert it into an ${\cal H}^{2}$-matrix, the complexity of which is not linear.

In this work, we propose a new algorithm to do the ${\cal H}^{2}$ matrix-matrix multiplication with controlled accuracy. The cluster bases are calculated instantaneously based on the prescribed accuracy during the computation of the matrix-matrix product. Meanwhile, we are able to keep the computational complexity to be linear for constant-rank ${\cal H}^{2}$. For variable-rank cases such as those in an electrically large analysis, the proposed MMP is also efficient since it only involves $O(2^{l})$ computations at level $l$, each of which costs $O(k_{l}^{3})$ only, where $k_{l}$ is the rank at tree level $l$. This algorithm can be used as a fundamental arithmetic in the error-controlled fast inverse, LU factorization, solution for many right hand sides, etc. Numerical experiments have demonstrated its accuracy and low complexity. In [12, 13], we present a fast algorithm to compute the product of two ${\cal H}^{2}$-matrices in controlled accuracy. However, unlike this work, the original cluster bases are not completely changed, but appended to account for the updates to the original matrix during the MMP. In [8], we present the basic idea of this work. However, it is a one-page abstract. In this paper, we provide a complete algorithm together with a comprehensive analysis of its accuracy and complexity, whose validity and performance are then demonstrated by abundant numerical examples.

II PRELIMINARIES

In an ${\cal H}^{2}$-matrix [1], the entire matrix is partitioned into multilevel admissible and inadmissible blocks, where inadmissible blocks are at the leaf level, noted as $\textbf{F}_{t,s}$. An admissible matrix block $\textbf{R}_{t,s}$ satisfies the following strong admissibility condition

[TABLE]

where $\Omega_{t}$ ($\Omega_{s}$) denotes the geometrical support of the unknown set $t$ ($s$), $diam\{\cdot\}$ is the Euclidean diameter of a set, $dist\{\cdot,\cdot\}$ denotes the Euclidean distance between two sets, and $\eta$ is a positive parameter that can be used to control the admissibility condition. An admissible matrix block in an ${\cal H}^{2}$-matrix is represented as

[TABLE]

where $\textbf{V}_{t}$ ($\textbf{V}_{s}$) is called cluster basis associated with cluster $t$ ($s$), $\textbf{S}_{t,s}$ is called coupling matrix. The cluster bases V in an ${\cal H}^{2}$-matrix has a nested property. This means the cluster basis for a non-leaf cluster $t$, $\textbf{V}_{t}$, can be expressed by its two children’s cluster bases, $\textbf{V}_{t_{1}}$ and $\textbf{V}_{t_{2}}$, as

[TABLE]

where $\textbf{T}_{t_{1}}$ and $\textbf{T}_{t_{2}}$ are called transfer matrices. Because of such a nested relationship, the cluster bases only need to be stored for leaf clusters. For non-leaf clusters, only transfer matrices need to be stored. The ${\cal H}^{2}$-matrix is stored in a tree structure, with the size of leaf-level clusters denoted by $leafsize$. The number of blocks formed by a single cluster at each tree level is bounded by a constant $C_{sp}$. In an ${\cal H}^{2}$-matrix, a large matrix block consisting of F and R is called a nonleaf block NL. As an example, a four-level block cluster ${\cal H}^{2}$-tree is illustrated in Fig. 1 (a), where the green link connects a row cluster with a column cluster, which form an admissible block, and the red links are for inadmissible blocks. The resultant ${\cal H}^{2}$-matrix is shown in Fig. 1 (b), where the admissible blocks are marked in green and the inadmissible blocks are marked in red.

III Proposed ${\cal H}^{2}$ Matrix-Matrix Product Algorithm—Leaf Level

To compute $\textbf{A}_{{\cal H}^{2}}\times\textbf{B}_{{\cal H}^{2}}=\textbf{C}_{{\cal H}^{2}}$, unlike the existing ${\cal H}^{2}$ formatted MMP [1], which is recursive, we propose to perform a one-way tree traversal from leaf level all the way up to the minimum level that has admissible blocks. Here, the tree is inverted with root level at level 0. While doing the multiplications at each level, we instantaneously compute the new row and column cluster bases based on prescribed accuracy to represent the product matrix accurately. We will use the ${\cal H}^{2}$-matrices shown in Fig. 2 to illustrate the proposed algorithm, but the algorithm is valid for any ${\cal H}^{2}$-matrix. The structures of $\textbf{A}_{{\cal H}^{2}}$, $\textbf{B}_{{\cal H}^{2}}$, and $\textbf{C}_{{\cal H}^{2}}$ matrices, i.e., which block is admissible and which is inadmissible, are determined based on the admissibility condition given in (1). During the product calculation, we will keep the structure of product $\textbf{C}_{{\cal H}^{2}}$ matrix while achieving prescribed accuracy. In this section, we detail proposed algorithm for leaf-level multiplications.

We start from leaf level $(l=L)$. Let F denote an inadmissible block, which is stored as a full matrix, and R be an admissible block. At leaf level, there are in total four matrix-matrix multiplication cases, i.e.,

• •

Case-1: $\textbf{F}^{\textbf{A}}$ $\times$ $\textbf{F}^{\textbf{B}}$

• •

Case-2: $\textbf{F}^{\textbf{A}}$ $\times$ $\textbf{R}^{\textbf{B}}$

• •

Case-3: $\textbf{R}^{\textbf{A}}$ $\times$ $\textbf{F}^{\textbf{B}}$

• •

Case-4: $\textbf{R}^{\textbf{A}}$ $\times$ $\textbf{R}^{\textbf{B}}$

The resulting matrix block in C is of two kinds: First, full matrix block, denoted by $\textbf{F}^{\textbf{C}}$, marked in red in Fig. 3(c); Second, admissible block of leaf size, which could be located at leaf level, denoted by $\textbf{R}^{\textbf{C},L}$ as marked in green in Fig. 3(c); which could also appear as a subblock in the non-leaf level $l$ as marked in blue in Fig. 3(c). The blue blocks in Fig. 3(c) are only for temporary storage, which will be changed to green admissible blocks during the upper level multiplication to preserve the structure of $\textbf{C}_{{\cal H}^{2}}$ matrix. The white blocks in Fig. 3 denote those blocks that are not involved in the leaf level multiplication. Next we show how to perform each matrix-matrix multiplication based on the two kinds of target blocks.

III-A Product is an inadmissible block (full matrix) in C

If the product matrix is a full block $\textbf{F}^{\textbf{C}}$, we can perform the four cases of multiplications exactly as they are by full matrix multiplications. For the admissible leaf blocks in four cases, we convert them into full matrices and then compute products. Since the size of these matrices is of $leafsize$, a user-defined constant, the computational cost is constant for each of such computations.

III-B Product is an admissible block in C

If the product is admissible in C whether it is a leaf-level block or a subblock of a non-leaf admissible block, case-4 can be performed as it is since the product matrix is obviously admissible, which also preserves the original row and column cluster bases. In other words, the row cluster basis of A is that of C; and the column cluster basis of B is kept in C. To see this point clearly, we can write

[TABLE]

where subscripts $i$, $j$, and $k$ denote cluster index, subscript $r$ denotes the corresponding cluster is a row cluster, whereas $c$ denotes the cluster is a column cluster. For example, $\textbf{V}_{i_{r}}^{\textbf{A}}$ denotes the cluster basis of row cluster $i$ in A, and $\textbf{V}_{k_{c}}^{\textbf{B}}$ denotes the cluster basis of column cluster $k$ in B. Eqn. (4) can be written in short as

[TABLE]

in which $\textbf{S}_{i,k}^{\textbf{C}}$ is the part in between the two cluster bases, which denotes the coupling matrix of the product admissible block in C. Clearly, this case of multiplication does not change the original row and column cluster bases.

For the other three cases, in existing MMP algorithms, a formatted multiplication is performed, which is done in the same way as case-4, i.e., using the original cluster bases of A and B or pre-assumed bases as the cluster bases of the product block. This obviously can be inaccurate since cases-1, 2, and 3, if performed as they are, would result in different cluster bases in the product matrix, which cannot be assumed. Specifically, case-1 results in a different row as well as column cluster bases in the product admissible block because

[TABLE]

case-2 yields a different row cluster basis since

[TABLE]

whereas case-3 results in a different column cluster basis in the product admissible block, because

[TABLE]

If we do not update the cluster bases in the product matrix, the accuracy of the multiplication is not controllable. Therefore, in the proposed algorithm, we update row and column cluster bases for multiplication cases 1, 2, and 3 based on prescribed accuracy. We also have to do so with the nested property taken into consideration so that the computation at nonleaf levels can be performed efficiently.

For case-1, both row and column cluster bases of the product block need to be updated. For case-2, we need to use $\textbf{F}_{i,j}^{\textbf{A}}\textbf{V}_{j_{r}}^{\textbf{B}}$ to update the original row cluster basis $\textbf{V}_{i_{r}}^{\textbf{A}}$. For case-3, we need to use $(\textbf{V}_{j_{c}}^{\textbf{A}})^{T}\textbf{F}_{j,k}^{\textbf{B}}$ to update column cluster basis $\textbf{V}_{k_{c}}^{\textbf{B}}$. Since there are many case-1, 2 and 3 products encountered at the leaf level for the same row or column cluster, we develop the following algorithm to systematically update the cluster bases. In this procedure, we also have to take the computation at all nonleaf levels into consideration so that the changed cluster bases at the leaf level can be reused at the nonleaf levels. To achieve this goal, when we update the cluster basis due to the case-1, 2, and 3 multiplications associated with this cluster, not only we consider the product admissible block in the leaf level, but also the admissible blocks at all nonleaf levels. In other words, when computing $\textbf{A}_{i,j}$ multiplied by $\textbf{B}_{j,k}$, if the $\textbf{C}_{i,k}$ block is part of a non-leaf admissible block, we will take the corresponding multiplication into account to update the cluster bases. The detailed algorithms are as follows.

III-C Computation of new cluster bases in matrix product $\textbf{C}_{{\cal H}^{2}}$

First, we show how to calculate the new row cluster bases of $\textbf{C}_{{\cal H}^{2}}$. Take an arbitrary row cluster $i$ as an example, let its cluster basis in C be denoted by $\textbf{V}^{\textbf{C}}_{i_{r}}$. This cluster basis is affected by both case-1 and case-2 multiplications, as analyzed in the above. We first find all the case-1 multiplications associated with cluster $i$, i.e., all $\textbf{F}_{i,j}^{\textbf{A}}\times\textbf{F}_{j,k}^{\textbf{B}}$ whose product block $\textbf{C}_{i,k}$ is admissible. Again, notice that the $\textbf{C}_{i,k}$ can be either admissible at leaf level or be part of a non-leaf admissible block. For any cluster $i$, the number of $\textbf{F}_{i,j}^{\textbf{A}}$ is bounded by constant $C_{sp}$, since the number of inadmissible blocks that can be formed by a cluster is bounded by $C_{sp}$. For the same reason, the number of $\textbf{F}_{j,k}^{\textbf{B}}$ for cluster $j$ is also bounded by constant $C_{sp}$. Hence, the total number of $\textbf{F}_{i,j}^{\textbf{A}}\times\textbf{F}_{j,k}^{\textbf{B}}$ multiplications is bounded by $C_{sp}^{2}$, thus also a constant. Then we calculate the Gram matrix sum of these products as:

[TABLE]

in which superscript $H$ denotes a Hermitian matrix. We also find all case-2 products associated with cluster $i$, which is the number of $\textbf{F}_{i,j}^{\textbf{A}}$ formed by cluster $i$ at leaf level in $\textbf{A}_{{\cal H}^{2}}$. This is also bounded by $C_{sp}$. Since in case-2 products, $\textbf{F}_{i,j}^{\textbf{A}}$ is multiplied by an admissible block in B, and hence $\textbf{V}_{j_{r}}^{\textbf{B}}$, we compute

[TABLE]

which incorporates all of the new cluster bases information due to case-2 products.

For case-3 and case-4 multiplications, the row cluster bases of $\textbf{A}_{{\cal H}^{2}}$ matrix are kept to be those of C. So we account for the contribution of $\textbf{V}_{i_{r}}^{\textbf{A}}$ as

[TABLE]

The column space spanning $\textbf{G}_{i_{r1}}^{\textbf{C},L}$, $\textbf{G}_{i_{r2}}^{\textbf{C},L}$ and $\textbf{G}_{i_{r3}}^{\textbf{C},L}$ would be the new cluster basis of $i$, since it takes both the original cluster basis and the change to the cluster basis due to matrix products into consideration. Since the magnitude of the three matrices may differ greatly, we normalize them before summing them up so that each component is captured. We thus obtain

[TABLE]

The  $\widehat{}$  above $\textbf{G}_{i_{r1}}^{\textbf{C},L}$, $\textbf{G}_{i_{r2}}^{\textbf{C},L}$ and $\textbf{G}_{i_{r3}}^{\textbf{C},L}$ denotes a normalized matrix. We then perform an SVD on $\textbf{G}_{i_{r3}}^{\textbf{C},L}$ to obtain the row cluster bases for cluster $i$ of $\textbf{C}_{{\cal H}^{2}}$ based on prescribed accuracy $\epsilon_{trunc}$. The singular vectors whose normalized singular values are greater than $\epsilon_{trunc}$ make the new row cluster basis $\textbf{V}_{i_{r}}^{\textbf{C}}$. It can be used to accurately represent the admissible blocks related to cluster $i$ in $\textbf{C}_{{\cal H}^{2}}$. Here, notice that the proposed algorithm for computing matrix-product cluster bases keeps nested property of $\textbf{V}_{i_{r}}^{\textbf{C}}$. This is because the Gram matrix sums in (9), (10) and (11) take the upper level admissible products into account.

To compute the column cluster bases in $\textbf{C}_{{\cal H}^{2}}$, the steps are similar to the row cluster basis computation. We account for the contributions from all the four cases of products to compute column cluster bases. As can be seen from (25) and (27), in case-1 and case-3 products, the column cluster bases are changed from the original ones; whereas in case-2 and case-4 products, the column cluster bases are kept the same as those in B.

Consider an arbitrary column cluster $k$ in $\textbf{C}_{{\cal H}^{2}}$. We find all of the case-1 products associated with $k$, which is $\textbf{F}_{i,j}^{\textbf{A}}\times\textbf{F}_{j,k}^{\textbf{B}}$ with target $\textbf{C}_{i,k}$ being admissible either at the leaf or non-leaf level. The number of such multiplications is bounded by $C_{sp}^{2}$. We then compute the sum of their Gram matrices as:

[TABLE]

Here, the superscript ∗ denotes a complex conjugate. We also find all of the case-3 products associated with $k$, which is $\textbf{R}_{i,j}^{\textbf{A}}\times\textbf{F}_{j,k}^{\textbf{B}}$ with target $\textbf{C}_{i,k}$ being admissible either at the leaf or non-leaf level. Hence, the new column cluster basis takes a form of $(\textbf{V}_{j_{c}}^{\textbf{A}})^{T}\times\textbf{F}_{j,k}^{\textbf{B}}$. The number of such multiplications is also bounded by $C_{sp}$. The sum of their Gram matrices can be computed as:

[TABLE]

For case-2 and case-4 products, the original column cluster bases of $\textbf{B}_{{\cal H}^{2}}$ are kept in $\textbf{C}_{{\cal H}^{2}}$, hence, we compute

[TABLE]

We also normalize these three Gram matrices $\textbf{G}_{k_{c1}}^{\textbf{C},L}$, $\textbf{G}_{k_{c2}}^{\textbf{C},L}$ and $\textbf{G}_{k_{c3}}^{\textbf{C},L}$ and sum them up as:

[TABLE]

We then perform an SVD on this $\textbf{G}_{k_{c}}^{\textbf{C},L}$ and truncate the singular values based on prescribed accuracy $\epsilon_{trunc}$ to obtain the column cluster bases $\textbf{V}_{k_{c}}^{\textbf{C}}$ for cluster $k$. Now this new column cluster basis $\textbf{V}_{k_{c}}^{\textbf{C}}$ can be used to accurately represent the admissible blocks formed by column cluster $k$ in $\textbf{C}_{{\cal H}^{2}}$.

III-D Computation of the four cases of multiplications with the product block being admissible

After computing the new row and column cluster bases of the product matrix, for the multiplication whose target is an admissible block described in Section III-B, the computation becomes the coupling matrix computation since the cluster bases have been generated. For the four cases of multiplications, their coupling matrices have the following expressions:

[TABLE]

The resulting admissible blocks in $\textbf{C}_{{\cal H}^{2}}$ are nothing but $\textbf{R}_{i,k}^{\textbf{C}}=\textbf{V}_{i_{r}}^{\textbf{C}}\times\textbf{S}_{i,k}^{\textbf{C}}\times(\textbf{V}_{k_{c}}^{\textbf{C}})^{T}$.

In (17), the $\textbf{B}_{j}$ is the cluster bases product, which is as shown below:

[TABLE]

Since it is only related to the original cluster bases, it can be prepared in advance before the MMP computation. Using the nested property of the cluster bases, $\textbf{B}_{j}$ can be computed in linear time for all clusters $j$, be $j$ a leaf or a non-leaf cluster.

In (17), the $(\textbf{V}_{i_{r}}^{\textbf{C}})^{H}\textbf{V}_{i_{r}}^{\textbf{A}}$ is simply the projection of the original row cluster basis of A onto the new cluster basis of the product matrix C. Similarly, $(\textbf{V}_{k_{c}}^{\textbf{B}})^{T}(\textbf{V}_{k_{c}}^{\textbf{C}})^{*}$ denotes the projection of the original column cluster basis of B onto the newly generated column cluster basis in C. The two cluster basis projections can also be computed for every leaf cluster after the new cluster bases have been generated. Hence, we compute

[TABLE]

for each leaf row cluster $i$, and each column leaf cluster $k$. In this way, it can be reused without recomputation for each admissible block formed by $i$ or $k$.

In (17), we can also see that the F block is front and back multiplied by cluster bases. It can be viewed as an F block collected based on the front (row) and back (column) cluster bases, which becomes a matrix of rank size. Specifically, in (17), there are three kinds of collected blocks

[TABLE]

which is used in case-1, 2, and 3 multiplication respectively.

As can be seen from (17), the case-1 multiplication with an admissible block being the target can be performed by first computing the full-matrix product, and then collecting the product onto the new row and column cluster bases of the product matrix. This collect operation is accurate because the newly generated row and column cluster bases have taken such a case-1 multiplication into consideration when being generated. As for the case-2 multiplication, as can be seen from (17), we can use the $\textbf{F}_{i,j}$ collected based on the new row cluster basis and the original column cluster basis, the size of which is rank, to multiply the coupling matrix of $\textbf{S}_{j,k}$, and then multiply the column basis projection matrix since the column bases have been changed. Similarly, for case-3, we use the collected block $(\textbf{F}_{j,k}^{\textbf{B}})_{coll.}$, and front multiply it by the coupling matrix of $\textbf{S}_{i,j}$, and then front multiply a row cluster basis transformation matrix. As for case-4, we multiply the coupling matrix of A’s admissible block by the cluster basis product, and then by the coupling matrix of B’s admissible block. Since the row and column cluster bases have been changed to account for the other cases of multiplications, at the end, we need to front and back multiply the cluster basis transformation matrices to complete the computation of case-4. Summarizing the aforementioned, the coupling matrix in (17) can be efficiently computed as

[TABLE]

III-E Summary of overall algorithm at leaf level

Here, we conclude all the operations related to leaf level computation when the target is an admissible block:

Prepare cluster bases product B; 2. 2.

Compute all the leaf-level row and column cluster bases of product matrix $\textbf{C}_{{\cal H}^{2}}$; 3. 3.

Collect the F blocks in $\textbf{A}_{{\cal H}^{2}}$ and $\textbf{B}_{{\cal H}^{2}}$ based on the new row and/or column cluster bases, also prepare cluster bases transformation matrix P ; 4. 4.

Perform four cases of multiplications.

After leaf level multiplications, we need to merge four coupling matrices at a non-leaf level admissible block, as shown by the blue blocks in Fig. 3 (c). These matrices correspond to the multiplication case of a nonleaf block NL multiplied by a nonleaf block NL generating an admissible block at next level. The merged block is the coupling matrix of this next-level admissible block. It will be used to update next level transfer matrices. The details will be given in next section.

IV Proposed ${\cal H}^{2}$ Matrix-Matrix Product Algorithm—Non-Leaf Level

After finishing the leaf level multiplication, we proceed to non-leaf level multiplications. In Fig. 4, we use level $L-1$ as an example to illustrate $\textbf{A}_{{\cal H}^{2}}^{L-1}$, $\textbf{B}_{{\cal H}^{2}}^{L-1}$, and $\textbf{C}_{{\cal H}^{2}}^{L-1}$.

At a nonleaf level $l$, there are also in total four matrix-matrix multiplication cases, i.e.,

• •

Case-1: $\textbf{NL}^{\textbf{A}}$ $\times$ $\textbf{NL}^{\textbf{B}}$

• •

Case-2: $\textbf{NL}^{\textbf{A}}$ $\times$ $\textbf{R}^{\textbf{B}}$

• •

Case-3: $\textbf{R}^{\textbf{A}}$ $\times$ $\textbf{NL}^{\textbf{B}}$

• •

Case-4: $\textbf{R}^{\textbf{A}}$ $\times$ $\textbf{R}^{\textbf{B}}$,

where NL denotes a non-leaf block. The resulting matrix block in C is also of two kinds: 1) non-leaf block NL at this level, marked in red in Fig. 4 (c), and 2) admissible block R, marked in green in Fig. 4 (c). Next we show how to perform each case of multiplications based on the two kinds of target blocks.

IV-A Product is an NL block in C

The NL target block would not exist for a case-1 multiplication, since if a case-1 multiplication results in an NL block, that computation should have been performed at previous level. As for the other three cases of multiplications, since at least one admissible block is present in the multipliers, the product must be an admissible block. Hence, we compute them as having an admissible block as the product, using the algorithm described in the following subsection, and associate the resulting admissible block with the NL block. After the computation is done at all levels, we perform a backward split operation to split the admissible block associated with each NL block to each leaf block of C based on its structure.

IV-B Product is an admissible block in C

Similar to the leaf level, if the product is an admissible block in C whether at the same non-leaf level or at an upper level, case-4 can be performed as it is since the product matrix is obviously admissible, which also preserves the original row and column cluster bases. We can write

[TABLE]

where T denotes a transfer matrix, and superscript $ch$ denotes the two children clusters of the non-leaf cluster $i$. If the cluster bases at leaf level and the transfer matrices at non-leaf levels are kept the same as before, then the computation of (IV-B) is to calculate the coupling matrix at level $l$, which is

[TABLE]

It is a product of three small matrices whose size is the rank at this tree level. Rewriting (IV-B) as

[TABLE]

If we exclude the children cluster bases in the front and at the back, we can see that T serves as the new cluster basis at this level. In other words, at a non-leaf level $l$, if we treat this level as the bottom level of the remaining tree, then the transfer matrix of the non-leaf cluster is nothing but the leaf cluster basis of the shortened tree.

Similar to the leaf-level computation, the other three cases of multiplications will result in a change of cluster basis in the matrix product. Specifically, case-1 results in a different row as well as column cluster bases in the product admissible block because

[TABLE]

case-2 yields a different row cluster basis since

[TABLE]

whereas case-3 results in a different column cluster basis in the product admissible block, because

[TABLE]

If we do not update the cluster bases in the product matrix, the accuracy of the multiplication is not controllable. However, if we update the cluster basis as they are, it is computationally very expensive since the matrix block size keeps increasing when we proceed from leaf level towards the root level. In addition to the cost of changing cluster bases, if we have to carry out the multiplications at each non-leaf level using the actual matrix block size, then the computation is also prohibitive. Therefore, the fast algorithm we develop here is to perform all computations using the rank size at each tree level, and meanwhile control the accuracy.

In the proposed algorithm, to account for the updates to the original matrix during the MMP procedure, the cluster bases of C are computed level by level, which are manifested by the changed leaf cluster bases and the transfer matrices at nonleaf levels. At a non-leaf level, its children-level cluster bases have already been computed, and they are different from the original ones in A and B. However, the new cluster bases have taken the upper-level multiplications into consideration. Hence, we can accurately represent the multiplication at the current non-leaf level using newly generated children cluster bases.

Take case-1 product as an example, where we perform $\textbf{NL}_{i,j}^{\textbf{A}}\times\textbf{NL}_{j,k}^{\textbf{B}}$ obtaining an admissible $\textbf{R}_{i,k}^{\textbf{C}}$. We can accurately represent this product using the children cluster bases of $i$ and $k$ as follows:

[TABLE]

in which

[TABLE]

This collected block, $(\textbf{NL}_{i,j}^{\textbf{A}}\textbf{NL}_{j,k}^{\textbf{B}})_{coll.}$, is actually the coupling matrix merged from the four small coupling matrices computed at previous level, when dealing with the multiplication case of having a target block as a subblock in the upper-level admissible block. It can be written as

[TABLE]

Each of the four coupling matrices has been obtained at previous level. From (IV-B), it is clear that using the nested property of the cluster bases, the collect operation does not need to start from leaf level, but using the four blocks obtained at previous one level.

For case-2 product, it can also be accurately expanded in the space of the children row cluster bases, and hence

[TABLE]

where

[TABLE]

which is ${\textbf{NL}_{i,j}^{\textbf{A}}}$ collected based on the children’s new row cluster bases in C and the original column cluster bases in B. From (IV-B), it can be seen that if excluding the children cluster bases, then ${\textbf{NL}_{i,j}^{\textbf{A}}}_{coll.}\textbf{T}_{j_{r}}^{\textbf{B}}$ resembles the $\textbf{F}_{i,j}^{\textbf{A}}\textbf{V}_{j_{r}}^{\textbf{B}}$ in the leaf level case-2 product. In other words, if we treat the current non-leaf level as the leaf level, then ${\textbf{NL}_{i,j}^{\textbf{A}}}_{coll.}$ is equivalent to a full matrix block, whereas T is the leaf cluster basis. An example of ${\textbf{NL}_{i,j}^{\textbf{A}}}_{coll.}$ block at level $(L-1)$ in $\textbf{A}_{{\cal H}^{2}}$ can be seen below:

[TABLE]

which consists of collected full matrices whose expressions are shown in (20), and projected coupling matrices of admissible blocks. Again, using the nested property of both new and original cluster bases, the collect operation does not need to start from leaf level, but using the four blocks obtained at previous one level. Each collect operation only costs $O(k_{l})^{3}$, where $k_{l}$ is the rank at level $l$.

Since the cluster bases at the previous level have been computed, for case-1 and case-2 products at a non-leaf level, we only need to compute the center block associated with the current non-leaf level, and this computation can be carried out in the same way as how we carry out leaf-level computation, if we treat the current non-leaf level as the leaf level of the remaining tree. The same is true to case-3 product, where we have

[TABLE]

in which

[TABLE]

We can see that $(\textbf{T}_{j_{c}}^{\textbf{A}})^{T}{\textbf{NL}_{j,k}^{\textbf{B}}}_{coll.}$ resembles the $(\textbf{V}_{j_{c}}^{\textbf{A}})^{T}\textbf{F}_{j,k}^{\textbf{B}}$ in the leaf level case-3 product. An example of collected NL block in $\textbf{B}_{{\cal H}^{2}}$ is given as follows

[TABLE]

which consists of collected full matrices whose expressions are shown in (20), and projected coupling matrices of admissible blocks.

Since the cluster bases have been changed at previous level, we also represent the case-4 product using the new children cluster bases of $i$ and $k$, thus

[TABLE]

and

[TABLE]

which can be written in short as

[TABLE]

where $ch$ denotes children. Here, there is a cluster basis transformation matrix in the front and at the back.

IV-C Computation of the new non-leaf level transfer matrices in C

If the target block is an admissible block at a nonleaf level, we need to represent it as $\textbf{R}_{t,s}=\textbf{T}_{t}\textbf{S}_{t,s}(\textbf{T}_{s})^{T}$ in controlled accuracy. Hence, we need to calculate new row and column transfer matrices T of product matrix $\textbf{C}_{{\cal H}^{2}}$. First, we introduce how to calculate the row transfer matrices. Similar to leaf level, case-1 and 2 products result in a change in the row cluster basis and hence row transfer matrix. Case-3 and 4 products do not require a change of transfer matrix if the cluster bases have not been changed at previous level. However, since the cluster bases have been changed at previous level, the transfer matrix requires an update as well.

For an arbitrary non-leaf cluster $i$, we first find all of the case-1 products associated with $i$. Each of such a product leads to a coupling matrix merged from the four coupling matrices obtained at previous level computation, denoted by $(\textbf{NL}_{i,j}^{\textbf{A}}\textbf{NL}_{j,k}^{\textbf{B}})_{coll.}$. Using them, we calculate the Gram matrix sum as:

[TABLE]

The second step is to take case-2 multiplications at a non-leaf level into consideration for row transfer matrix calculation of product matrix $\textbf{C}_{{\cal H}^{2}}$. We find all the collected nonleaf blocks ${\textbf{NL}_{i,j}^{\textbf{A}}}_{coll.}$ of cluster $i$ at level $l$ in $\textbf{A}_{{\cal H}^{2}}$ matrix and multiply them with corresponding transfer matrices $\textbf{T}_{j_{r}}^{\textbf{B}}$ from $\textbf{B}_{{\cal H}^{2}}$ matrix. And we calculate the Gram matrix sum as

[TABLE]

Finally, we count the contributions from case-3 and case-4 products by computing

[TABLE]

Again, we normalize these three Gram matrices and obtain

[TABLE]

We then calculate an SVD of this $\textbf{G}_{i_{r}}^{\textbf{C},l}$ and truncate the singular values based on prescribed accuracy $\epsilon_{trunc}$ to obtain row transfer matrix $\textbf{T}_{i_{r}}^{\textbf{C}}$ for cluster $i$ at nonleaf level.

Similarly, we can compute the new column transfer matrices for non-leaf cluster $k$, which is $\textbf{T}_{k_{c}}^{\textbf{C}}$. The first part is

[TABLE]

The second part is

[TABLE]

The third part is

[TABLE]

Then we normalize the three Gram matrices and sum them up as

[TABLE]

After we perform an SVD on $\textbf{G}_{k_{c}}^{\textbf{C},l}$ matrix and truncate the singular values based on prescribed accuracy $\epsilon_{trunc}$, we get new column transfer matrix $\textbf{T}_{k_{c}}^{\textbf{C}}$.

IV-D Computation of the four cases of multiplications with the product block being admissible

Now we obtain both row and column transfer matrices for product matrix $\textbf{C}_{{\cal H}^{2}}$, hence, the four multiplications become the computation of the coupling matrices, so that the admissible block at the current level has a form of $\textbf{R}_{t,s}=\textbf{T}_{t}\textbf{S}_{t,s}(\textbf{T}_{s})^{T}$. The coupling matrix S’s calculation is similar to that of leaf level in (17), which has the following expressions:

[TABLE]

Again, we should prepare some matrix products in advance so that we can achieve linear complexity MMP for constant rank ${\cal H}^{2}$-matrix. For nonleaf levels, the cluster bases product $\textbf{B}_{j}$ can be readily calculated using children’s cluster bases based on the nested property. For example, given a nonleaf cluster $j$, we can generate $\textbf{B}_{j}$ by using the cluster bases product of its children clusters $j_{1}$ and $j_{2}$, which is shown as:

[TABLE]

Besides, since the cluster bases product $\textbf{B}_{j}$ only involve original cluster bases in $\textbf{A}_{{\cal H}^{2}}$ and $\textbf{B}_{{\cal H}^{2}}$ matrices, we can prepare the above $\textbf{B}_{j}$ for all leaf and nonleaf clusters before MMP algorithm. In addition, the nonleaf level cluster bases projection (transformation) can also be calculated using children’s ones as shown in (19). The formulas are given below:

[TABLE]

We also compute the collected NL matrix block in $\textbf{A}_{{\cal H}^{2}}$ and $\textbf{B}_{{\cal H}^{2}}$ at current level $l$ by the following equation,

[TABLE]

where superscript $l$ denotes tree level. After we prepare the matrix products in (49), (50), and (51), we can proceed to calculate the coupling matrices in (48) efficiently as:

[TABLE]

All the coupling matrices calculation are performed in rank size $k_{l}$. So the computational cost is $O(k_{l}^{3})$. After coupling matrices calculation in (52), all the admissible products at this nonleaf level multiplication can be represented as $\textbf{R}_{i,j}^{\textbf{C}}=\textbf{T}_{i_{r}}^{\textbf{C}}\textbf{S}_{i,j}^{\textbf{C}}\textbf{T}_{j_{c}}^{\textbf{C}}$.

IV-E Summary of overall algorithm at each non-leaf level

The cluster bases products $\textbf{B}_{j}$ have been computed for all clusters $j$ before the MMP starts, since they are only related to the original cluster bases.

At each non-leaf level, we do the following:

Collect four blocks in an NL block in $\textbf{A}_{{\cal H}^{2}}$ to a block of $O(k_{l+1})$ size, using the newly generated children row cluster bases of C (transfer matrices if children are not at the leaf level) and the original column cluster bases of B (or transfer matrices). This is to generate the $(\textbf{NL}_{i,j}^{\textbf{A}})_{coll.}$, shown in (32). 2. 2.

Collect four blocks in an NL block in $\textbf{B}_{{\cal H}^{2}}$ to a block of $O(k_{l+1})$ size, using the original row cluster bases of A (or transfer matrices) and the new children column cluster bases of C (transfer matrices if children are not at the leaf level). This is to generate $(\textbf{NL}_{j,k}^{\textbf{B}})_{coll.}$, shown in (35). 3. 3.

Merge four blocks in an R block in $\textbf{C}_{{\cal H}^{2}}$. This corresponds to the $(\textbf{NL}_{i,j}^{\textbf{A}}\textbf{NL}_{j,k}^{\textbf{B}})_{coll.}^{(l)}$ in (52). 4. 4.

Calculate new row and column transfer matrices of product matrix $\textbf{C}_{{\cal H}^{2}}$ at this level. 5. 5.

Prepare cluster bases projections $\textbf{P}_{i}^{\textbf{A}}$, $\textbf{P}_{k}^{\textbf{B}}$, and perform an NL block collect shown in (51); 6. 6.

Perform four cases of multiplications shown in (52).

After we finish one-way bottom-up tree traversal to calculate block matrix products at all the levels, i.e. from leaf level all the way up to minimal admissible level, we need to perform a post-processing for the coupling matrices associated with the NL blocks in $\textbf{C}_{{\cal H}^{2}}$. They exist because of the multiplications cases described in Section IV-A. This could be efficiently done by performing one-way top-down split process, the same as the matrix backward transformation shown in [1]. This post processing stage is to split the coupling matrices in NL to lower level admissible or inadmissible blocks.

V Accuracy and Complexity Analysis

In this section, we analyze the accuracy and computational complexity of the proposed algorithm to compute ${\cal H}^{2}$-matrix-matrix products.

V-A Accuracy

Different from existing formatted ${\cal H}^{2}$-matrix-matrix products [1], in the proposed new algorithm, the accuracy of the product is directly controlled by $\epsilon_{trunc}$. No formatted multiplications are performed, and the cluster bases are changed to represent the updates to the original matrix accurately. This makes each operation performed in the proposed MMP controlled by accuracy or exact. When generating an ${\cal H}^{2}$-matrix to represent the original dense matrix, the accuracy is controlled by $\epsilon_{{\cal H}^{2}}$, which is the same as in [11].

V-B Time and Memory Complexity

The proposed MMP involves $O(L)$ levels of computation. At each level, there are $2^{l}$ clusters. For each cluster, the cost of changing the cluster bases at the leaf level due to four cases of multiplications is to perform $O(C_{sp})^{2}$ multiplications, and each of which has a constant cost, as can be seen from (13), (14), and (15). The cost of changing the cluster bases at the non-leaf level due to the four cases of multiplications is also to perform $O(C_{sp})^{2}$ multiplications for each cluster, and each of which has a cost of $O(k_{l})^{3}$, as can be seen from (40), (41), (42). Notice that the NL blocks in A and B are collected level by level, at each level, there are $2^{l}O(C_{sp})$ NL block, and each collect operation also costs $O(k_{l})^{3}$ only. Other auxiliary matrices are generated using a similar computational cost.

As for the computation of the four cases of multiplications at each level, each case involves $O(C_{sp})^{2}$ multiplications for each cluster, and each of which costs $O(k_{l})^{3}$ at the non-leaf level and ${O(leafsize)}^{3}$ at the leaf level as can be seen from (21), and (52).

Hence, the time complexity of the proposed MMP can be found as

[TABLE]

And the storage for each block is $O(k_{l}^{2})$, with each cluster having $C_{sp}$ blocks. So the memory complexity is

[TABLE]

Recall $k_{l}$ is the rank at tree level $l$. Hence, (53) and (54) show that the overall complexity is a function of rank $k_{l}$. Taking into account the rank’s growth with electrical size as shown in [5], we can get the time and memory complexity of proposed MMP for different rank scaling. For constant-rank ${\cal H}^{2}$-matrices, since $k_{l}$ is a constant irrespective of matrix size, the complexity of the proposed direct solution is strictly $O(N)$ in both CPU time and memory consumption, as shown below.

For constant $k_{l}$:

[TABLE]

[TABLE]

For electrodynamic analysis, to ensure a prescribed accuracy, the rank becomes a function of electrical size, and thereby tree level. Different ${\cal H}^{2}$-matrix representations can result in different complexities, because their rank’s behavior is different. Using a minimal-rank ${\cal H}^{2}$-representation, as shown by [5], the rank grows linearly with electrical size for general 3-D problems. In a VIE, $k_{l}$ is proportional to the cubic root of matrix size at level $l$, because this is the electrical size at level $l$. Hence for a VIE, (53) and (54) become

For $k_{l}$ linearly growing with electrical size:

[TABLE]

[TABLE]

So the time complexity of the proposed MMP algorithm for 3D electrodynamic analysis is $O(NlogN)$, and the memory complexity is $O(N)$.

VI Numerical Results

In order to demonstrate the accuracy and low computational complexity of the proposed fast ${\cal H}^{2}$-matrix-matrix multiplication for general ${\cal H}^{2}$-matrices, we use ${\cal H}^{2}$-matrices resulting from large-scale capacitance extraction and volume integral equation (VIE) solvers for electromagnetic analysis as examples. The capacitance extraction matrix is shown in [7]. The VIE formulation is based on [9] with SWG vector bases for expanding electric flux density in each tetrahedral element. A variety of large-scale examples involving over one million unknowns are simulated on a single CPU core to examine the accuracy and complexity of the proposed MMP algorithm. The ${\cal H}^{2}$-matrix for each example is constructed based on the method described in [10, 11]. The capacitance matrix is used to demonstrate the proposed MMP algorithm performance for constant-rank ${\cal H}^{2}$-matrices. We also simulate large scale 2- and 3-D scattering examples to examine the time and memory complexity of the proposed MMP for variable rank cases. The computer used has an Intel(R) Xeon(R) CPU E5-2690 v2 running at 3 GHz, and only a single core is employed to carry out the computation.

VI-A Two-layer Cross Bus

The first example is the capacitance extraction of a 2-layer cross bus structure. In each layer, there are $m$ conductors, and each conductor has a dimension of $1\times 1\times(2m+1)$ m3. We simulate a suite of such structures with 16, 32, 64, 128, and 256 buses respectively. The parameters used in the ${\cal H}^{2}$-matrix construction are $leafsize=30$, admissibility condition [1] $\eta=1.0$, and $\epsilon_{{\cal H}^{2}}=10^{-4}$. For the proposed ${\cal H}^{2}$ MMP, the $\epsilon_{trunc}$ is chosen to be $10^{-2}$, $10^{-4}$ and $10^{-6}$ respectively to examine the error controllability. As shown in Fig. 5, the proposed MMP exhibit clear linear complexities in time and memory regardless of the choice of $\epsilon_{trunc}$. Certainly, the smaller the $\epsilon_{trunc}$, the larger the computational cost.

The accuracy of the proposed MMP is assessed by using the following criterion:

[TABLE]

where $\textbf{A}_{{\cal H}^{2}}\times(\textbf{B}_{{\cal H}^{2}}\times x)$ is used as the reference solution, since given an ${\cal H}^{2}$ matrix, a matrix-vector product can be carried out without any approximation as shown in [1]. In generating the reference solution, we first compute $y=\textbf{B}_{{\cal H}^{2}}\times x$, and then compute $\textbf{A}_{{\cal H}^{2}}\times y$, both of which are done in exact arithmetic. The proposed solution is generated by first computing an MMP of $\textbf{A}_{{\cal H}^{2}}\textbf{B}_{{\cal H}^{2}}$ to obtain $\textbf{C}_{{\cal H}^{2}}$, and then compute $\textbf{C}_{{\cal H}^{2}}x$. From Table I, we can see the accuracy of the proposed MMP is good, and it is also controllable.

VI-B Large-scale Dielectric Slab Scattering

We then simulate a dielectric slab with $\epsilon_{r}=2.54$ at 300 MHz. The thickness of the slab is fixed to be $0.1\lambda_{0}$. The width and length are simultaneously increased from $4\lambda_{0}$, $8\lambda_{0}$, $16\lambda_{0}$, to $28\lambda_{0}$. With a mesh size of $0.1\lambda_{0}$, the resultant $N$ ranges from 22,560 to 1,098,720 for this suite of slab structures. The parameters used in the ${\cal H}^{2}$-matrix construction are $leafsize=40$, admissibility condition [1] $\eta=2.0$, and $\epsilon_{{\cal H}^{2}}=10^{-3}$. For the proposed ${\cal H}^{2}$ MMP, the $\epsilon_{trunc}$ is chosen to be $10^{-2}$, $10^{-3}$, $10^{-4}$, $10^{-5}$and $10^{-6}$ respectively, to examine the computational complexity and error controllability of the proposed MMP. Based on [5], the rank’s growth rate with electrical size for 2-D slab is lower than linear, and being a square root of the log-linear of the electric size. Substituting such a rank’s growth into the complexity analysis in (53) and (54), we obtain linear complexity in both memory and time.

In Fig. 6 (a), we plot the MMP time with respect to $N$, for all different choices of $\epsilon_{trunc}$. It is clear that the smaller $\epsilon_{trunc}$ value, the larger the MMP time. However, the complexity remains the same as linear regardless of the choice of $\epsilon_{trunc}$. The memory cost is plotted in Fig. 6 (b). Obviously, it scales linearly with the number of unknowns. The error of the proposed MMP is measured in the same way as shown in (59). In Table II, we list the error as a function of $\epsilon_{trunc}$. Excellent accuracy can be observed in the entire unknown range. Furthermore, the accuracy can be controlled by $\epsilon_{trunc}$, and overall smaller $\epsilon_{trunc}$ results in better accuracy.

VI-C Scattering from Large-scale Array of Dielectric Cubes

Next, we simulate a large-scale array of dielectric cubes at 300 MHz. The relative permittivity of the cube is $\epsilon_{r}=4.0$. Each cube is of size $0.3\lambda_{0}\times 0.3\lambda_{0}\times 0.3\lambda_{0}$. The distance between adjacent cubes is kept to be $0.3\lambda_{0}$. The number of the cubes is increased along the $x$-, $y$-, and $z$- directions simultaneously from 2 to 16, thus producing a 3-D cube array from $2\times 2\times 2$ to $16\times 16\times 16$ elements. The number of unknowns $N$ is respectively 3,024, 24,192, 193,536, and 1,548,288 for these arrays. During the construction of ${\cal H}^{2}$-matrix, we set $leafsize=20$, $\eta=1$ and $\epsilon_{{\cal H}^{2}}=10^{-2}$. For the proposed ${\cal H}^{2}$ MMP, the $\epsilon_{trunc}$ is chosen as $10^{-2}$, $10^{-3}$ and $10^{-4}$.

For a cubic growth of unknowns in 3-D problems, we observe that constant $C_{sp}$ is quite different for different unknowns, as can be seen from Table III. It is thus important to analyze the performances of the proposed MMP as $Memory/C_{sp}$ and $Multiplication\ time/C_{sp}^{2}$ respectively to examine the true scaling rate. In Fig. 7 (a) and Fig. 7 (b), we plot the ${\cal H}^{2}$-matrix-matrix multiplication time divided by $C_{sp}^{2}$, and the storage cost normalized with $C_{sp}$ with respect to N. As can be seen, their scaling rate with $N$ agrees very well with our theoretical complexity analysis. For the over one-million unknown case which is a $16\times 16\times 16$ cube array having thousands of cube elements, the error is still controlled to be as small as $0.809\%$ using $\epsilon_{trunc}=10^{-4}$. The error of the proposed MMP is listed in Table IV for this example, which again reveals excellent accuracy and error controllability of the proposed MMP.

We also compare the accuracy of the proposed MMP with existing MMP [1] using this 3-D example. As shown in Table IV, the proposed MMP has much better accuracy, and also it is controllable.

VII Conclusions

In this paper, we develop a fast accuracy-controlled algorithm to compute ${\cal H}^{2}$-matrix-matrix products for general ${\cal H}^{2}$-matrices. This proposed algorithm not only has explicitly controlled accuracy, but also generates a rank-minimized representation of the product matrix based on prescribed accuracy. The row and column cluster bases are instantaneously changed so that the new matrix content generated during the MMP can be accurately represented. This ensures that each multiplication performed in the proposed MMP is well controlled by accuracy. Meanwhile, we retain the complexity to be linear for constant-rank ${\cal H}^{2}$-matrices. The proposed algorithm has been applied to calculate ${\cal H}^{2}$-matrix-matrix products for large-scale capacitance extraction matrices whose kernel is static and real-valued and electrically large VIEs whose kernel is oscillatory and complex-valued. For constant-rank ${\cal H}^{2}$-matrices, the proposed MMP has an $O(N)$ complexity in both time and memory. For rank growing with the electrical size linearly, the proposed MMP has an $O(NlogN)$ complexity time and $O(N)$ complexity in memory. ${\cal H}^{2}$-matrix products with millions of unknowns are simulated on a single core CPU in fast CPU run time. Comparisons with existing ${\cal H}^{2}$-matrix-matrix product algorithm have demonstrated clear advantages of the proposed new MMP algorithm.