A Hierarchical Butterfly LU Preconditioner for Two-Dimensional Electromagnetic Scattering Problems Involving Open Surfaces

TL;DR

This paper presents a hierarchical butterfly LU preconditioner that significantly accelerates the solution of large-scale 2D electromagnetic scattering problems involving open surfaces, with scalable efficiency and high accuracy.

Contribution

It introduces the H-IDBF-LU preconditioner using interpolative decomposition butterfly factorization for efficient iterative solutions of EFIEs on open surfaces.

Findings

Memory and computational cost scale as O(N log^2 N) per iteration.

Number of iterations remains bounded or grows slowly with problem size.

Successfully applied to problems with up to 100 million unknowns.

Abstract

This paper introduces a hierarchical interpolative decomposition butterfly-LU factorization (H-IDBF-LU) preconditioner for solving two-dimensional electric-field integral equations (EFIEs) in electromagnetic scattering problems of perfect electrically conducting objects with open surfaces. H-IDBF-LU leverages the interpolative decomposition butterfly factorization (IDBF) to compress dense blocks of the discretized EFIE operator to expedite its application; this compressed operator also serves as an approximate LU factorization of the EFIE operator leading to an efficient preconditioner in iterative solvers. Both the memory requirement and computational cost of the H-IDBF-LU solver scale as $O(N\log^2 N)$ in one iteration; the total number of iterations required for a reasonably good accuracy scales as $O(1)$ to $O(\log^2N)$ in all of our numerical tests. The efficacy and accuracy of the proposed preconditioned iterative solver are demonstrated via its application to a broad range of scatterers involving up to $100$ million unknowns.

Tables

Table 1. Table 1: Measured solution error for different geometry shapes.

shape	semicircle	corner	spiral	strips	square	cup
$N$=5E3	2.24E-06	9.51E-06	8.13E-06	7.12E-05	2.28E-05	1.60E-05
$N$=5E4	1.11E-05	9.84E-06	3.82E-05	6.45E-04	2.24E-04	1.84E-04
$N$=5E5	5.86E-06	3.85E-06	4.01E-05	9.86E-04	2.38E-04	3.56E-04
$N$=5E6	1.10E-05	8.17E-06	1.37E-04	3.74E-04	2.33E-04	6.16E-04

Equations

$$K\approx U^L{U}^{L-1}\cdots U^h{S}^h{V}^h\cdots V^{L-1}{V}^L,$$

$$E_z^{inc}(\rho )=\frac{\kappa {\eta }_0}{4}{\int }_S{J}_z({\rho }^{\prime })H_0^{(2)}(\kappa |\rho -{\rho }^{\prime }|)ds,\forall \rho \in S.$$

$$Ax=b,$$

$$A_{ij}$$

$$\tilde{L}(i,j)$$

$$\tilde{U}(i,j)$$

$${\tilde{L}}^{-1}A{\tilde{U}}^{-1}y={\tilde{L}}^{-1}b,$$

$$\tilde{U}x=y.$$

$$L=\left(\begin{array}{cc}{F}_{11}& 0\\ F_{21}& F_{22}\end{array}\right);$$

$$U=\left(\begin{array}{cc}{F}_{11}& F_{12}\\ 0& F_{22}\end{array}\right);$$

$$A=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]=\left[\begin{array}{cc}{L}_{11}& \\ A_{21}{U}_{11}^{-1}& I\end{array}\right]\left[\begin{array}{cc}I& \\ & S_{22}\end{array}\right]\left[\begin{array}{cc}{U}_{11}& L_{11}^{-1}{A}_{12}\\ & I\end{array}\right],$$

$$A_{ij}$$

$$H_0^{(2)}(x)=\sqrt{\frac{2}{\pi x}}\frac{1}{\sqrt{{\alpha }^{\prime }(x)}}{e}^{-i\alpha (x)},$$

$$B_{ij}$$

$$C_{ij}$$

$$x\in [x_0-\frac{d_x}{2N},x_0+\frac{d_x}{2N}]\text{and}y\in [y_0-\frac{d_y}{2N},y_0+\frac{d_y}{2N}]$$

$$d_x{d}_y\le M\text{and}(x_0-\frac{d_x}{2N})-(y_0+\frac{d_y}{2N})\ge O(\frac{M}{N}).$$

$$e^{-i\alpha (\kappa |\rho (x)-\rho (y)|)}=\beta (x)e^{-i{\partial }_{xy}\alpha (\kappa |\rho (\bar{x})-\rho (\bar{y})|)(x-x_0)(y-y_0)}\gamma (y),$$

$$D_{ij}$$

$${\partial }_x\alpha (\kappa |\rho (x)-\rho (y)|)=\kappa {\alpha }^{\prime }(\kappa |\rho (x)-\rho (y)|)\frac{({\rho }_1(x)-{\rho }_1(y)){\rho }_1^{\prime }(x)+({\rho }_2(x)-{\rho }_2(y)){\rho }_2^{\prime }(x)}{|\rho (x)-\rho (y)|},$$

$${\partial }_{xy}\alpha (\kappa |\rho (x)-\rho (y)|)$$

$$|{\partial }_{xy}\alpha (\kappa |\rho (\bar{x})-\rho (\bar{y})|)|\le O({\kappa }^2{\alpha }^{\prime \prime }(\kappa |\rho (\bar{x})-\rho (\bar{y})|)+\frac{Nk}{M}{\alpha }^{\prime }(k|\rho (\bar{x})-\rho (\bar{y})|))),$$

$$\underset{x\to \infty }{\lim }{\alpha }^{\prime }(x)=1.$$

$$\underset{x\in [a_1,\infty )}{\max }|{\alpha }^{\prime }(x)|\le b_1,$$

$$|\frac{N\kappa }{M}{\alpha }^{\prime }(\kappa |\rho (\bar{x})-\rho (\bar{y})|))|\le O(\frac{N\kappa }{M}).$$

$${\alpha }^{\prime \prime }(x)\sim \frac{-1}{4x^3}+\frac{25}{32x^5}+\frac{-3219}{512x^7}+\dots$$

$$|{\kappa }^2{\alpha }^{\prime \prime }(\kappa |\rho (\bar{x})-\rho (\bar{y})|)|=O(\frac{1}{\kappa |\rho (\bar{x})-\rho (\bar{y}){|}^3})=O(\frac{N^2}{M^3}),$$

$$|(x-x_0)(y-y_0)|=O(\frac{d_x{d}_y}{N^2})\le O(\frac{M}{N^2}).$$

$$|\omega (x-x_0)(y-y_0)|$$

