# Hierarchical Matrices Method and Its Application in Electromagnetic Integral Equations

Authors: Han Guo; Jun Hu; Hanru Shao; Zaiping Nie
Journal: International Journal of Antennas and Propagation
 (2012)

---

## Abstract

Hierarchical (ℋ-) matrices method is a general mathematical framework providing a highly compact representation and efficient numerical arithmetic. When applied in integral-equation- (IE-) based computational electromagnetics, ℋ-matrices can be regarded as a fast algorithm; therefore, both the CPU time and memory requirement are reduced significantly. Its kernel independent feature also makes it suitable for any kind of integral equation. To solve ℋ-matrices system, Krylov iteration methods can be employed with appropriate preconditioners, and direct solvers based on the hierarchical structure of ℋ-matrices are also available along with high efficiency and accuracy, which is a unique advantage compared to other fast algorithms. In this paper, a novel sparse approximate inverse (SAI) preconditioner in multilevel fashion is proposed to accelerate the convergence rate of Krylov iterations for solving ℋ-matrices system in electromagnetic applications, and a group of parallel fast direct solvers are developed for dealing with multiple right-hand-side cases. Finally, numerical experiments are given to demonstrate the advantages of the proposed multilevel preconditioner compared to conventional “single level” preconditioners and the practicability of the fast direct solvers for arbitrary complex structures.


---

## Body

## 1. Introduction

Integral equation (IE) method [1] is widely used in electromagnetic analysis and simulation. Compared to finite difference time domain (FDTD) method [2] and Finite element (FE) method [3], IE method is generally more accurate and avoids annoying numerical dispersion problems as well as complex boundary conditions. Despite the many advantages, IE method is often involved in enormous computational consumption, since its numerical discretization leads to dense linear system. With the efforts for years, various fast algorithms had developed focusing on reducing the computational complexity for IE method such as adaptive integral method (AIM) [4], fast multipole method (FMM) [5], IE-Fast Fourier transform (IE-FFT) [6], and fast low-rank compression methods [7], Though these fast algorithms are based on different theories, they use the same idea to reduce CPU time and memory usage complexity, which is to compute and store the major entries of the dense system matrix indirectly, and employ iterative methods instead of direct methods (such as LU decomposition) to solve the system. Preconditioners are often used to accelerate the convergence of iterative methods. However, iterative methods cannot always guarantee a reasonable solution with high precision; therefore, in some complex cases, we expect to employ powerful preconditioners to obtain visible acceleration of convergence or even use direct methods to avoid this problem entirely. But to implement these ideas in traditional fast algorithms encounters some difficulties that are, once a fast algorithm is applied, the major entries of the system matrix are computed and stored indirectly, and only a few entries can be accessed to form the sparse pattern which is essential to construct the preconditioners. This condition confines the flexibility of the preconditioning; therefore, a powerful preconditioner is generally hard to get. For direct methods, accessibleness of the whole matrix entries is also prerequisite, and employing them in fast algorithms is difficult as well as flexible preconditioning.Hierarchical (ℋ-) matrices method [8, 9] is a general mathematical framework providing a highly compact representation and efficient numerical arithmetic. Applying ℋ-matrices in IE method can reduce both CPU time and memory usage complexity significantly, so it can be regarded as a fast algorithm to IE method. Unlike the traditional fast algorithms mentioned above, any entry of ℋ-structured system matrix can be recovered easily, though the major entries are still implicitly expressed. This quality enables us to construct more efficient preconditioners, and moreover, a fast direct solver is available since its unique format. In this paper, a novel sparse approximate inverse (SAI) preconditioner in multilevel fashion is proposed to accelerate the convergence rate of Krylov iterations for solving ℋ-matrices system in electromagnetic applications, which is mentioned in [10], and a group of parallel fast direct solvers are developed for dealing with multiple right-hand-side cases.This paper is organized as follows. Section2 gives a brief review of the IE method and basic conception of ℋ-matrices. In Section 3, we elaborate the construction of the proposed multilevel SAI preconditioner and the implementation of parallel fast direct solvers. Numerical experiments are given in Section 4 to demonstrate the advantages of the proposed multilevel preconditioner compared to conventional “single level” preconditioners and the practicability of the fast direct solvers for arbitrary complex structures. Finally, some conclusions are given in Section 5.


## 2.ℋ-Matrices Representation for IE Method

