# Efficient and Memory Saving Method Based on Pseudoskeleton Approximation for Analysis of Finite Periodic Structures.

1. Introduction

Periodic structures have recently found wide applications in the electromagnetic engineering such as antenna arrays and metamaterials with negative permittivity and negative permeability. Hence, the accurate and efficient analysis of periodic structures becomes quite essential. If the periodic structure is an infinite array, simple methods can be applied based on Floquet's theorem  or periodic Green's function , where only a single cell of the periodic structure is the domain of interest.

However, all periodic structures have finite size in real-life problems, although the size may be very large. Therefore, the numerical algorithms which accurately consider the mutual coupling between all cells should be used if the accurate results are required. The method of moments (MoM)  and its fast algorithms such as fast multipole method (FMM) , adaptive cross approximation (ACA) , and FFT-based methods  are flexible approaches to study the surface problems. However, the efficiency of numerical methods is still limited since the periodic property is not used in the algorithm framework. Recently, a hybrid method combining the accurate MoM and periodic method of moment (PMM)  has been proposed which can gain the balance between the two methods. Moreover, some physically based entire-domain basis functions  have been developed to reduce the number of unknowns. Further, the FMM and FFT techniques are integrated to accelerate the calculation .

Compared to the ACA method, pseudoskeleton approximation (PSA)  is also an efficient low-rank-based algebraic fast algorithm which makes it a really competitive alternative. In this paper, we propose an efficient method with low-memory requirement based on PSA to perform the analysis for finite periodic structures effectively and accurately. In consideration of the accuracy of the mutual interactions  and the simplicity of implementation, our proposed method uses the formulations derived from the local basis functions instead of macro basis function (MBF) [11, 12]. In this paper, PSA is not only used to accelerate the matrix-vector product (MVP) inside the single unit but also adopted to decrease the calculation burden of the coupling between the different cells. Moreover, the number of decomposed coupling matrices is minimized due to the displacement invariance of the periodic property. With these improvements, an efficient method with low-memory usage of finite periodic objects can be achieved. Several numerical examples are given to show the priority of the proposed method compared to the conventional multilevel fast multipole algorithm (MLFMA)  for periodic structures.

2. MoM and PSA Formulation

In this section, the basis principle of MoM and PSA is briefly introduced at first. Then, the choice of arguments in PSA is discussed.

2.1. MoM Equations and Its Fast Algorithms. Consider a time-harmonic [e.sup.-j[omega]t] electromagnetic wave scattering or radiation problem of an arbitrary perfect electrically conducting (PEC) object. The object is excited by an incident electric field [E.sup.inc](r), then the electric field integral equation (EFIE) associated with the surface equivalent currents J(r) can be expressed by