$$K\approx U^L{U}^{L-1}\cdots U^h{S}^h{V}^h\cdots V^{L-1}{V}^L,$$

Full text

Keywords. Preconditioned iterative solver, interpolative decomposition butterfly factorization, LU factorization, electric-field integral equation (EFIE), scattering.

AMS subject classifications: 44A55, 65R10 and 65T50.

1 Introduction

Iterative and direct surface integral equation (IE) techniques are popular tools for the scattering analysis involving electrically large perfect electrically conducting (PEC) objects. In iterative techniques, fast matvec algorithms including multilevel fast multipole algorithms (MLFMA) [1, 2, 3, 4, 5, 6, 7], directional compression algorithms [8, 9, 10, 11], special function transforms [12, 13, 14, 15, 16], and butterfly factorizations (BFs) [17, 18, 19, 20, 21, 22] can be used to rapidly apply discretized IE operators to trial solution vectors. These techniques typically require $O({N_{i}^{\alpha}}N\log^{\beta}N)$ (${\alpha,}~{}\beta=1$ or $2$) CPU and memory resources, where $N$ is the dimension of the discretized IE operator and ${N_{i}}$ is the number of iterations required for convergence. The widespread adoption and success of fast iterative methods for solving real-world electromagnetic scattering problems can be attributed wholesale to their low computational costs when ${N_{i}}$ is small. However, in many applications when scatterers support resonances or are discretized using multi-scale/dense meshes, the corresponding linear system becomes ill-conditioned especially for first-kind integral equation formulations and ${N_{i}}$ can be prohibitively large. This motivates a significant amount of research devoted to the design of efficient semi-analytic or analytic preconditioners for integral equations for both open and closed surfaces [23, 24, 25, 25, 26].

Direct solvers do not suffer (to the same degree) from the aforementioned drawbacks, since they construct a compressed representation of the inverse of the discretized IE operator directly. Most state-of-the-art direct methods apply low rank (LR) approximations to compress judiciously selected off-diagonal blocks of the discretized IE operator and its inverse [27, 28, 29, 30, 31, 32, 33]. LR compression schemes are highly efficient and lead to nearly linear scaling direct solvers for electrically small [28, 34], or specially-shaped structures [35, 36, 37, 38]. However, for electrically large and arbitrarily shaped scatterers, the numerical rank of blocks of the discretized IE operators and its inverse is no longer small. As a result, there is no theoretical guarantee of low computational costs for LR schemes in this high-frequency regime; experimentally their CPU and memory requirements have been found to scale as $O(N^{\alpha}\log^{\beta}N)$ ($\alpha\in[2.0,3.0]$, $\beta\geq 1$) and $O(N^{\alpha}\log N)$ ($\alpha\in[1.3,2.0]$), respectively. More recently, butterfly factorizations have been applied to construct reduced-complexity direct solvers in the high-frequency regime [39, 40, 41]. These solvers construct butterfly factorizations with constant ranks for blocks (that are LR less compressible) in the discretized IE operator and its inverse and rely on fast randomized algorithms to construct the inverse in $O(N^{1.5}\log N)$ operations.

The lack of quasi-linear complexity direct solvers in the high-frequency regime motivates this work to develop a hierarchical butterfly-compressed algebraic preconditioner for solving EFIEs with $O(N\log^{2}N)$ CPU and memory complexity per iteration and up to $O(\log^{2}N)$ total iterations, for analyzing scattering from electrically large two-dimensional PEC objects with open surfaces. First, the interpolative decomposition butterfly factorization (IDBF) algorithm [20] is used to compress off-diagonal blocks of the discretized IE operator, leading to $O(N\log^{2}N)$ construction and application algorithms with the hierarchical IDBF (H-IDBF) of the IE operator. Second, we construct an approximate butterfly-compressed LU factorization of H-IDBF inspired by the work in [40]. In contrast to [40] that computes the LU factorization via randomized butterfly algebra, here the lower and upper triangular parts of the H-IDBF can directly serve as its approximate LU factorization. This is justified by the observation that the discretized IE operator and its LU factors exhibit similar oscillation patterns after a proper row-wise/column-wise ordering. This approximate LU factorization permits construction of a quasi-linear complexity preconditioner for EFIEs when applied to a wide range of scatterers. Compared to the analytic preconditioners in [25, 26], the proposed H-IDBF-LU preconditioner is capable of analyzing scattering from electrically large objects involving up to $100$ million unknowns.

2 Interpolative Decomposition Butterfly Factorization (IDBF)

Since the IDBF will be applied repeatedly in this paper, we briefly sumarize the $O(N\log N)$ IDBF algorithm proposed in [20] for a complementary low-rank matrix $K\in\mathbb{C}^{M\times N}$ with $M\approx N$.

Let $X$ and $\Omega$ be the row and column index sets of $K$. Two trees $T_{X}$ and $T_{\Omega}$ of the same depth $L=O(\log N)$, associated with $X$ and $\Omega$ respectively, are constructed by dyadic partitioning with approximately equal node sizes with leaf node sizes no larger than $n_{0}$. Denote the root level of the tree as level [math] and the leaf one as level $L$. Such a matrix $K$ of size $M\times N$ is said to satisfy the complementary low-rank property if for any level $\ell$, any node $A$ in $T_{X}$ at level $\ell$, and any node $B$ in $T_{\Omega}$ at level $L-\ell$, the submatrix $K_{A,B}$, obtained by restricting $K$ to the rows indexed by the points in $A$ and the columns indexed by the points in $B$, is numerically low-rank. See Figure 9 for an illustration of low-rank submatrices in a complementary low-rank matrix of size $16n_{0}\times 16n_{0}$.

Given $K$, or equivalently an $O(1)$ algorithm to evaluate an arbitrary entry of $K$, IDBF aims at constructing a data-sparse representation of $K$ using the ID of low-rank submatrices in the complementary low-rank structure (see Figure 9) in the following form:

[TABLE]

where the depth $L=O(\log N)$ is assumed to be even without loss of generality, $h=L/2$ is a middle level index, and all factors are sparse matrices with $O(N)$ nonzero entries. Storing and applying IDBF requires only $O(N\log N)$ memory and time.

The construction of the butterfly factorization was usually expensive though the application is cheap. The IDBF in [20] applies a linear scaling interpolation decomposition technique with carefully selected matrix skeleton of $K$ to achieve $O(N\log(N))$ factorization time, making the butterfly factorization more attractive in large-scale applications. Hence, the H-IDBF-LU to be introduced in the next section also admits nearly linear scaling factorization time.