We first proceed with a description of the IE method for solving electromagnetic scattering problems from 3D perfectly electric conductor (PEC). For concise introduction, only the electric field integral equation (EFIE) [1] is considered. The EFIE is written as(1)t̂⋅∫s′ds′G¯(r,r′)J(r′)=4πikηt̂⋅Ei(r).Discretize integral equation (1) by expanding J with local basis functions(2)J(r)=∑n=1Nxnf(r),
where N is the number of unknowns, denoting the vector basis functions, and xn is the unknown expansion coefficients. Applying Galerkin’s method results in a matrix equation(3)A⋅x=b,
where(4)Amn=∫sdsfm(r)⋅∫s′ds′G¯(r,r′)⋅fn(r′),bm=4πikη∫sdsfm(r)⋅Ei(r).
If matrix A is in ℋ-matrices formatted, we use Aℋ to denote it.Considering that the construction of proposed multilevel SAI preconditioner is built upon the structure ofℋ-matrices, some of its basic concepts need be reviewed at first. The basis functions set indexed by ℐ:={1,2,…,N} can be constructed as a tree Tℐ=(V,E) with vertex set V and edge set E. For each vertex v∈V, we define the set of son of v as Sℐ(v):={w∈V∣(v,w)∈E}. The tree Tℐ is called a cluster tree if the following conditions hold:(5)IistherootofTI,∀v∈V:ifSI(v)≠∅,wehavev=⋃w∈SI(v)w,∀w1,w2∈SI(v)withw1≠w2:w1∩w2=∅.Each vertexv in Tℐ is called a cluster, representing a bunch of basis functions. Typically, each vertex has two sons. If the amount of basis functions contained in a cluster is less than a certain number of nmin, that cluster has no sons, also called the leaves of Tℐ. By rearranging the basis functions by their indices, they are numbered consecutively in each cluster, and, moreover, they are concentrated geometrically. Precisely, let Ωv denote the domain of basis functions represented by cluster v, and it is bounded by the cardinality of v:(6)diam⁡(Ωv)≤cg|v|
in which diam(·) is the Euclidean diameter of a set. cg is a real constant; let the inequality valid for all clusters.The electromagnetic interaction of any two clusters, including self-interaction, maps certain subblock of the system coefficient matrix. Practically, most of these subblocks can be approximated by low-rank matrices with high-accuracy. Therefore, a systemic and appropriate partitioning procedure is needed for the coefficient matrix. Based on the cluster treeTℐ which contains hierarchy of partitions of ℐ, we are able to construct the block-cluster tree Tℐ×ℐ describing a hierarchical partition of ℐ×ℐ by the following maps:(7)SI×I(t×s)={∅,ifSI(t)=∅orSI(s)=∅,∅,ift×sisadmissible,SI(t)×SI(s),otherwise,
in which the definition of Sℐ×ℐ(·) is similar to Sℐ(·) that denotes the sons of certain block-cluster t×s. Admissible is a criterion which will be elaborated later and decides whether the subblock should be approximated by low-rank matrix, or split and combined by its sons, suspending for further transaction. Finally, the whole system matrix is segmented into pieces of subblocks by the procedure above. Corresponding to the block-cluster tree, these subblocks are called the leaves of Tℐ×ℐ, denoted by ℒ(Tℐ×ℐ).If two clusters are well separated geometrically, the Green function which connects the interaction of them is barely varying in their domains. That means, only a few patterns of interactional vectors can represent the whole mode; therefore, the subblock representing their interaction is rank deficit. In order to discerning these subblocks appropriately, we introduce the admissibility condition [9]:(8)max⁡{diam⁡(Ωt),diam⁡(Ωs)}≤ηdist⁡(Ωt,Ωs),
in which dist(·,·) is the Euclidean distance of two sets, Ωt and Ωs denote the domains of basis functions of cluster t, and s and η are a controllable real parameter, which is also illustrated in Figure 1. If the statement (8) is true, sub-block t×s is admissible, otherwise it is nonadmissible. The leaves of block-cluster tree ℒ(Tℐ×ℐ) are either admissible or nonadmissible. For admissible subblocks, we substitute them by low-rank representation through specific algorithms such as adaptive cross approximation (ACA) [11]. For nonadmissilbe subblocks, we calculate the entries and store them directly. The total computational complexity for constructing is close to O(NlogN), in which N regards to the number of unknowns, and we give a numerical investigation Part A of Section 4.Figure 1
Two basis function clusters domains and their distance, describe the definition ofadmissible condition.


## 3. Multilevel SAI Preconditioning and Fast Direct Methods

### 3.1. Multilevel SAI Preconditioning

If the Krylov iterative methods are used to solve the linear system, we always expect to find a high-performance preconditioner to accelerate the convergence. Generally, instead of solving the linear system of the formAℋx=b, we solve it by the following forms:
(9)MAHx=MborAHMy=b,x=My.
In which M is a sparse matrix satisfying the condition M≈Aℋ-1 to a certain degree. Therefore, MAℋ≈I or AℋM≈I can make the iterative solver more efficient. Because of limited space, we only discuss the left preconditioning case below. The right preconditioning case can be analyzed similarly. Traditional preconditioner has a fixed form that we could only execute “single level” preconditioning. Here, we elaborate how to construct a multilevel preconditioner under ℋ-matrices format.For a block-cluster treeTℐ×ℐ in ℋ-matrices, let L(Tℐ×ℐ) denotes its depth and ℒ(Tℐ×ℐ) denotes all the leaves of the cluster tree. Any leaf block-cluster b=(t×s)∈ℒ(Tℐ×ℐ) belongs to certain level of Tℐ×ℐ; hence, 1≤level(b)≤L(Tℐ×ℐ). We define that he finest level is the 1st level, where the smallest block cluster is located. A substructure of Aℋ denoted by Aℋ(l) can be defined as follows:(10)AH(l)(i,j)={AH(i,j),(i,j)∈B(l),0,otherwise
in which B(l) is a subset of ℒ(Tℐ×ℐ), (11)B(l)={b∣b∈L(TI×I),level(b)≤l},
so Aℋ(l) is sparse regarded as a standard matrix, which could be made use of as the primary data to form sparse approximate inverse. A trivial example of Aℋ(l) is that(12)AH(L(TI×I))=AH.The intuitive grasp ofAℋ(l) is shown in Figure 2. Consider an IE linear system matrix in ℋ-matrices form, the nonzero entries distribution of Aℋ(1),Aℋ(2), and Aℋ(3) are marked as shadow area corresponding to the whole matrix. Obviously high level of Aℋ(l) contains more useful data of the system.The data distribution of the finest 3 levels ofAℋ(l). This ℋ-matrices-formatted system is derived from an aircraft model which will be introduced in numerical experiments of Section 4.
(a)
Level 1(b)
Level 2(c)
Level 3Correspondingly, a set of preconditioning matrices{M(l)∣l=1,2,…L(Tℐ×ℐ)} can be generated through Aℋ(l), which satisfied the conditions below:(13)‖I-M(1)AH‖≥‖I-M(2)AH‖≥⋯≥‖I-M(L(TI×I))AH‖≈0.
In practice, M(l) can be solved by the sparse approximate inverse technique, which aims to minimize ∥I-M(1)Aℋ(l)∥F [12]. ∥·∥F is the Frobenius norm of a matrix, and M(l) is subjected to a certain sparsity pattern, expressed as(14)min⁡M(l)∈G(l)‖I-M(l)AH(l)‖F2=∑i=1nmin⁡mi(l)∈Gi(l)‖ei-mi(l)AH(l)‖22,
where ei and mi(l) are the row vectors of the matrices I and M(l), respectively. G(l) represents the sparsity pattern for different M(l). Each of the subminimization problems in the right hand of (14) can be solved independently.Becausemi(l) in (14) is constrained to a certain sparse pattern G, it has many zero entries that force the corresponding rows of Aℋ(l) to be zero in matrix-vector multiply. Let m̃i(l) denotes the subvector containing the nonzero entries of mi(l), and its corresponding rows of Aℋ(l) are denoted by Aℋ(l)i. Besides, since Aℋ(l) is sparse, the submatrix Aℋ(l)i has many columns that are identically zero. By removing the zero columns, we have a much smaller submatrix Ãℋi(l). The individual minimization problem of (14) is reduced to a least squares problem:(15)min⁡m̃i(l)‖ẽi-m̃i(l)ÃHi(l)‖2,i=1,2,…n.Aℋ(l)i is just the submatrix of Aℋ(l) that we implement the QR decomposition for solving (15).


### 3.2. Hierarchical Fast Direct Methods

If the system matrix is severely ill-conditioned and even the iterative solver with powerful preconditioner cannot obtain acceptable results, fast direct solvers are good alternative for IEℋ-matrices. According to the arithmetic of partition matrices, the inversion of a matrix can be calculated by operating its submatrices. Considering the hierarchical structure of ℋ-matrices is indeed a nested quaternary submatrices structure, and a potential direct method can be made to solve the linear system. If a system matrix A can be partitioned as(16)A=[A11A12A21A22],
then its inversion can be written as(17)A-1=[A11-1+A11-1A12C-1A21A11-1-A11-1A12C-1-C-1A21A11-1C-1],
where C=A22-A21A11-1A12. Applying this arithmetic into ℋ-matrices, we can obtain a hierarchical inversing procedure as elaborated in Algorithm 1.Algorithm 1
Recursive  procedure H_Inverse(m,x)If the matrix m is partition thenH_Inverse(m11,x11)x12=-x11⊙m12x21=-m12⊙x11m22=m22⊕m21⊙x12H_Inverse(m22,x22)m12=x12⊙m22m11=m11⊕m12⊙x21m21=m22⊙x21elseDirect_Inverse(m11)end ifThe operator⊙ and ⊕ manipulate matrices multiplication and addition which are specific for ℋ-matrices. Reference [9] gives detailed information about these operators. By implementing submatrices partition, aggregation, and truncation of singular value decomposition (SVD), the hierarchical structure can be maintained after the manipulation of these operators, and, more importantly, both CPU time and memory cost are saved compared to conventional matrices arithmetic. Consequently, the complexity of hierarchical direct solvers is O(N2) for CPU time and O(N1.5) for memory usage, in contrast with O(N3) and O(N2) of conventional direct solver. This leads to realizing kinds of fast direct solving processes including the hierarchical inversing algorithm described above.For numerical implementation, hierarchical inversing is not as fast as hierarchical LU decomposition, which is another fast direct method based on matrices decomposition [9, 13]. Considering LU decomposition for partition matrix A:(18)A=[A11A12A21A22]=[L11L21L22][U11U12U22]
matrices L and U can be figured out by the procedures:(1)
computingL11 and U11 from the LUDL11U11=A11;(2)
computingU12from L11U12=A12;(3)
computngL21from L21U11=A21;(4)
computingL22 and U22 from the LUDL22U22=A22-L21U12.Process (2) and process (3) can be refined as subprocesses of partitioned LU decomposition; therefore, by recursive applying of this procedure, the whole decomposition can be achieved. In the process of recursion, when the partitioning reaches the finest level in which submatrices are explicit expressed, conventional LU decomposition is employed. After the decomposition done, the primary L and U are still in ℋ-form and overwrite the original Aℋ. Then, we can use partitioned backward/forward processes to solve the linear system which is similar to linear equation solving by conventional LU decomposition.  For EFIE, we can utilize its symmetric property and then construct a fast LLT decomposition solver which is even much faster and memory saving. Considering LLT decomposition for partition matrix A:(19)A=[A11A21TA21A22]=[L11L21L22][L11TL21TL22T],
matrices L can be figure out by the procedures:(1)
computingL11 from the LLTL11L11T=A11;(2)
computngL21 from L21L11T=A21;(3)
computingL22from the LLTDL22L22T=A22-L21L21T.Comparing with the procedure of hierarchical LU decomposition, one step is removed because we can obtainL12T from trivial transposition. Therefore, the implementing of LLT decomposition is approximately one time faster than that of LU decomposition, in a recursive way.


## 3.1. Multilevel SAI Preconditioning

If the Krylov iterative methods are used to solve the linear system, we always expect to find a high-performance preconditioner to accelerate the convergence. Generally, instead of solving the linear system of the formAℋx=b, we solve it by the following forms:
(9)MAHx=MborAHMy=b,x=My.
In which M is a sparse matrix satisfying the condition M≈Aℋ-1 to a certain degree. Therefore, MAℋ≈I or AℋM≈I can make the iterative solver more efficient. Because of limited space, we only discuss the left preconditioning case below. The right preconditioning case can be analyzed similarly. Traditional preconditioner has a fixed form that we could only execute “single level” preconditioning. Here, we elaborate how to construct a multilevel preconditioner under ℋ-matrices format.For a block-cluster treeTℐ×ℐ in ℋ-matrices, let L(Tℐ×ℐ) denotes its depth and ℒ(Tℐ×ℐ) denotes all the leaves of the cluster tree. Any leaf block-cluster b=(t×s)∈ℒ(Tℐ×ℐ) belongs to certain level of Tℐ×ℐ; hence, 1≤level(b)≤L(Tℐ×ℐ). We define that he finest level is the 1st level, where the smallest block cluster is located. A substructure of Aℋ denoted by Aℋ(l) can be defined as follows:(10)AH(l)(i,j)={AH(i,j),(i,j)∈B(l),0,otherwise
in which B(l) is a subset of ℒ(Tℐ×ℐ), (11)B(l)={b∣b∈L(TI×I),level(b)≤l},
so Aℋ(l) is sparse regarded as a standard matrix, which could be made use of as the primary data to form sparse approximate inverse. A trivial example of Aℋ(l) is that(12)AH(L(TI×I))=AH.The intuitive grasp ofAℋ(l) is shown in Figure 2. Consider an IE linear system matrix in ℋ-matrices form, the nonzero entries distribution of Aℋ(1),Aℋ(2), and Aℋ(3) are marked as shadow area corresponding to the whole matrix. Obviously high level of Aℋ(l) contains more useful data of the system.The data distribution of the finest 3 levels ofAℋ(l). This ℋ-matrices-formatted system is derived from an aircraft model which will be introduced in numerical experiments of Section 4.
(a)
Level 1(b)
Level 2(c)
Level 3Correspondingly, a set of preconditioning matrices{M(l)∣l=1,2,…L(Tℐ×ℐ)} can be generated through Aℋ(l), which satisfied the conditions below:(13)‖I-M(1)AH‖≥‖I-M(2)AH‖≥⋯≥‖I-M(L(TI×I))AH‖≈0.
In practice, M(l) can be solved by the sparse approximate inverse technique, which aims to minimize ∥I-M(1)Aℋ(l)∥F [12]. ∥·∥F is the Frobenius norm of a matrix, and M(l) is subjected to a certain sparsity pattern, expressed as(14)min⁡M(l)∈G(l)‖I-M(l)AH(l)‖F2=∑i=1nmin⁡mi(l)∈Gi(l)‖ei-mi(l)AH(l)‖22,
where ei and mi(l) are the row vectors of the matrices I and M(l), respectively. G(l) represents the sparsity pattern for different M(l). Each of the subminimization problems in the right hand of (14) can be solved independently.Becausemi(l) in (14) is constrained to a certain sparse pattern G, it has many zero entries that force the corresponding rows of Aℋ(l) to be zero in matrix-vector multiply. Let m̃i(l) denotes the subvector containing the nonzero entries of mi(l), and its corresponding rows of Aℋ(l) are denoted by Aℋ(l)i. Besides, since Aℋ(l) is sparse, the submatrix Aℋ(l)i has many columns that are identically zero. By removing the zero columns, we have a much smaller submatrix Ãℋi(l). The individual minimization problem of (14) is reduced to a least squares problem:(15)min⁡m̃i(l)‖ẽi-m̃i(l)ÃHi(l)‖2,i=1,2,…n.Aℋ(l)i is just the submatrix of Aℋ(l) that we implement the QR decomposition for solving (15).


## 3.2. Hierarchical Fast Direct Methods

If the system matrix is severely ill-conditioned and even the iterative solver with powerful preconditioner cannot obtain acceptable results, fast direct solvers are good alternative for IEℋ-matrices. According to the arithmetic of partition matrices, the inversion of a matrix can be calculated by operating its submatrices. Considering the hierarchical structure of ℋ-matrices is indeed a nested quaternary submatrices structure, and a potential direct method can be made to solve the linear system. If a system matrix A can be partitioned as(16)A=[A11A12A21A22],
then its inversion can be written as(17)A-1=[A11-1+A11-1A12C-1A21A11-1-A11-1A12C-1-C-1A21A11-1C-1],
where C=A22-A21A11-1A12. Applying this arithmetic into ℋ-matrices, we can obtain a hierarchical inversing procedure as elaborated in Algorithm 1.Algorithm 1
Recursive  procedure H_Inverse(m,x)If the matrix m is partition thenH_Inverse(m11,x11)x12=-x11⊙m12x21=-m12⊙x11m22=m22⊕m21⊙x12H_Inverse(m22,x22)m12=x12⊙m22m11=m11⊕m12⊙x21m21=m22⊙x21elseDirect_Inverse(m11)end ifThe operator⊙ and ⊕ manipulate matrices multiplication and addition which are specific for ℋ-matrices. Reference [9] gives detailed information about these operators. By implementing submatrices partition, aggregation, and truncation of singular value decomposition (SVD), the hierarchical structure can be maintained after the manipulation of these operators, and, more importantly, both CPU time and memory cost are saved compared to conventional matrices arithmetic. Consequently, the complexity of hierarchical direct solvers is O(N2) for CPU time and O(N1.5) for memory usage, in contrast with O(N3) and O(N2) of conventional direct solver. This leads to realizing kinds of fast direct solving processes including the hierarchical inversing algorithm described above.For numerical implementation, hierarchical inversing is not as fast as hierarchical LU decomposition, which is another fast direct method based on matrices decomposition [9, 13]. Considering LU decomposition for partition matrix A:(18)A=[A11A12A21A22]=[L11L21L22][U11U12U22]
matrices L and U can be figured out by the procedures:(1)
computingL11 and U11 from the LUDL11U11=A11;(2)
computingU12from L11U12=A12;(3)
computngL21from L21U11=A21;(4)
computingL22 and U22 from the LUDL22U22=A22-L21U12.Process (2) and process (3) can be refined as subprocesses of partitioned LU decomposition; therefore, by recursive applying of this procedure, the whole decomposition can be achieved. In the process of recursion, when the partitioning reaches the finest level in which submatrices are explicit expressed, conventional LU decomposition is employed. After the decomposition done, the primary L and U are still in ℋ-form and overwrite the original Aℋ. Then, we can use partitioned backward/forward processes to solve the linear system which is similar to linear equation solving by conventional LU decomposition.  For EFIE, we can utilize its symmetric property and then construct a fast LLT decomposition solver which is even much faster and memory saving. Considering LLT decomposition for partition matrix A:(19)A=[A11A21TA21A22]=[L11L21L22][L11TL21TL22T],
matrices L can be figure out by the procedures:(1)
computingL11 from the LLTL11L11T=A11;(2)
computngL21 from L21L11T=A21;(3)
computingL22from the LLTDL22L22T=A22-L21L21T.Comparing with the procedure of hierarchical LU decomposition, one step is removed because we can obtainL12T from trivial transposition. Therefore, the implementing of LLT decomposition is approximately one time faster than that of LU decomposition, in a recursive way.


## 4. Numerical Experiments

### 4.1. Complexity of Constructingℋ-Matrices

To investigate the complexity of constructing IEℋ-matrices for electromagnetic analysis, there are two themes, namely, discretization density varying and wave frequencies varying. The former one means we change the number of unknowns along with varying the mesh density regarding electromagnetic wavelength. And the latter one means we change the number of unknowns along with varying wave frequencies, but the mesh density regarding wavelength is fixed. Here, we use a PEC sphere model of 1.0λ (wavelength) radius to test both these two themes, as shown in Figures 3 and 4.Figure 3
Computational complexity of constructingℋ-matrices by fixed incident wave at 300.0 MHz.Figure 4
Computational complexity of constructingℋ-matrices by fixed discretization mesh density at 8 points per wavelength.From the investigation, we can easily see that if incident wave frequency is fixed, the complexity of both CPU time and memory usage approach toO(NlogN), where N regards to the number of unknowns. If we fix the discretization density and change the frequency of incident wave, the complexity curve is close to O(N4/3logN). This is because the object domain contains more phase information when incident by higher frequency wave; therefore, the compression ratio of low rank submatrices in ℋ-matrices is lower in higher frequency cases.


### 4.2. Iterative Solving with Multilevel SAI Preconditioner

Firstly, a conducting sphere is used to demonstrate the improvement of the spectrum characteristics of the linear coefficient matrix by employing multilevel-SAI preconditioner (ML-SAI). SupposingM(l) is the lth level preconditioner and xrnd(i) is the ith randomized vector, we define the regression index cr(p)of lth preconditioning level as(20)cr(p)(l)=1p∑i=1p‖xrnd(i)-M(l)AHxrnd(i)‖2‖xrnd(i)‖2,
if M(l)Aℋ is close to identity matrix I, then cr(p)(l)→0 is expected. Particularly, M(0)=I, which means no preconditioner is imposed. A sphere of 2λ diameter, discretized by Rao-Wilton-Gliso (RWG) basis with different densities, is used to obtain the linear system matrix Aℋ. The different regression indices of different preconditioning levels are presented in Table 1 and compared to analogous indices of conventional SAI preconditioning in MLFMA.Table 1
The regression index of different levels of ML-SAI and conventional SAI preconditioning.
Regression indexNumber of unknowns9573,97216,47365,892ℋ-Mat.cr(100)(0)5.5682.2971.3711.042cr(100)(1)0.13110.24400.31550.4237cr(100)(2)0.089250.20340.25230.2816cr(100)(3)1.116E-60.11120.18280.2165ML-FMANone5.7232.4151.3550.989SAI0.10520.22370.28440.3125From Table1, we can see that with the preconditioning level of ML-SAI increased, the regression index is decreased, which means the effect of preconditioning is getting better. And the preconditioning effect of conventional SAI in MLFMA is better than the low level ML-SAI but not as good as high level ML-SAI. It should be noted that the cr(100)(3) in case of 957 unknowns is close to zero because the ℋ-matrices structure of this model is only 3 levels, which means M(3) is very close to the exact inverse of the linear system.Next, a40λ diameter PEC sphere is solved by ℋ-matrices method. Combine field integral equation (CFIE) is used, and the number of unknowns is 426,827. To set up the hierarchical system matrix, the memory cost is 42.6 GBytes and 26,387 seconds are used by 8 core parallel computing, on a workstation of Intel Xeon X5460 processor with 64 GBytes RAM. Generalized minimal residual (GMRES) [14] algorithm is employed as the Krylov iterative solver, and we use 1~3 level SAI preconditioners and diagonal block preconditioners (DBPs) as reference to accelerate the convergence of iteration. From Figure 5, we can see that the RCS curve conform to the Mie series identically. And the iterative history in Figure 6 shows that ML-SAIs have a distinct impact on accelerating the convergence of iteration.Figure 5
The bistatic RCS of a conducting sphere with 40λ.Figure 6
The iterative history of solving the PEC sphere cases. GMRES(30) is used and accelerated by different preconditioners.Another example is an aircraft model which is50.75λ long, 29.20λ wide, and 13.57λ high, excited by a plane wave incoming from its nose with an upper 45° angle to the center axis of its main body. The CFIE is used, and there are 412,737 unknowns to simulate the exciting surface currents as shown in Figure 7. By solving the linear system by ℋ-matrices method, there are 14 levels in the block-cluster tree. To set up the hierarchical system matrix, the memory cost is 38.4 GBytes and 25,634 seconds are used by 8 core parallel computing, on a workstation of Intel Xeon X5460 processor with 64 GBytes RAM. The data of Aℋ(l) in the finest 3 levels are to construct preconditioners M(1),M(2), and M(3), respectively.Figure 7
The surface current distribution of the aircraft model obtained by solving the ML-SAI preconditionedℋ-matrices method. The unit of surface current is A/m2.GMRES(90) is used to solve the linear system. The convergence histories of iteration with different preconditioners are presented in Figure8. This is a little tough case because of complicated structure and relatively high frequency.  If no preconditioner or just DBP is employed, the convergence is very poor. By employing ML-SAI, the iteration converges quickly under satisfied residual error, and the convergence is much better with the higher level preconditioning processed. The conventional SAI with MLFMA is also able to achieve the request, but its convergence is not as good as higher level ML-SAI. Table 2 shows the solving time and memory cost of ML-SAI, conventional SAI with MLFMA, and diagonal blocks preconditioner. The memory cost is mainly constituted by the storage of preconditioner M and Krylov subspace vectors of GMRES. In the fact that the restart number of GMRES, namely, the number of subspace vectors, is chosen the same for all cases, so the memory cost can be regarded as the scale of preconditioner M.Table 2
The time and memory cost of iterative solving aircraft model with different preconditioners.
Time (seconds)Memory (MBytes)ML-SAI level 1582.8325.4ML-SAI level 2234.5464.5ML-SAI level 3184.3786.5SAI with MLFMA643.7312.3DBP1134.2213.9Figure 8
The convergence histories of GMRES(90) with no preconditioner, DB preconditioner, and ML-SAI preconditioners of different levels.


### 4.3. Hierarchical Fast Direct Solvers

In this part, we give some numerical examples solved by fastℋ-matrices direct solvers. We use hierarchical LU decomposition to solve CFIE and use hierarchical LLT decomposition to solve EFIE, since utilizing its symmetric property. First, a PEC sphere model is used to investigate the accuracy of fast direct solvers. From Figure 9, we can see that the bistatic RCSs obtained by both fast LU and LLT decomposition solvers have very high precision regarding to analytic Mie series. The EFIE result has a slight bigger relative error than CFIE because the system matrix of EFIE is not diagonal dominated, consequently, the numerical error accumulated through LLT decomposition is more than that in CFIE.Figure 9
The bistatic RCS solved by hierarchical fast direct solvers. CFIE is solved by LU decomposition solver, and EFIE is solved byLLT decomposition.Next, the computational complexity of fastLLT decomposition solver is tested via the PEC sphere model. The number of unknowns is varying along with the change of frequency of the incident wave, and the mesh discretization density is fixed at 10.0 points per wavelength. The test results are shown in Figures 10 and 11, respectively, in which we can observe that the complexity of CPU time is close to O(N2) and memory usage is close to O(N1.5). This is much more efficient than conventional direct solvers.Figure 10
CPU time complexity of hierarchicalLLT decomposition.Figure 11
Memory usage complexity of hierarchicalLLT decomposition.The aircraft model shown in the previous part is also used here to test fastLLT decomposition solver by employing EFIE. The incident wave is 10.0 GHz, and mesh discretization density is 10.0 points per wave number. The number of unknowns is 119,202, and 8 levels are set for ℋ-matrices. This example is parallelly solved by a workstation equipped with dual Intel Xeon X5460 processor and 64 GBytes of RAM. 2,256 seconds is used for building up the hierarchical system matrix, and the memory cost is 5,128.4 MBytes. LLT decomposition time is 8,415 seconds, and 6,437.2 MBytes are used after decomposition. Figure 12 presents the surface current solved by LLT decomposition.Figure 12
The surface current distribution of the aircraft model obtained by solving EFIE with hierarchicalLLT decomposition. The frequency of incident wave is at 10.0 GHz. The unit of surface current is A/m2.The most notable advantage of fast direct solver is the high efficiency of handling multiple right-hand-side cases. After LU orLLT decomposition, solving multiple right-hand vectors is only one-round process of forward/backward substitution, so that it is much faster than iterative solvers, which gives a full-round iteration for each vector step by step. In IE electromagnetic analysis, monostatic RCS is the typical multiple right-hand-side problem. We here use LLT decomposition to solve the monostatic RCS problem of the aircraft model incident by 3.0 GHz electromagnetic wave. The range of sweep angles is 180°, and 181 samples are calculated. The total solving time is 468.3 seconds, contrast to MLFMA combined with GMRES iterative solver, which costs 4215.7 seconds, almost ten times of the former one. In order to compare to LLT decomposition, here, we use EFIE for iterative MLFMA. Figure 13 shows the comparing results between LLT decomposition and iteration of MLFMA, we can observe that direct solved curve is smoother than that of iterative solved and there is a big difference between them. We believe that the result of LLT decomposition is more accurate and the reason is elaborated as follow. Theoretically, the monostatic RCS curve of this intermediate-frequency example should be smooth, because the sampling is sufficiently continuous comparing to the wavelength. However, for EFIE iteration, the convergence is poor due to the large condition number of its system matrix. Therefore, a loose relative residual error threshold 0.005 is set here, for the purpose of getting iterative result with a reasonable CPU time cost, and this causes the iterative result not very accurate.Figure 13
The monostatic RCS of the aircraft model. The incident wave is 3.0 GHz, HH-polarization.


## 4.1. Complexity of Constructingℋ-Matrices

To investigate the complexity of constructing IEℋ-matrices for electromagnetic analysis, there are two themes, namely, discretization density varying and wave frequencies varying. The former one means we change the number of unknowns along with varying the mesh density regarding electromagnetic wavelength. And the latter one means we change the number of unknowns along with varying wave frequencies, but the mesh density regarding wavelength is fixed. Here, we use a PEC sphere model of 1.0λ (wavelength) radius to test both these two themes, as shown in Figures 3 and 4.Figure 3
Computational complexity of constructingℋ-matrices by fixed incident wave at 300.0 MHz.Figure 4
Computational complexity of constructingℋ-matrices by fixed discretization mesh density at 8 points per wavelength.From the investigation, we can easily see that if incident wave frequency is fixed, the complexity of both CPU time and memory usage approach toO(NlogN), where N regards to the number of unknowns. If we fix the discretization density and change the frequency of incident wave, the complexity curve is close to O(N4/3logN). This is because the object domain contains more phase information when incident by higher frequency wave; therefore, the compression ratio of low rank submatrices in ℋ-matrices is lower in higher frequency cases.


## 4.2. Iterative Solving with Multilevel SAI Preconditioner

Firstly, a conducting sphere is used to demonstrate the improvement of the spectrum characteristics of the linear coefficient matrix by employing multilevel-SAI preconditioner (ML-SAI). SupposingM(l) is the lth level preconditioner and xrnd(i) is the ith randomized vector, we define the regression index cr(p)of lth preconditioning level as(20)cr(p)(l)=1p∑i=1p‖xrnd(i)-M(l)AHxrnd(i)‖2‖xrnd(i)‖2,
if M(l)Aℋ is close to identity matrix I, then cr(p)(l)→0 is expected. Particularly, M(0)=I, which means no preconditioner is imposed. A sphere of 2λ diameter, discretized by Rao-Wilton-Gliso (RWG) basis with different densities, is used to obtain the linear system matrix Aℋ. The different regression indices of different preconditioning levels are presented in Table 1 and compared to analogous indices of conventional SAI preconditioning in MLFMA.Table 1
The regression index of different levels of ML-SAI and conventional SAI preconditioning.
Regression indexNumber of unknowns9573,97216,47365,892ℋ-Mat.cr(100)(0)5.5682.2971.3711.042cr(100)(1)0.13110.24400.31550.4237cr(100)(2)0.089250.20340.25230.2816cr(100)(3)1.116E-60.11120.18280.2165ML-FMANone5.7232.4151.3550.989SAI0.10520.22370.28440.3125From Table1, we can see that with the preconditioning level of ML-SAI increased, the regression index is decreased, which means the effect of preconditioning is getting better. And the preconditioning effect of conventional SAI in MLFMA is better than the low level ML-SAI but not as good as high level ML-SAI. It should be noted that the cr(100)(3) in case of 957 unknowns is close to zero because the ℋ-matrices structure of this model is only 3 levels, which means M(3) is very close to the exact inverse of the linear system.Next, a40λ diameter PEC sphere is solved by ℋ-matrices method. Combine field integral equation (CFIE) is used, and the number of unknowns is 426,827. To set up the hierarchical system matrix, the memory cost is 42.6 GBytes and 26,387 seconds are used by 8 core parallel computing, on a workstation of Intel Xeon X5460 processor with 64 GBytes RAM. Generalized minimal residual (GMRES) [14] algorithm is employed as the Krylov iterative solver, and we use 1~3 level SAI preconditioners and diagonal block preconditioners (DBPs) as reference to accelerate the convergence of iteration. From Figure 5, we can see that the RCS curve conform to the Mie series identically. And the iterative history in Figure 6 shows that ML-SAIs have a distinct impact on accelerating the convergence of iteration.Figure 5
The bistatic RCS of a conducting sphere with 40λ.Figure 6
The iterative history of solving the PEC sphere cases. GMRES(30) is used and accelerated by different preconditioners.Another example is an aircraft model which is50.75λ long, 29.20λ wide, and 13.57λ high, excited by a plane wave incoming from its nose with an upper 45° angle to the center axis of its main body. The CFIE is used, and there are 412,737 unknowns to simulate the exciting surface currents as shown in Figure 7. By solving the linear system by ℋ-matrices method, there are 14 levels in the block-cluster tree. To set up the hierarchical system matrix, the memory cost is 38.4 GBytes and 25,634 seconds are used by 8 core parallel computing, on a workstation of Intel Xeon X5460 processor with 64 GBytes RAM. The data of Aℋ(l) in the finest 3 levels are to construct preconditioners M(1),M(2), and M(3), respectively.Figure 7
The surface current distribution of the aircraft model obtained by solving the ML-SAI preconditionedℋ-matrices method. The unit of surface current is A/m2.GMRES(90) is used to solve the linear system. The convergence histories of iteration with different preconditioners are presented in Figure8. This is a little tough case because of complicated structure and relatively high frequency.  If no preconditioner or just DBP is employed, the convergence is very poor. By employing ML-SAI, the iteration converges quickly under satisfied residual error, and the convergence is much better with the higher level preconditioning processed. The conventional SAI with MLFMA is also able to achieve the request, but its convergence is not as good as higher level ML-SAI. Table 2 shows the solving time and memory cost of ML-SAI, conventional SAI with MLFMA, and diagonal blocks preconditioner. The memory cost is mainly constituted by the storage of preconditioner M and Krylov subspace vectors of GMRES. In the fact that the restart number of GMRES, namely, the number of subspace vectors, is chosen the same for all cases, so the memory cost can be regarded as the scale of preconditioner M.Table 2
The time and memory cost of iterative solving aircraft model with different preconditioners.
Time (seconds)Memory (MBytes)ML-SAI level 1582.8325.4ML-SAI level 2234.5464.5ML-SAI level 3184.3786.5SAI with MLFMA643.7312.3DBP1134.2213.9Figure 8
The convergence histories of GMRES(90) with no preconditioner, DB preconditioner, and ML-SAI preconditioners of different levels.


## 4.3. Hierarchical Fast Direct Solvers

In this part, we give some numerical examples solved by fastℋ-matrices direct solvers. We use hierarchical LU decomposition to solve CFIE and use hierarchical LLT decomposition to solve EFIE, since utilizing its symmetric property. First, a PEC sphere model is used to investigate the accuracy of fast direct solvers. From Figure 9, we can see that the bistatic RCSs obtained by both fast LU and LLT decomposition solvers have very high precision regarding to analytic Mie series. The EFIE result has a slight bigger relative error than CFIE because the system matrix of EFIE is not diagonal dominated, consequently, the numerical error accumulated through LLT decomposition is more than that in CFIE.Figure 9
The bistatic RCS solved by hierarchical fast direct solvers. CFIE is solved by LU decomposition solver, and EFIE is solved byLLT decomposition.Next, the computational complexity of fastLLT decomposition solver is tested via the PEC sphere model. The number of unknowns is varying along with the change of frequency of the incident wave, and the mesh discretization density is fixed at 10.0 points per wavelength. The test results are shown in Figures 10 and 11, respectively, in which we can observe that the complexity of CPU time is close to O(N2) and memory usage is close to O(N1.5). This is much more efficient than conventional direct solvers.Figure 10
CPU time complexity of hierarchicalLLT decomposition.Figure 11
Memory usage complexity of hierarchicalLLT decomposition.The aircraft model shown in the previous part is also used here to test fastLLT decomposition solver by employing EFIE. The incident wave is 10.0 GHz, and mesh discretization density is 10.0 points per wave number. The number of unknowns is 119,202, and 8 levels are set for ℋ-matrices. This example is parallelly solved by a workstation equipped with dual Intel Xeon X5460 processor and 64 GBytes of RAM. 2,256 seconds is used for building up the hierarchical system matrix, and the memory cost is 5,128.4 MBytes. LLT decomposition time is 8,415 seconds, and 6,437.2 MBytes are used after decomposition. Figure 12 presents the surface current solved by LLT decomposition.Figure 12
The surface current distribution of the aircraft model obtained by solving EFIE with hierarchicalLLT decomposition. The frequency of incident wave is at 10.0 GHz. The unit of surface current is A/m2.The most notable advantage of fast direct solver is the high efficiency of handling multiple right-hand-side cases. After LU orLLT decomposition, solving multiple right-hand vectors is only one-round process of forward/backward substitution, so that it is much faster than iterative solvers, which gives a full-round iteration for each vector step by step. In IE electromagnetic analysis, monostatic RCS is the typical multiple right-hand-side problem. We here use LLT decomposition to solve the monostatic RCS problem of the aircraft model incident by 3.0 GHz electromagnetic wave. The range of sweep angles is 180°, and 181 samples are calculated. The total solving time is 468.3 seconds, contrast to MLFMA combined with GMRES iterative solver, which costs 4215.7 seconds, almost ten times of the former one. In order to compare to LLT decomposition, here, we use EFIE for iterative MLFMA. Figure 13 shows the comparing results between LLT decomposition and iteration of MLFMA, we can observe that direct solved curve is smoother than that of iterative solved and there is a big difference between them. We believe that the result of LLT decomposition is more accurate and the reason is elaborated as follow. Theoretically, the monostatic RCS curve of this intermediate-frequency example should be smooth, because the sampling is sufficiently continuous comparing to the wavelength. However, for EFIE iteration, the convergence is poor due to the large condition number of its system matrix. Therefore, a loose relative residual error threshold 0.005 is set here, for the purpose of getting iterative result with a reasonable CPU time cost, and this causes the iterative result not very accurate.Figure 13
The monostatic RCS of the aircraft model. The incident wave is 3.0 GHz, HH-polarization.


## 5. Conclusion

The hierarchical matrice methods presented in this paper is embedded in electromagnetic IE method. Due to its special structure, we can construct a multilevel SAI preconditioner to accelerate the convergence of iterative solving, and even kinds of fast direct solvers can be made, which is not viable for the traditional IE fast algorithm. The multilevel SAI preconditioner proposed here is more efficient than conventional “single level” preconditioners, and hierarchical fast direct solvers are good alternatives to iterative solvers, very suitable for ill-conditioned system and multiple right-hand-side problems. Furthermore, the kernel independence feature of hierarchical matrices method is adapted to varied electromagnetic problems without being limited to integral equation with free-space Green function.



---

*Source: 756259-2012-02-01.xml*