- j/[omega][mu] [[??] * [E.sup.inc](r)] = [??] * [[integral].sub.S'] [1 + 1/[k.sup.2] [nabla][nabla]] J(r')G(r,r')dr', (1)

where [??] is the tangential unit vector on the surface S' of the object and [omega] and [mu] stand for the angular frequency and permeability, respectively. G(r,r') = [e.sup.-jk[absolute value of r-r']]/(4[pi] [absolute value of r-r]) is 3D scalar Green's function. The linear system of MoM is obtained by discretizing the unknown vector J(r) with Rao-Wilton-Glisson (RWG)  basis functions and applying Galerkin's testing method. Let ZI = V represent the above EFIE matrix system. The MVP process can be accelerated by fast algorithms which can be written as follows:

ZI = [Z.sub.near]I + [Z.sub.far]I, (2)

where [Z.sub.near] is the matrix of near field interactions which are directly computed and stored and [Z.sub.far] stands for couplings between the far-field interactions which will be accelerated together with [Z.sub.far]I in the iterative solving process.

2.2. Basic PSA Frameworks. According to the low-rank decomposition, the far-field interaction matrix [Z.sub.far] with rank-deficient property can be approximated by a product of two much smaller submatrices U and V:

[Z.sup.mxn.sub.far] [approximately equal to] [U.sup.mxn][V.sup.rxn] (3)

where r is the effective rank and satisfies r [much less than] m, n. While in the skeleton approximation (SA) theory , there is a nonsingular r x r submatrix [??] in [Z.sub.far]. Denote the rectangular matrices as C and R which contain the columns and rows of [??], respectively, then [Z.sub.far] is expressed as

[Z.sub.far] = C[[??].sup.-1]R, (4)

where C and R have dimensions with m x r and r x n, respectively. In the PSA method, [Z.sub.far] is reevaluated as

[Z.sub.far] [approximately equal to] CGR, (5)

where G is the pseudoinverse of [??] with dimensions of p x p(p > r), then p columns and p rows are chosen from [Z.sub.far] to get the C and R. Three aspects need to be noted here: (i) [[??].sup.-1] is the inverse of [??] and the computation of inverse is very expensive; (ii) the determination of the value of p is a balance between accuracy and efficiency, where p is a number large enough so that r most important bases will be embedded; (iii) how to choose those working columns and rows of C and R.

For the (ii) and (iii) problems, we will discuss them in the next subsection. For the first problem, assuming that [??] can be decomposed via singular value decomposition (SVD) as

[??] = P[summation][Q.sup.*], (6)

where P and Q are p x p unitary matrices, [SIGMA] is a diagonal p x p matrix with nonnegative real numbers, and [Q.sup.*] is the complex conjugate transpose of Q. In the actual implementation, the dimension of P, Q, and [SIGMA] can be further decreased by a preset threshold . Let [??], [??], and [??] represent the reduced submatrix of P, [SIGMA], and Q, respectively. Then, calculation of the pseudoinverse of [??] (i.e., G) is straightforward:

[mathematical expression not reproducible]. (7)

By combining (5) and (7), the original far-field interaction matrix [Z.sub.far] can be decomposed as

[mathematical expression not reproducible]. (8)

2.3. Choice of Arguments in PSA. As mentioned in the previous subsection, the choice of the value of p and the specific sampling rows and columns determine the performance of PSA. In the randomized PSA (RPSA) , p is equal to 2r. Then, the problem of (ii) transforms into how to estimate the rank r of the original matrix. In , Chai and Jiao give the approximation of rank of the 3D EM problem by [Rank.sub.3D]~O([k.sub.0]), where [k.sub.0] is the studying wave number. In this paper, the rank is approximated by

[Rank.sub.3D] [approximately equal to] [k.sub.0]d + [alpha] ln ([pi] + [k.sub.0]d), (9)

where d is the diameter of the bounding box corresponding to the octree structure and [alpha] is a preset positive parameter. The larger the [alpha] is, the more accurate the matrix decomposition is. In this paper, when [alpha] is set as 3, the satisfactory accuracy can be guaranteed.

For the problem of (iii), instead of using random numbers in RPSA, we use a strategy analogous to ACA in this paper. Firstly, initialize from the 0th row as the first row index. Then, find the largest entry in this row, and the corresponding column value in which this entry is located is chosen as the next column index. Similarly, find the largest entry in the current column and get the next row index which should be different from all previous row indexes. This process is carried out iteratively until p rows and columns are found and stored. Please refer to  for more information.

3. Proposed Method for Periodic Structures

We consider a case of arbitrarily shaped PEC patch, for example, refer in Figure 1. The total impedance matrix contains two parts: self-coupling blocks and mutual coupling blocks. Therefore, both the blocks are analyzed and decomposed by PSA to gain better efficiency and low-memory usage.

3.1. PSA for Decomposition of Self-Coupling Matrix. In the previous mentioned periodic algorithms, the self-coupling block elements are calculated by a direct MoM method which is not efficient for a large array unit. Therefore, PSA is used to decompose the self-coupling matrix. As the same to all fast algorithms of MoM, the low-rank decompositions based on PSA are performed on the far-field groups while the near-field elements are directly computed and saved. Moreover, some operations are taken out to further improve the efficiency of PSA. Firstly, the pseudoinverse of the Z (i.e., G) will not be directly calculated since the complexity of SVD is very high. Instead, the LU decomposition will be performed when the dimension grows up. Different from the previously proposed PSA, C, R, and the LU decomposition of [??] will be stored. Moreover, C and R in (8) are further decomposed by ACA technique. Finally, the original far-field coupling matrix can be decomposed into six submatrices which is showed by the following formula:

[Z.sub.far] = [U.sub.C][V.sub.C] [([L.sub.[??]][U.sub.[??]]).sup.-1][U.sub.R][V.sub.R], (10)

where [U.sub.C][V.sub.C] (or [U.sub.R][V.sub.R]) is the ACA decomposition of C (or R) and [L.sub.[??]][U.sub.[??]] is the LU decomposition of [??]. Note that the inverse of LU matrix of [??] is not directly computed. Actually, the LU back substitute is performed after the MVP based on the matrix [U.sub.R][V.sub.R] in each MVP process.

3.2. PSA for Decomposition of Mutual Coupling Matrices. In the traditional fast algorithms, the matrix decompositions are only performed on the far-field groups especially for the physically related methods such as MLFMA. However, the mutual interactions between two different cells (off-diagonal blocks) are all computed by PSA in our implementation. This treatment has been also used in  and demonstrated to be more efficient but losing little accuracy. Moreover, the displacement invariance property is explored in our implementation. Since the mutual coupling remains the same when the distance between two groups' centers is not changed, the corresponding matrix will be calculated and stored only once under the index of relative distance. In the MVP process, the stored coupling matrices may be used several times when the relative distances are the same. By using this scheme, the number of mutual coupling matrices can be decreased dramatically. For example, the original number of off-diagonal matrix blocks of Figure 1 is 20 x 20 - 20 = 380 while the reduced number will be 62 if the displacement invariance is used.

3.3. Preconditioner Considerations. The preconditioning stage should be also considered carefully while dealing with complex structures. In this paper, the preconditioning matrix M for impedance matrix Z is built based on the self-coupling matrix [M.sub.0].

M = diag [[M.sub.0], [M.sub.0], ..., [M.sub.0]]. (11)

Moreover, the inner-outer iteration scheme can be also applied when facing ill-conditioned matrix systems such as EFIE for radiation analysis. In our implementation, the inner iteration is the solution of the self-coupling matrices. Certainly, the inner solution area can be also extended to improve the convergence rate of outer iteration.

4. Numerical Results

In this section, several numerical experiments are presented to demonstrate the efficiency and low-memory requirement of the proposed method. For all the simulations, the mesh sizes are no less than 0.1 [lambda] ([lambda] represents the wave length). The biconjugate gradient stabilized method (BiCGSTAB) is adopted as the iterative solver for the matrix equations, and the threshold of the iterative stopping criterion for residuum is set to 0.001. The sparse approximate inverse preconditioner is used in all simulations. All the computations were carried out on a workstation with four Xeon E5-4620 CPU and 256 GB of RAM with OpenMP technique, and the digits were stored in double precision.

4.1. Scattering Simulation. In order to verify the accuracy and efficiency of the proposed method, we first consider a periodic structure consisting of [N.sub.0] = 4x5 = 20 element cells (also shown in Figure 1), where the unit cell is a 0.5 m PEC sphere. The working frequency is set to be 1 GHz which leads to 20286 unknowns (degrees of freedom) in each patch. Hence, the total number of unknowns is N = 20 x 20286 = 416520. The gap between two unit cells is set to be 0.5 m. The radar cross sections (RCS) computed by MLFMA, the proposed method (periodic PSA), and the commercial software (FEKO) are illustrated in Figure 2(a). It could be clearly seen that the numerical results from the proposed method have excellent agreement with both the conventional MLFMA and FEKO. While as in Table 1, both the computation time and peak memory usage by periodic PSA are less than the MLFMA method. What is more, another case is considered when the gap between two unit cells is decreased to be 0.1 m. The RCS results and the use of computer resources are shown in Figure 2(b) and Table 1, respectively. Although the rank of mutual coupling matrix gets bigger due to the tight coupling interactions, the proposed method still shows the good efficiency and low-memory usage with losing little accuracy.

Furthermore, the scattering of a large-scale periodic structure consisting of [N.sub.0] = 32 x 32 = 1024 element cells is considered. The working frequency is set to be 1 GHz. Different radii of the sphere are considered which lead to about 0.81, 3.38, and 5.46 million numbers of unknowns corresponding to 0.1 m, 0.2 m, and 0.25 m models, respectively. The distance between two centers of neighboring sphere is set to be a fixed value of 1.0 m. Table 2 shows the information of models and the numerical performance between MLFMA and the periodic PSA. It can be seen that for a large-scale problem, the proposed PSA method still needs less computation time and much less memory requirement. It should be noted that in the 3rd example in Table 2 when the sphere radius is set to 0.25 m, the finest box in MLFMA has to be set to 0.20A to avoid low-frequency breakdown. Figure 3 shows that the RCS results of 0.25 m sphere models computed from the periodic PSA still agree well with FEKO and MLFMA.

4.2. Radiation Simulation. Lastly, we consider a radiation problem, which contains [N.sub.0] = 10 x 10 = 100 antenna cells with delta-gap excitations. Each cell has 1145 unknowns, and the total number is N = 100 x 1145 = 114500. The gap between two cells in this model is 0.005 m which is smaller than 0.1 [lambda]. Since EFIE is the only choice for radiation problems, the inner-outer iteration scheme is adopted to accelerate the calculation process. Figure 4 shows the current density distribution of the antenna array. The far-field gain patterns computed by the proposed method and MLFMA are shown in Figure 5.

Obviously, the two methods give nearly the same results. However, as in Table 3, the computation time has been reduced from 419 seconds to 172 seconds in solving matrix equations for the radiation problem. Moreover, the memory usage for MLFMA is 8841 MB while the value is only 358 MB for the proposed method. Again, the proposed method performs much better than MLFMA. The first reason is that only 360 off-diagonal blocks are need to be calculated and stored in our method while the original value is 100 x 100 - 100 = 9900. Secondly, the mesh size is about 0.01 [lambda] for the antenna cell (dense mesh). As it is known, the lowest-level box size is limited to be no less than 0.20 [lambda] since MLFMA suffers from the low-frequency breakdown, which results in a heavy near-field matrix consumption.

5. Conclusion

In this paper, an efficient and memory saving method based on pseudoskeleton approximation (PSA) for the effective and accurate analysis of finite periodic structures is presented. The appropriate choice of sampling rows and columns as well as sampling number in PSA is discussed and confirmed. Based on the algebraic fast algorithm PSA, we have carefully handled the self-coupling and mutual-coupling blocks of the original impedance matrix. Moreover, the number of decomposed coupling matrices is minimized by the employment of the displacement invariance of the periodic property. The simulation results of large-scale scattering and radiation examples show that the proposed periodic PSA method needs less computation time and much less memory compared to conventional MLFMA. Hence, numerical results demonstrate the efficiency and superiority of the proposed method.

Data Availability

No additional data are available.

https://doi.org/10.1155/2018/1612498

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China under Grant 61171035.

References

 F. Xu, Y. Zhang, W. Hong, K. Wu, and T. J. Cui, "Finite-difference frequency-domain algorithm for modeling guided-wave properties of substrate integrated waveguide," IEEE Transactions on Microwave Theory and Techniques, vol. 51, no. 11, pp. 2221-2227, 2003.

 K. Yasumoto and K. Yoshitomi, "Efficient calculation of lattice sums for free-space periodic Green's function," IEEE Transactions on Antennas and Propagation, vol. 47, no. 6, pp. 1050-1055, 1999.

 R. F. Harrington, Field Computation by Moment Methods, Wiley-IEEE Press, 1993.

 C.-C. Lu and W. C. Chew, "A multilevel algorithm for solving a boundary integral equation of wave scattering," Microwave and Optical Technology Letters, vol. 7, no. 10, pp. 466-470, 1994.

 M. Bebendorf, "Approximation of boundary element matrices," Numerische Mathematik, vol. 86, no. 4, pp. 565-589, 2000.

 E. Bleszynski, M. Bleszynski, and T. Jaroszewicz, "Aim: adaptive integral method for solving large-scale electromagnetic scattering and radiation problems," Radio Science, vol. 31, no. 5, pp. 1225-1251, 1996.

 J. Su, X. Xu, and B. Hu, "Hybrid PMM-MoM method for the analysis of finite periodic structures," Journal of Electromagnetic Waves and Applications, vol. 25, no. 2-3, pp. 267-282, 2011.

 W. B. Lu, T. J. Cui, Z. G. Qian, X. X. Yin, and W. Hong, "Accurate analysis of large-scale periodic structures using an efficient sub-entire-domain basis function method," IEEE Transactions on Antennas and Propagation, vol. 52, no. 11, pp. 3078-3085, 2004.

 W. B. Lu, T. J. Cui, and H. Zhao, "Acceleration of fast multipole method for large-scale periodic structures with finite sizes using sub-entire-domain basis functions," IEEE Transactions on Antennas and Propagation, vol. 55, no. 2, pp. 414-421, 2007.

 S. A. Goreinov, N. L. Zamarashkin, and E. E. Tyrtyshnikov, "Pseudo-skeleton approximations by matrices of maximal volume," Mathematical Notes, vol. 62, no. 4, pp. 515-519, 1997.

 E. Suter and J. R. Mosig, "A subdomain multilevel approach for the efficient MoM analysis of large planar antennas," Microwave and Optical Technology Letters, vol. 26, no. 4, pp. 270-277, 2000.

 D. Gonzalez-Ovejero, F. Mesa, and C. Craeye, "Accelerated macro basis functions analysis of finite printed antenna arrays through 2D and 3D multipole expansions," IEEE Transactions on Antennas and Propagation, vol. 61, no. 2, pp. 707-717, 2013.

 J. Song, C.-C. Lu, and W. C. Chew, "Multilevel fast multipole algorithm for electromagnetic scattering by large complex objects," IEEE Transactions on Antennas and Propagation, vol. 45, no. 10, pp. 1488-1493, 1997.

 S. Rao, D. Wilton, and A. Glisson, "Electromagnetic scattering by surfaces of arbitrary shape," IEEE Transactions on Antennas and Propagation, vol. 30, no. 3, pp. 409-418, 1982.

 X. Zhu and W. Lin, "Randomised pseudo-skeleton approximation and its application in electromagnetics," Electronics Letters, vol. 47, no. 10, pp. 590-592, 2011.

 W. Chai and D. Jiao, "Theoretical study on the rank of integral operators for broadband electromagnetic modeling from static to electrodynamic frequencies," IEEE Transactions on Components, Packaging and Manufacturing Technology, vol. 3, no. 12, pp. 2113-2126, 2013.

 Y. Zhang and H. Lin, "Localized pseudo-skeleton approximation method for electromagnetic analysis on electrically large objects," Progress In Electromagnetics Research Letters, vol. 57, pp. 103-109, 2015.

 A. Heldring, J. M. Tamayo, C. Simon, E. Ubeda, and J. M. Rius, "Sparsified adaptive cross approximation algorithm for accelerated method of moments computations," IEEE Transactions on Antennas and Propagation, vol. 61, no. 1, pp. 240-246, 2013.

Chunbei Luo (iD), Yong Zhang, and Hai Lin (iD)

State Key Laboratory of CAD&CG, Zhejiang University, Hangzhou 310058, China

Correspondence should be addressed to Hai Lin; lin@cad.zju.edu.cn

Received 5 April 2018; Revised 11 June 2018; Accepted 24 June 2018; Published 22 July 2018

Academic Editor: Paolo Baccarelli

Caption: Figure 1: Periodic structure with a finite size of 4x5.

Caption: Figure 2: Radar cross sections (V-V polarization) of the 4x5 sphere array with (a) 0.5 m and (b) 0.1 m gap computed by FEKO, MLFMA, and periodic PSA. The incident plane wave angles satisfy [theta] = 0[degrees] and [phi] = 0[degrees]. The scattering angles work with [theta] varies from 0[degrees] to 360[degrees] and [phi] = 0[degrees].

Caption: Figure 3: Radar cross sections (V-V polarization) of the 32 x 32 sphere array with 0.5 m gap computed by FEKO, MLFMA, and periodic PSA. The incident plane wave angles satisfy [theta] = 0[degrees] and [phi] = 0[degrees]. The scattering angles work with [theta] varies from 0[degrees] to 360[degrees] and [phi] = 0[degrees].

Caption: Figure 4: Current density distribution of 10 x 10 patch array.

Caption: Figure 5: Radiation patterns of the 10 x 10 antenna array with 0.005 m gap computed by MLFMA and periodic PSA. The farfield angles work with [theta] varies from 0[degrees] to 360[degrees] and [phi] = 0[degrees].
```Table 1: Computation time and memory usage of MLFMA and
periodic PSA for scattering of 4 x 5 elements.

Gap       Methods      Computation time   Memory usage

0.5 m      MLFMA            210 s           8726 MB
Periodic PSA        109 s           2145 MB

0.1 m      MLFMA            222 s           8801 MB
Periodic PSA        131 s           2305 MB

Table 2: Computation time and memory usage of MLFMA and
periodic PSA for scattering of 32 x 32 elements.

Radius    Gap      Methods      Computation time   Memory usage

0.1 m    0.8 m      MLFMA            369 s           18229 MB
Periodic PSA         113s            707 MB
0.2 m    0.6 m      MLFMA            1551 s          71455 MB
Periodic PSA        895 s            3740MB
0.25 m   0.5 m      MLFMA            2341s           96906MB
Periodic PSA        2067 s           6136MB

Table 3: Computation time and memory usage of MLFMA and
periodic PSA for radiation of 10 x 10 antenna array.

Methods        Computation time   Memory usage

MLFMA               419 s           8841 MB
Periodic PSA        172 s            358 MB
```
COPYRIGHT 2018 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.