3 H-IDBF-LU

3.1 Overview of the H-IDBF-LU preconditioner

Here we briefly describe the motivation and workflow of H-IDBF-LU in the context of EFIE for 2D $\mathrm{TM_{z}}$ scattering. Let $S$ denote a PEC cylindrical shell residing in free space. A current $J_{z}(\rho)$ is induced on $S$ by an incident electric field $E^{inc}_{z}(\rho)$. Enforcing a vanishing total electric field on $S$ results in the following EFIE:

[TABLE]

Here, ${\kappa}=2\pi/\lambda_{0}$ is the wavenumber, $\lambda_{0}$ represents the free-space wavelength, $\eta_{0}$ is the intrinsic impedance of free space, and $H^{(2)}_{0}$ denotes the zeroth-order Hankel function of the second kind. Upon discretization of the current with pulse basis functions and point matching the above EFIE, the following system of linear equations is attained [17]

[TABLE]

where $A$ is the impedance matrix with

[TABLE]

$b_{i}=E^{inc}_{z}(\rho_{i})$ and $x_{i}=J_{z}(\rho_{i})$ where $e\approx 2.718$, $\gamma\approx 1.781$, $w_{i}$ and $\rho_{i}$ are the length and center of the line segment associated with $i$-th pulse basis function. In the proposed solver, the linear system in (2) is rescaled to $aAx=ab$, where we choose $a:=\frac{1}{\max_{i}\{|A(i,1)|\}}$ or $a:=\left\lVert A\right\rVert_{2}$. Here $\left\lVert A\right\rVert_{2}$ can be rapidly computed via randomized SVD algorithms. Rescaling is a key step for the proposed LU preconditioner to avoid arithmetic overflow in solution vectors as the iteration count increases. Note that for second-kind integral equations the rescaling may not be necessary. With the abuse of notations, we will still use $Ax=b$ to denote the linear system after rescaling.

It was validated in [17, 39] that $A$ admits a hierarchical complementary low-rank property, i.e., off-diagonal blocks (only well-separated ones if necessary) are complementary low-rank matrices (See Figure 2 (left) for an example). [17] originally proposed a hierarchical compression technique with $O(N\log^{2}N)$ operations to compress the impedance matrix $A$; the IDBF in [20] proposed a more stable algorithm with the same accuracy to compress $A$. After compression, it requires $O(N\log^{2}N)$ operations to apply the impedance matrix and makes it possible to design efficient linear solvers for (2). The algorithm in [20] for compression of the impedance matrix is referred to as hierarchical IDBF (H-IDBF) in this paper.

Since the H-IDBF of the impedance matrix $A$ leads to an $O(N\log^{2}N)$ fast matvec algorithm in a purely algebraic fashion, an immediate question is how to construct a direct solver or an efficient preconditioner for the linear system in (2). It was observed (without proof) in [40] that the $L$ and $U$ factors in the LU factorization of $A$ are also hierarchically complementary low-rank matrices. Hence, [40] proposed an $O(N^{1.5}\log N)$ factorization algorithm based on randomized butterfly algebra to construct hierarchical compressions of the $L$ and $U$ factors, followed by an $O(N\log^{2}N)$ block triangular solution algorithm for rapidly application of the matrix inverse to a vector. The lack of a quasi-linear LU factorization algorithm as a direct solver in [40] motivates this work to construct an approximate LU factorization of $A$ in $O(N\log^{2}N)$ operations as an efficient preconditioner.

The main observation of the proposed approximate LU factorization of $A$ is that, the lower and upper triangular parts of $A$ have similarly oscillation patterns as the $L$ and $U$ factors of its LU factorization, when the scatterer consists of open surfaces that do not support high-Q resonance. It is also critical that the rows and columns of $A$ are properly ordered such that $A$ exhibits only $O(1)$ discontinuities with respective to $N$ in each row and column. Figure 3 and 4 show two examples of such scatterers; one is a spiral scatterer and the other one consists of two parallel strips. To be precise, define triangular matrices $\tilde{L}$ and $\tilde{U}$ as follows:

[TABLE]

and

[TABLE]

Let $L$ and $U$ be the lower and upper triangular matrices of the LU factorization of $A$, we observe $L\approx\tilde{L}$, $U\approx\tilde{U}$, and $\tilde{L}^{-1}A\tilde{U}^{-1}\approx I$ numerically, where $I$ is an identity matrix. See Figure 3 and 4 for the visualization of the approximation $\tilde{L}^{-1}A\tilde{U}^{-1}\approx I$ and the eigenvalue distribution of $\tilde{L}^{-1}A\tilde{U}^{-1}\approx I$. Take a spiral scatterer for example, the real (and imaginary) parts of $A$ (Figure 3(d)) and its LU factors (Figure 3(e)) exhibit very similar oscillation patterns. Before preconditioning, there are eigenvalues of $A$ clustered at the origin (Figure 3(b)), while the preconditioned system $\tilde{L}^{-1}A\tilde{U}^{-1}$ has no eigenvalues near the origin (Figure 3(c)). In fact, $\tilde{L}^{-1}A\tilde{U}^{-1}$ exhibits vanishingly small off-diagonal values (Figure 3(f)). The H-IDBFs of $\tilde{L}$ and $\tilde{U}$ are referred to as the H-IDBF-LU of the impedance matrix $A$ in this paper.

Note that $\tilde{L}$ and $\tilde{U}$ are the lower and upper triangular parts of $A$, whose H-IDBFs are readily available from the H-IDBF of $A$ without any extra cost. Following the block triangular solve algorithm for H-matrices [42], it only requires $O(N\log^{2}N)$ operations to apply $\tilde{L}^{-1}$ and $\tilde{U}^{-1}$ to a vector. Hence, the approximate LU factorization results in an $O(N\log^{2}N)$ preconditioner for the linear system in (2). To be more specific, instead of directly solving (2) using traditional iterative solvers, we solve a preconditioned linear system

[TABLE]

and

[TABLE]

Since the eigenvalues of $\tilde{L}^{-1}A\tilde{U}^{-1}$ are well separated from [math] and gathered around $1$ (see Figure 3 and 4 for two examples), the condition number of $\tilde{L}^{-1}A\tilde{U}^{-1}$ is well controlled and the number of iterations for solving (5) is small. In all of our examples, the iteration number is $O(1)$, $O(\log N)$, or $O(\log^{2}N)$ for scatterers consisting of open surfaces. As the construction, application and triangular solve algorithms all require $O(N\log^{2}N)$ operations, the overall preconditioner is also quasi-linear.

3.2 ** Algorithm description **

The remaining task in this section is to introduce the fast construction, application and solution phases of the H-IDBF-LU solver.

3.2.1 Construction and application of H-IDBF

As we have seen in the previous subsection, the construction of the proposed preconditioner, i.e., the H-IDBFs of $\tilde{L}$ and $\tilde{U}$, is an immediate result of the H-IDBF of $A$. Hence, it is sufficient to introduce the construction and application of the H-IDBF of $A$. In fact, these resemble standard techniques for hierarchical matrices [43, 44, 45], and can be summarized in Algorithm 1 and 2, respectively. For illustration purpose, we assume $N$ is a power of $2$ and omit the parameters of IDBF in all algorithms. Furthermore, we assume all off-diagonal blocks of $A$ are butterfly compressible (i.e., weak admissibility). However, the proposed solver trivially extends to arbitrary problem sizes and strong admissibility. Let $n_{0}$ denote the predefined leaf-level size parameter for all IDBFs, algorithm 1 constructs a IDBF of $L=\log_{2}N-\log_{2}n_{0}-v$ levels for each off-diagonal block of dimension $N/2^{v}$. Once constructed, algorithm 2 can apply the compressed $A$ to arbitrary vectors by accumulating the matvec results of each submatrix.

3.2.2 Solution of H-IDBF-LU

The previous subsection has introduced $O(N\log^{2}N)$ algorithms for constructing and applying the H-IDBF of $A$ and forming $\tilde{L}$ and $\tilde{U}$ as a preconditioner. Next, we will introduce $O(N\log^{2}N)$ algorithms for applying the inverse of $\tilde{L}$ and $\tilde{U}$ to a vector in this subsection, i.e, solving triangular systems in a format of H-IDBFs. Similarly to solving a triangular system in a format of H-matrix, the triangular solver for H-IDBF can also be done recursively as in Algorithm 3 and 4.

3.3 Algorithm Analysis

Here, we briefly analyze the effectiveness of the proposed LU preconditioner and the compressibility of off-diagonal blocks in the impedance matrix and its LU factors.

3.3.1 Oscillation patterns of $A$ and its LU factors

We observed that the oscillation pattern of $A$ resembles that of its LU factors. More specifically, $\tilde{L}\approx c_{l}L$ and $\tilde{U}\approx c_{u}U$, where $\tilde{L}$/$\tilde{U}$ are introduced in (3)/(4), $L$ and $U$ are the LU factors of $A$, and $c_{l}$ and $c_{u}$ are scalar constants with $c_{l}c_{u}=1$. In this subsection, we provide some physical explanations of this observation assuming that 1) the scatterer consists of smooth non-resonant surfaces; 2) the LU factorization follows a geometrical elimination order defined as follows.

Assume that the scatterers consist of one simple contour with an analytic boundary. Let $\rho:[0,1)\rightarrow\mathbb{R}^{2}$ be the analytic map that maps $[0,1)$ on to the boundary and discretize $[0,1)$ with a uniform grid $x_{i}=\frac{i-1}{N}$ for $i=1,\dots,N$. The geometrical order is the order that corresponds to the natural order of the discretized parameter $x_{i}$. When the scatterers consist of multiple simple contours, the ordering between contours is not studied here and the optimal ordering remains an open question.

For simplicity, it is assumed that $A$ has been rescaled such that $A(i,i)\approx 1$. We want to show $\tilde{L}\approx L$ and $\tilde{U}\approx U$, respectively. Note that the LU factorization process at the step ${p}\in\{1,...,N\}$ can be rewritten as

[TABLE]

where $A_{11}=L_{11}U_{11}$ represents the completed LU factorization of the leading submatrix of sizes ${p\times p}$, and $S_{22}=A_{22}-A_{21}U_{11}^{-1}L_{11}^{-1}A_{12}$ denotes the updated Schur complement to be factorized in the next step. It suffices to show $L_{11}^{-1}A_{12}\approx A_{12}$, $A_{21}U_{11}^{-1}\approx A_{21}$, and $S_{22}\approx A_{22}$. As the factorization follows the geometrical order, let $c_{1}$ and $c_{2}$ denote the subscatterer corresponding to the eliminated (i.e., $A_{11}$) and remaining (i.e., $A_{22}$) rows. See the curves highlighted in red and blue in Figure 5 (a)-(c) for illustrations.

First note that $S_{22}$ represents the numerical Green’s function that accounts for the curve $c_{1}$ as the background medium. As $c_{1}$ is smooth and $c_{2}$ resides on one side of $c_{1}$, the scattering component $A_{21}U_{11}^{-1}L_{11}^{-1}A_{12}$ does not dominate over the free-space interaction component $A_{22}$. Physically, this means that the current induced on $c_{1}$ by sources on $c_{2}$, namely $U_{11}^{-1}L_{11}^{-1}A_{12}$, has weak Green’s function interaction in the directions towards $c_{2}$. Therefore we can use an aggressive approximation $S_{22}\approx A_{22}$.

We then show $L_{11}^{-1}A_{12}\approx A_{12}$ and $A_{21}U_{11}^{-1}\approx A_{21}$ by induction. Note that ${p}=1$ is trivial and assume that this holds true for step ${p}-1$, i.e., $L_{11}\approx\tilde{L}_{11}$ and $U_{11}\approx\tilde{U}_{11}$ at step ${p}$. Let $b_{j}$ denote the $j$th column of $A_{12}$, $j=1,...,N-{p}$, the vector $v=A_{11}^{-1}b_{j}$ is the current solution on $c_{1}$ residing in free space excited by a dipole source residing on $c_{2}$. Note that $c_{2}$ is located on one side of $c_{1}$. For the linear system $A_{11}v=b_{j}$, we can define the forward operator $\tilde{U}_{11}$ and backward operator $\tilde{L}_{11}$ which only propagates the current away from and towards the source point, respectively. As $c_{1}$ is smooth and non-resonant, the linear system can be approximated as $\tilde{U}_{11}v\approx b_{j}$ by replacing the integral operator $A_{11}$ by the forward operator $\tilde{U}_{11}$ which only propagates the current away from the source point. For example, using the forward operator, the current $v(i_{1})$ does not contribute to the field $b_{j}(i_{2})$ given that $i_{1}<i_{2}$ (i.e., $x_{i_{1}}$ represents a point that are further away from the source when compared to $x_{i_{2}}$). Note that the approximation $A_{11}v\approx\tilde{U}_{11}v$ does not guarantee that $A_{11}\approx\tilde{U}_{11}$. For more detailed explanation of this approximation, the readers are referred to Section III of [46]. From such approximation, it is easy to see that $b_{j}=A_{11}v=L_{11}U_{11}v\approx\tilde{L}_{11}\tilde{U}_{11}v\approx\tilde{L}_{11}b_{j}$. Note that the first approximate equality holds via induction. In other words, $\tilde{L}_{11}^{-1}b_{j}\approx b_{j}$ or equivalently ${L}_{11}^{-1}A_{12}\approx\tilde{L}_{11}^{-1}A_{12}\approx A_{12}$. Similarly, it can be shown that $A_{21}U_{11}^{-1}\approx A_{21}$.

The analysis above explains the observation $\tilde{L}\approx L$, $\tilde{U}\approx U$ and suggests that $\tilde{L}\tilde{U}$ can be used as an efficient preconditioner. We would like to highlight a few related works here. Our analysis is closely related to the fix-point iteration methods such as the forward-backward method [46, 47] or the symmetrical Gauss-Seidel iteration method [48] for solving integral equations with a certain type of excitation vectors. When these methods are used as a preconditioner in a Krylov subspace iterative solver for arbitrary excitations, our physical explanations above show effectiveness of such preconditioner. Another related work is the sweeping preconditioner for differential equation solvers [49] which requires a geometrical elimination order such that the Schur complement represents well-compressible numerical Green’s functions. Here we show that for surface integral equation solvers, a geometrical order can induce Schur complement that resembles the free-space Green’s function for smooth, non-resonant scatterers. Both algorithms rely on the fact that the geometry points corresponding to the remaining unknowns reside approximately on one side of the geometry points corresponding to the the already eliminated unknowns.

3.3.2 Butterfly compressibility of off-diagonal blocks

Here we analyze the complementary low-rank property of the off-diagonal blocks of $A$, and it follows from the arguments in Section 3.3.1 that the LU factors of $A$ are also butterfly compressible. Note that the compressibility of off-diagonal blocks resulting from strong-admissibility is well-studied in [9, 17], here we extend the compressibility analysis to the off-diagonal blocks from weak-admissibility.

To simplify the analysis, we assume that the scatterer has an analytic boundary and its diameter and length are of order $1$. The parametrization of the boundary leads to an impedance matrix of size $N$ by $N$ and the frequency parameter ${\kappa}$ is of order $N$. Let $\rho:[0,1)\rightarrow\mathbb{R}^{2}$ be the analytic map that maps $[0,1)$ on to the boundary of the scatterer and discretize $[0,1)$ with a uniform grid $x_{i}=\frac{i-1}{N}$ for $i=1,\dots,N$. Let $\rho_{i}=\rho(x_{i})$, then the impedance matrix $A$ is

[TABLE]

It is sufficient to show the butterfly compressibility of the off-diagonal block for $1\leq j<\frac{N}{2}<i\leq N$ assuming $N$ is an even number. The analysis for other off-diagonal blocks is similar. By the theory of non-oscillatory phase functions for special functions [50, 51, 52], the Hankel function $H_{0}^{(2)}(x)$ for $x\in(0,\infty)$ admits a representation

[TABLE]

where $\alpha(x)$ is an analytic non-oscillatory phase function. For more details about the definition, computation, and asymptotic expansion of the non-oscillatory phase function, the reader is referred to Section 2.5 of [52], where the notation $\alpha(t)$ is used to represent the phase function as well. If we introduce a non-oscillatory matrix $B$ and a purely oscillatory matrix $C$ as follows,

[TABLE]

and

[TABLE]

then $A$ is the Hadamard product of $B$ and $C$.

Note that ${\kappa}|\rho_{i}-\rho_{j}|\geq O(1)$ for all $i\neq j$. Hence, $M(x_{i},x_{j}):=\sqrt{\frac{2}{\pi{\kappa}|\rho(x_{i})-\rho(x_{j})|}}\frac{1}{\sqrt{\alpha^{\prime}({\kappa}|\rho(x_{i})-\rho(x_{j})|)}}$ is analytic and non-oscillatory in $(x_{i},x_{j})$ leading to the following lemma.

Lemma 3.1**.**

Suppose the map $\rho$ associated to the boundary of the scatterer is analytic and nearly isomatric, i.e., $c_{1}|x-y|\leq|\rho(x)-\rho(y)|\leq c_{2}|x-y|$ for some positive constants $c_{1}$ and $c_{2}$. Suppose $B$ of size $N$ by $N$ is defined in (8), then the rank of the off-diagonal block of $B$ for $1\leq j<\frac{N}{2}<i\leq N$ is $O(\log(N))$ with a prefactor depending on a relative approximation error $\epsilon$ and independent of $N$.

Proof.

Consider the range $x\in[x_{0}-\frac{d_{x}}{2},x_{0}+\frac{d_{x}}{2}]$ and $y\in[y_{0}-\frac{d_{y}}{2},y_{0}+\frac{d_{y}}{2}]$ with $(x_{0}-\frac{d_{x}}{2})-(y_{0}+\frac{d_{y}}{2})\geq\frac{1}{N}$. Then as discussed previously, $M(x,y)$ is analytic and non-oscillatory in $(x,y)$. Hence, as long as $d_{x}$ and $d_{y}$ are sufficiently small depending on $M(x,y)$ and are still of order $1$, a low-rank separation approximation for $M(x,y)$ with a relative error $\epsilon$ can be constructed via Taylor expansion or Chebyshev interpolation. The number of terms in the low-rank separation approximation depends only on $\epsilon$ and is independent of $N$. The conclusion is an immediate result of this fact.

For example, let us assume low-rank approximations are valid as long as $d_{x}\leq\frac{1}{4}$ and $d_{y}\leq\frac{1}{4}$. Then the off-diagonal block of $B$ is partitioned hierarchically as in Figure 6 (left), the submatries of which have rank $O(1)$ independent of $N$. For an off-diagonal matrix of size $O(N)$ by $O(N)$, there are $O(\log(N))$ levels of submatrices. Hence, the total rank of the whole off-diagonal matrix in Figure 6 (left) is $O(\log(N))$. ∎

The next lemma below analyzes the low-rank structure of the purely oscillatory matrix $C$.

Lemma 3.2**.**

Suppose the map $\rho$ associated to the boundary of the scatterer is analytic and nearly isomatric, i.e., $c_{1}|x-y|\leq|\rho(x)-\rho(y)|\leq c_{2}|x-y|$ for some positive constants $c_{1}$ and $c_{2}$. Consider an off-diagonal contiguous submatrix of $C$ in (9) of size $d_{x}$ by $d_{y}$ with $d_{x}d_{y}\leq M$ and the submatrix is at least $O(M)$ entries away from the diagonal of $C$. Then the submatrix is rank $O(1)$ with a prefactor depending on a relative approximation error $\epsilon$ and independent of $N$, when $M$ is larger than an $O(1)$ constant independent of $N$.

Proof.

Without loss of generality, we consider

[TABLE]

with

[TABLE]

Recall that $x$ and $y$ in $[0,1]$ are discretized with $N$ uniform grid points to form the matrix $C$. Hence, any submatrix in Lemma 3.2 is from the discretization of the conditions above. See Figure 6 (middle) for an illustration of the submatrix of $C$ corresponding to the range in (10).

By the Taylor expansion of $\alpha({\kappa}|\rho(x)-\rho(y)|)$ at the point $(x_{0},y_{0})$, there exists a function $\beta(x)$ and a function $\gamma(y)$ such that

[TABLE]

with $\bar{x}\in[x_{0}-\frac{d_{x}}{2N},x_{0}+\frac{d_{x}}{2N}]$ and $\bar{y}\in[y_{0}-\frac{d_{y}}{2N},y_{0}+\frac{d_{y}}{2N}]$. Hence, the separability of the function $e^{-i\alpha({\kappa}|\rho(x)-\rho(y)|)}$ is equivalent to the separability of $e^{-i\partial_{xy}\alpha({\kappa}|\rho(\bar{x})-\rho(\bar{y})|)(x-x_{0})(y-y_{0})}$, i.e., to show Lemma 3.2 for the matrix $C$, it is sufficient to show the same property for the matrix $D$ defined below.

[TABLE]

where $\omega=\partial_{xy}\alpha({\kappa}|\rho(\bar{x})-\rho(\bar{y})|)$ and $x_{i}=y_{i}=\frac{i-1}{N}$ for $1\leq i\leq N$.

Let $\rho(x)=[\rho_{1}(x),\rho_{2}(x)]$. Then

[TABLE]

and

[TABLE]

Hence, for the range considered in (10), we have $|\bar{x}-\bar{y}|\geq O(\frac{M}{N})$, then

[TABLE]

since $\rho$ is nearly isometric. Note that $\alpha(x)$ is non-oscillatory and analytic with the property (see Equations (57) in [52]) that

[TABLE]

Hence, there exists a positive constant $b_{1}$ independent of $N$ such that

[TABLE]

where $a_{1}$ is a positive number such that $a_{1}\leq{\kappa}|\rho(x_{i})-\rho(x_{j})|$ as long as $i\neq j$, where $x_{i}=\frac{i-1}{N}$ for $i=1,\dots,N$, i.e., $a_{1}=O(1)$. Hence,

[TABLE]

By the asymptotics of Bessel functions in Equations (63) in Section 2.5 of [52], the following expansion is valid for all $x>0$:

[TABLE]

By a finite term truncation of the above expansion, there exist positive constants $c_{2}$ and $b_{2}$ such that $|\alpha^{\prime\prime}(x)|\leq\frac{c_{2}}{4x^{3}}$ for $x\in[b_{2},\infty)$. Hence, as long as $M$ is sufficiently large, e.g, larger than $O(b_{2})=O(1)$, ${\kappa}|\rho(\bar{x})-\rho(\bar{y})|\geq b_{2}$, and hence

[TABLE]

since ${\kappa}=O(N)$ and $|\rho(\bar{x})-\rho(\bar{y})|=O(|\bar{x}-\bar{y}|)\geq O(\frac{M}{N})$. Note that

[TABLE]

Hence,

[TABLE]

which means that the submatrix of $D$ corresponding to the range in (10) is non-oscillatory and hence is rank $O(1)$ with a prefactor depending on a relative approximation error $\epsilon$ and independent of $N$. By (11), this conclusion is also true for the matrix $C$ and hence, we have completed the proof of Lemma 3.2.

∎

An immediate result of Lemma 3.2 is that the off-diagonal block of $C$ is a hierarchically complementary low-rank matrix, e.g., see Figure 6 (right) for an illustration. This conclusion is also true for the impedance matrix $A$. Theoretically, the fast matvec of $A$ can be performed via hierarchically decomposition $A$ as in Figure 6 (right) and applying the IDBF for each submatrix. This leads to a fast matvec with a complexity $O(N\log^{2}(N))$. Numerically, a more convenient way is to directly compress the whole off-diagonal block with IDBF, which also leads to a fast matvec with a numerical complexity $O(N\log^{2}(N))$.

4 Numerical results

This section demonstrates the efficiency and accuracy of the proposed preconditioner via its application to wide classes of open surfaces: a semicircle, a corrugated corner reflector, a spiral line, two parallel strips, a cup-shaped cavity, and an array of open arcs. In all examples, the surfaces are discretized with approximately 20 pulse basis functions per wavelength. The leaf sizes in the dyadic trees of each butterfly-compressed block are set to approximately $n_{0}=200$. The compression tolerance in linear scaling ID is set to $\epsilon$=1.0E-4, and the oversampling parameter $t$ in the linear scaling ID is set to $1$. For the aforementioned surfaces, the maximum butterfly ranks k among all blocks of the impedance matrices are respectively 7, 11, 10, 7, 14, 13, which are almost independent of the matrix sizes $N$. The matrix entries are scaled with a scalar such that the largest diagonal entry has unit magnitude. A transpose-free quasi-minimal residual (TFQMR) iterative solver is utilized with a convergence tolerance 1.0E-5. All experiments are performed on the Cori Haswell machine at NERSC, which is a Cray XC40 system and consists of 2388 dual-socket nodes with Intel Xeon E5-2698v3 processors running 16 cores per socket. The nodes are configured with 128 GB of DDR4 memory at 2133 MHz.

First, the accuracy of the proposed preconditioner is demonstrated by changing the matrix size $N$ from 5000 to 5,000,000. Let $x^{t}$ be a randomly generated $N\times 1$ vector, the right hand side is generated as $b=Ax^{t}$. Note that $A$ is replaced by its H-IDBF compression $F$ when $N>10000$. Let $x^{a}$ denote the solution vector computed from the preconditioned linear system. The solution error is defined as $\epsilon^{a}=\left\lVert x^{a}-x^{t}\right\rVert_{2}/\left\lVert x^{t}\right\rVert_{2}$. From Table 1, reasonably good accuracy has been observed for all classes of surfaces. However, a slight degradation in accuracy when $N$ increases is also observed.

Next, the computational efficiency is demonstrated by investigating the computation time and memory as matrix size $N$ increases. The construction time, iterative solution time, and overall memory are plotted in Figure 7 for all test surfaces. It can be easily validated that all three quantities scale as $O(N\log^{2}N)$ as predicted. It is worth mentioning that for the semicircle example, the proposed preconditioner permits rapid solution for $N=100$ million, which is very competitive to the state-of-the-art MLFMA-based iterative solvers. We also observe that the solution time (i.e., application and triangular solve of H-IDBF-LU in the iterative solver) dominates the overall computation time especially for surfaces that requires constant but a high number of iterations.

Finally, we compare the iteration counts of the proposed preconditioner and an iterative solver without any preconditioner (see Figure 8). First, we investigate the dependence on $N$ for different surface shapes. The iteration counts required by the iterative solver without preconditioner grow rapidly for all surfaces (dashed lines in Figure 8(a)). In contrast, for surfaces including the semicircle, spiral lines, and corrugated corners, the iteration counts using the proposed preconditioner stay as small constant (typically less than 30); for the other surfaces, the observed iteration counts scale at most as $O(\log^{2}N)$. It’s worth mentioning that for complicated geometries such as highly resonant cavities or closed surfaces, the iteration counts can grow much faster. To see this, we then investigate the dependence on the degree of resonance or multiple reflections using an open arc with varying angles (Figure 8(b)), an Archimedean spiral with varying rotation angles (Figure 8(c)), (d) a 1-D array of open arcs with different element counts (Figure 8(d)). For the arc and spiral, the surface supports a higher degree of resonance as the angle increases (especially beyond $2\pi$), therefore the proposed preconditioner becomes less effective but still shows significant improvement compared to the iterative solver without any preconditioner. For the array of arcs, the scatterers support a dominant direction of propagation (along the direction of repetition) as the element increases, therefore the proposed preconditioner becomes very effective.

5 Conclusion and discussion

This paper has introduced a simple and efficient algorithm, the hierarchical interpolative decomposition butterfly-LU factorization (H-IDBF-LU) preconditioner for solving two-dimensional electric-field integral equations (EFIEs) in electromagnetic scattering problems of perfect electrically conducting objects with open surfaces. H-IDBF-LU consists of two main parts: the first part applies the newly developed interpolative decomposition butterfly factorization (IDBF) to compress dense blocks of the discretized EFIE operator to expedite its application; the second part treats the lower and upper triangular part of the IDBF as an approximate LU factorization of the EFIE operator leading to an efficient preconditioner in iterative solvers.

Both the memory requirement and computational cost of the H-IDBF-LU solver scale as $O(N\log^{2}N)$ in one iteration; the total number of iterations required for a reasonably good accuracy scales as $O(1)$ to $O(\log^{2}N)$ in all of our numerical tests. Our algorithm is simple to implement, automatically adapts to different structures with open surfaces, and is competitive with state-of-the-art MLFMA-based algorithms. A user-friendly MATLAB package, ButterflyLab (https://github.com/ButterflyLab/ButterflyLab), and a distributed parallel Fortran/C++ package, ButterflyPack (https://github.com/liuyangzhuan/ButterflyPACK), are freely available online.

The lower and upper triangular part of the EFIE operator can serve as an approximate LU factorization of the EFIE operator is an interesting observation and deserves much theoretical attention in the future. This idea could also be applied to three-dimensional EFIE’s with open surfaces and we will explore this in future work.

Acknowledgments. H. Yang was partially supported by Grant R-146-000-251-133 in the Department of Mathematics at the National University of Singapore, by the Ministry of Education in Singapore under the grant MOE2018-T2-2-147, and the start-up grant of the Department of Mathematics at Purdue University. Y. Liu thanks the support from the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program through the FASTMath Institute under Contract No. DE-AC02-05CH11231 at Lawrence Berkeley National Laboratory. The authors also thank Dr. Xiaoye Sherry Li for useful discussions and suggestions on the paper.

Appendix: IDBF

5.1 Overview

Since the IDBF will be applied repeatedly in this paper, we briefly review the $O(N\log N)$ IDBF algorithm proposed in [20] for a complementary low-rank matrix $K\in\mathbb{C}^{M\times N}$ with $M\approx N$ for the purpose of completeness. Let $X$ and $\Omega$ be the row and column index sets of $K$. Two trees $T_{X}$ and $T_{\Omega}$ of the same depth $L=O(\log N)$, associated with $X$ and $\Omega$ respectively, are constructed by dyadic partitioning with approximately equal node sizes with leaf node sizes no larger than $n_{0}$. Denote the root level of the tree as level [math] and the leaf one as level $L$. Such a matrix $K$ of size $M\times N$ is said to satisfy the complementary low-rank property if for any level $\ell$, any node $A$ in $T_{X}$ at level $\ell$, and any node $B$ in $T_{\Omega}$ at level $L-\ell$, the submatrix $K_{A,B}$, obtained by restricting $K$ to the rows indexed by the points in $A$ and the columns indexed by the points in $B$, is numerically low-rank. See Figure 9 for an illustration of low-rank submatrices in a complementary low-rank matrix of size $16n_{0}\times 16n_{0}$.

Given $K$, or equivalently an $O(1)$ algorithm to evaluate an arbitrary entry of $K$, IDBF aims at constructing a data-sparse representation of $K$ using the ID of low-rank submatrices in the complementary low-rank structure (see Figure 9) in the following form:

[TABLE]

where the depth $L=O(\log N)$ is assumed to be even without loss of generality, $h=L/2$ is a middle level index, and all factors are sparse matrices with $O(N)$ nonzero entries. Storing and applying IDBF requires only $O(N\log N)$ memory and time.

In what follows, uppercase letters will generally denote matrices, while the lowercase letters $c$, $p$, $q$, $r$, and $s$ denote ordered sets of indices. For a given index set $c$, its cardinality is written as $|c|$. Given a matrix $A$, $A_{pq}$, $A_{p,q}$, or $A(p,q)$ is the submatrix with rows and columns restricted to the index sets $p$ and $q$, respectively. We also use the notation $A_{:,q}$ to denote the submatrix with columns restricted to $q$. $s:t$ is an index set containing indices $\{s,s+1,s+2,\dots,t-1,t\}$. For the sake of simplicity, we assume that $N=2^{L}n_{0}$, where $n_{0}=O(1)$ is the number of column or row indices in a leaf in the dyadic trees of row and column spaces, i.e., $T_{X}$ and $T_{\Omega}$, respectively. In practical numerical implementations, this assumption is not required.

5.2 Linear scaling Interpolative Decompositions

Suppose $A$ has rank $k$, a rank revealing QR decomposition to $A_{s,:}$ gives

[TABLE]

where $s$ is an index set containing $tk$ carefully selected rows of $A$ with $t$ as an oversampling parameter, $Q\in\mathbb{C}^{t{k\times{k}}}$ is an orthogonal matrix, $R\in\mathbb{C}^{{k}\times n}$ is upper trapezoidal, and $\Lambda\in\mathbb{C}^{n\times n}$ is a carefully chosen permutation matrix such that $R_{1}\in\mathbb{C}^{{k\times k}}$ is nonsingular. These $t{k}$ rows can be chosen from the Mock-Chebyshev grids of the row indices as in [22, 53, 54]. Let

[TABLE]

then

[TABLE]

which is equivalent to the traditional form of a column ID,

[TABLE]

where $q$ is the complementary set of $p$, ∗ denotes the conjugate transpose of a matrix, and $V$ is the column interpolation matrix.. It can be easily shown that all the steps above require only $O({k}^{2}n)$ operations and $O({k}n)$ memory.

In practice, the true rank of $A$ is not available i.e., ${k}$ is unknown. As is in standard randomized algorithms, we could choose to fix a test rank ${k}\leq n$ or fix the approximation accuracy $\epsilon$ and find a numerical rank ${{k}_{\epsilon}}$ such that

[TABLE]

with $T\in\mathbb{C}^{{{k}_{\epsilon}\times(n-{k}_{\epsilon})}}$ and $V\in\mathbb{C}^{{{k}_{\epsilon}}\times n}$. We refer to this linear scaling column ID with an accuracy tolerance $\epsilon$ and a rank parameter ${k}$ as $(\epsilon,{k})$-cID ($(\epsilon,{k})$-cID for short). For convenience, we will drop the term $(\epsilon,{k})$ when it is not necessary to specify it.

Similarly, a row ID for the matrix $A\in\mathbb{C}^{m\times n}$

[TABLE]

can be attained by performing cID on $A^{*}$ with $O({k}^{2}m)$ operations and $O({k}m)$ memory. We refer to this linear scaling row ID as ${(\epsilon,{k})}$-rID and $U$ as the row interpolation matrix.

5.3 Leaf-root complementary skeletonization (LRCS)

Assume that at the leaf level of the row (and column) dyadic trees, the row index set $r$ (and the column index set $c$) of $A$ are divided into leaves $\{r_{i}\}_{1\leq i\leq m}$ (and $\{c_{i}\}_{1\leq i\leq m}$) in the following way:

[TABLE]

with $|r_{i}|=n_{0}$ (and $|c_{i}|=n_{0}$) for all $1\leq i\leq m$, where $m=2^{L}=\frac{N}{n_{0}}$, $L=\log_{2}N-\log_{2}n_{0}$, and $L+1$ is the depth of the dyadic trees $T_{X}$ (and $T_{\Omega}$). ${\it rID}$ is applied to each $A_{r_{i},:}$ to compute the row interpolation matrix and denote it as $U_{i}$; the associated skeleton indices is denoted as $\hat{r}_{i}\subset r_{i}$. Let $\hat{r}=[\hat{r}_{1},\hat{r}_{2},\cdots,\hat{r}_{m}]$, then $A_{\hat{r},:}$ is the important skeleton of $A$ and we can arrange all the small ID factors into a larger matrix factorization as follows:

[TABLE]

Similarly, ${\it cID}$ is applied to each $A_{\hat{r},c_{j}}$ to obtain the column interpolation matrix $V_{j}$ and the skeleton indices $\hat{c}_{j}\subset c_{j}$. Then finally we form the LRCS of $A$ as

[TABLE]

For a concrete example, Figure 10 visualizes the non-zero pattern of the LRCS in (21).

It is noteworthy that we only generate and store the skeleton of row and column index sets corresponding to $M$ and $S$, instead of computing $M$ and $S$ explicitly. Hence, it only takes $O(\frac{{k}^{3}}{n_{0}}N)$ operations and $O(\frac{{k}^{2}}{n_{0}}N)$ memory to generate and store the factorization in (21), since there are $2m=\frac{2N}{n_{0}}$ IDs in total.

5.4 Matrix splitting with complementary skeletonization (MSCS)

A complementary low-rank matrix $A$ (with row and column dyadic trees $T_{X}$ and $T_{\Omega}$ of depth $L+1$ and with $m=2^{L}$ leaves) can be split into a $2\times 2$ block matrix, $A=\begin{pmatrix}A_{11}&A_{12}\\ A_{21}&A_{22}\end{pmatrix}$, according to the nodes of the second level of the dyadic trees $T_{X}$ and $T_{\Omega}$ (right next to the root). As a result, each $A_{ij}$ is also a complementary low-rank matrix with its row and column dyadic trees of $L-1$ levels. For example, $A$ in Figure 9 is a five-level complementary low-rank matrix and $A_{11}$ is a three-level complementary low-rank matrix (see the highlighted submatrices in Figure 9).

Suppose $A_{ij}\approx U_{ij}S_{ij}V_{ij}$, for $i,j=1,2$, is the LRCS of $A_{ij}$. Then $A\approx USV$, where

[TABLE]

The factorization in (22) is referred as the MSCS. Recall that the middle factor $S$ is not explicitly computed, resulting in a linear scaling algorithm for forming (22). Figure 11 visualizes the MSCS of a complementary low-rank matrix $A$.

5.5 Recursive MSCS

Now we apply MSCS recursively to get the full IDBF. Suppose we have the first level of MSCS as $A\approx U^{L}S^{L}V^{L}$ with

[TABLE]

By construction, we know $S^{L}_{ij}$ are complementary low-rank. Next, apply MSCS to each $S^{L}_{ij}$:

[TABLE]

where

[TABLE]

Finally, organizing (24) forms $S^{L}\approx U^{L-1}S^{L-1}V^{L-1}$ (see its visualization in Figure 12), where

[TABLE]

[TABLE]

leading to a second level factorization of $A$ as $A\approx U^{L}U^{L-1}S^{L-1}V^{L-1}V^{L}$.

Similarly, we can apply MSCS recursively to each $S^{\ell}$ and assemble matrix factors hierarchically for $\ell=L$, $L-1$, $\dots$, $L/2$ to obtain

[TABLE]

where $h=L/2$. In the entire computation procedure, linear IDs only require $O(1)$ operations for each low-rank submatrix, and hence at most $O(N)$ for each level of factorization, and $O(N\log N)$ for the whole IDBF.