On the Application of the Incomplete QR Algorithm to the Analysis of Microstrip Antennas

In this paper, we provide some insight into the usage of fast, iterative, method-of-moments (MoM) solution of integral equations (IE) describing antennas and other metallic structures immersed in a planar multilayered environment. Based on the form of multilayered media Green’s functions, we extract free-space terms, associated with direct rays within the analyzed structure, reducing the number of significant interactions required to describe the rest of MoM matrix. Next, we show that it is possible to construct a hybrid algorithm, where the fast multipole method (FMM) is used to the free-space matrix part, while the reduced rank incomplete QR (iQR) decomposition is applied to the remaining portion of the MoM matrix. This HM-iQR (hybrid multipole incomplete QR) method is applied to a relatively large (in terms o f the number of unknowns) problem of plane wave scattering by a finite array of rectangular microstrip patches printed on a grounded dielectric slab. Computation results from the new algorithm are compared to literature data and to the results of the pure low rank IE-QR method.


Introduction
Analysis of metal structures placed in a planar multilayered environment has been of interest to the community of researchers dealing with computational electromagnetics for decades.Although the integral equations/method-ofmoments approach to that problem is well established [1], the issue of long computation times for electrically large objects is still relevant.The problem, as in the free-space situation, is associated with the large number of unknowns N in the final linear equation set.
In free-space problems, limitations due to the computational complexity of order O(N 3 ), when applying direct methods like LU decomposition, are usually overcome by us-ing Krylov-subspace iterative methods with techniques like fast multipole method (FMM) [2] or low rank IE-QR algorithm [3] applied to speed-up computation of matrix-vector product.Both techniques allow for reducing the cost of a single iteration from O(N 2 ) to O(N 3/2 ) in the case of the single-level method, or to O(N log N ) in the case of its multilevel counterpart, in the case of FMM called multilevel fast multipole algorithm (MFLMA) (cf.[2,3]).
In multilayered media, the application of the FMM is limited, as the algorithm is based on the clever representation of the free-space Green's function [2].One solution is to represent multilayered Green's function with set freespace-like functions, which is called discrete complex image method (DCIM) [2,4].However, for structures with large lateral dimensions, the physical phenomena are often dominated by surface waves, for which discrete images are not well-suited [5].On the other hand, the low rank IE-QR algorithm (in this paper called the incomplete QR -iQR) is based on purely algebraic MoM matrix operations, and -as suchdoes not depend on the form of Green's functions.However, as shown below, for some configurations it can be made faster by applying the hybridization with FMM.
The rest of this paper is organized as follows: first we formulate the representative problem concerning a finite array of rectangular microstrip patches printed on a grounded dielectric slab.Then, we outline the conventional method of analysis, as well as FMM and iQR algorithms, and the idea of the hybrid method is proposed, which enables analysis of the properties of the iQR when applied to planar stratified media problems.Finally, we show computational examples validating the proposed idea and making it possible to draw some general conclusions.
It is to be stated that by combining the two methods with computational complexity of (in the case of the single level implementation) order O(N 3/2 ), we do not expect to obtain algorithm with lower complexity.Instead, we expect to reduce the constant factor standing in front of N

Formulation
Let's consider a finite array of microstrip patches depicted in Fig. 1.The patches are made of the perfect electric conductor (PEC).The structure is illuminated by a plane wave, coming from the specified direction from the upper half-space, which is filled with air.The quantity to be determined is the radar cross section (RCS).
The total electric field at any point of space is the superposition of the incident and the scattered fields.On the surface of the PEC patches, it satisfies the boundary condition: (E s (J) where E i indicates the incident field, while sub-script t stands for tangential field component.The scattered field E s is generated by the surface current on the patches: with the magnetic vector potential A and scalar potential Φ defined as: In ( 3) and (4) S represents total surface of the patches.Vectors r i r denote observation and source points, respectively, µ 0 and 0 are permeability and permittivity of free space, J denotes the surface current, and K A (dyadic) and K Φ (scalar) are the Green's functions [1], which can be expressed as Sommerfeld integrals of their spectral domain counterparts KA and KΦ .For the problem at hand, the currents are assumed planar, and also only horizontal field components are taken into account in (1).Thus, only GA vv and KΦ from the whole set described in [1] are necessary, which are expressed by means of so-called transmission line Green's functions (TLGFs) V h i and V e i , as [1]: TLGFs V p i for the case depicted in Fig. 1, expressed in the notation of [1], have the general form: The first term in (7) represents the direct ray between the source and the field points, while the second term represents the rays that undergo partial reflections at the upper and lower slab boundaries before reaching the observation point [1].The details description of all coefficients used in (7) can be found in [1].Here, we only note that all the coefficients depend on the configuration and properties of the multilayered environment, frequency of interest and z coordinates of source and observation points.The index n refers to the section of the equivalent transmission line that represents layered medium in the spectral domain (cf.[1]).As the microstrip patches lie on the interface separating different media, we have a freedom of choice of the section to which the interface belongs -total formulas remain the same irrespective of that choice.So, let's assume that the interface (which means both z and z ) belongs to the upper half-space filled with air, with which we associate index n = 0.Then, substituting (7) into (5), we get: where ρ , k 0 is the free-space wavenumber, and F A (k ρ , z, z ) comes from the second term in (7), containing the summation.Substituting ( 7) into (6), we can transform KΦ into the similar form: where again F Φ (k ρ , z, z ) corresponds to phenomena associated with layered structure of the medium under investigation.
In the spacial domain, the Green's functions become [1]: where Above, J n is Bessel function of order n, while ρ denotes distance between the source and observation points.Note that in the above formulas, we dropped out the z dependence, as all source and observation points are confined to the single value of z.
Putting ( 8), ( 9) into ( 10) and (11), respectively, and using Sommerfeld identity [6]: where R = |r − r | denotes distance between source (r ) and observation (r) points, we finally get: where we bear in mind that for the analyzed case z = z , so R present in ( 13) is equal to the lateral distance ρ.Details regarding integration of Sommerfeld integrals may be found in [7].
Please note that when the multilayered medium is replaced by free-space, the corresponding transmission lines in the spectral domain become infinite ones with constant parameters.Then, F A (k ρ , z, z ) and F Φ (k ρ , z, z ) terms disappear, and the above formulas reduce to their well known free-space counterparts.
After suitable triangulation of S, the application of RWG basis functions [8], and after the usual testing procedure, equation ( 1) is transformed into the matrix formula: where Z is a dense N × N impedance matrix, while I and V are column vectors of length N describing, respectively, the approximation of the current distribution and the (tested) incident field on S.
Taking into account the form of ( 14) and ( 15), we can also conclude that in the general, multilayer case, the matrix elements may be expressed as the sum of two terms, one coming from the direct free-space rays, and the second associated with the presence of the multilayered structure.So, also the MoM matrix can be written as the sum of two matrices: where Z 0 is the "free-space" MoM matrix, while the elements of Z m (m for "multilayered") are computed with the use of numerical integration of Sommerfeld integrals associated with F A and F Φ .Apart from the way, how MoM matrix entries are computed, putting Z into the form of ( 17) has also consequences for the solution of the linear system, which will become apparent in the subsequent discussion.
Finally, when ( 16) is solved, the obtained approximation of the current allows to compute resulting scattered fields at arbitrary observation point.In the subsequent examples we use far-field approximate formulas for multilayared medium [9].

FMM and Incomplete QR
As outlined in the Introduction, when an iterative method is applied to solve (16), either FMM or incomplete QR may be used to speed-up the associated matrix-vector multiply, depending on the type of the problem.
The fast multipole method is a numerical algorithm proposed by Greengard and Rokhlin for solving Helmholtz and Laplace problems [2,10].Nowadays, FMM and its multilevel extension MLFMA are applied to seek rapid solution to integral equations of scattering.The detailed description of FMM can be found for example in [10,11], so we will omit here the detailed description.
In contrast to the above, some details about the incomplete QR method are required to follow the idea shown in the paper, so we decided to include here a short outline.In the method, the first step is to divide the N basis functions into G localized groups, each supporting M = N/G basis functions.The good choice is G ∼ N 0.5 .The significant information between two groups is limited if the distance between them is greater than their physical size [3].Then the number of independent vectors needed for expression of interactions between these two groups is much smaller than number of basis functions.The incomplete QR algorithm consists of four steps: column index selection, updating the Q matrix, row index selection, and computing the R matrix.
The algorithm begins with the arbitrary selection of a receiver from the receiver group.With respect to this receiver, the most dominant transmitter is chosen from the transmitter group.The most dominant transmitter will then in turn determine the most dominant receiver, which is not the duplicate of the chosen one.With these two chosen receivers, the second most dominant transmitter is selected.This process continues until the remaining receivers and transmitters contribute no significant information [12].In this process Modified Gram-Schmidt (MGS) algorithm is applied to determine the index of the column vector in Z, that is the most linearly independent on the previously selected columns.In this case local impedance matrix Z m×n i j is factorized into Q m×r i j and R r×n i j , where i and j are receiving and transmitting group, respectively, and m and n are the number of basis functions, and r is the approximate rank of Z m×n i j [3].This can be written as: Here z p is the column vector due to the pth basis function in the transmitting group, and v q is the row vector due to the qth basis function in the receiving group [13].This low rank IE-QR factorization reduces the complexity matrixvector product from O(mn) to O(r (n + m)), where r is much less than m and n for well separated groups.This leads to speed up solution time for general electromagnetic radiation and scattering problems.Detailed description of incomplete QR algorithm can be found in [3].
Note that in both FMM and iQR algorithms, the speedup is associated with well-separated groups.Self interactions, and interactions with neighboring groups are computed as in the standard MoM solutions.

The Hybrid Method
Taking into account (17), we can see that each matrixvector multiplication in the iterative algorithm may be computed via separate multiplications of the two matrices involved, namely Z 0 and Z m .To speed-up the first multiplication, we can use FMM, which in free space can be faster than incomplete QR [14].For the second multiplication, the iQR is applied.Here however, this algorithm should run faster than in the case "pure" iQR method, i.e. the method without extraction of Z 0 , as the number of significant relations between receivers and transmitters is expected to be smaller.Also, when distances between groups increase, the number of significant relations becomes smaller [15].We check those expectations in the next section.

Numerical Results
In order to check the correctness and efficiency of the algorithms, they were implemented in Fortran, and run on the Intel®Xeon®X5472 3.00 GHz workstation with 64 GB of RAM.No parallelization was used, so the computations were performed on a single core.
As computational examples, we analyzed antenna arrays described in [16], and shown in Fig. 1.Different numbers of patches (3 × 3, 7 × 7, and 11 × 11) were taken into account.The quantity of interest was monostatic radar cross section (RCS) for a θ polarized plane wave incident from broadside.The results obtained by iQR and HM-iQR methods are compared with literature data [16,17].
The results (RCS vs. frequency) are shown in Figs. 2, 3 and 4. It can be seen that both iQR and HM-iQR methods give results consistent with literature data, so the methods are correct.Let's take a look at the number of significant interactions in off-diagonal MoM matrix blocks in both iQR and HM-iQR (Z m portion).In our case, the groups defined  within the algorithms are associated with patches with the same number of basis/testing functions, so the MoM matrix blocks corresponding to group-group interactions are all square, and of the same size.However, the number r in (18) is different, depending on the mutual position of the two interacting patches.Therefore, we decided to give the mean value of r, which may be used as the overall measure of the iQR efficiency (the lower r, the better).The values of r are given in Tab. 1.We can see that r, and, what follows, the number of multiplications performed in the iterative algorithm, for HM-iQR is about 1/3 lower than for "pure" iQR.QR decomposition is done for well separated groups, e.g. for groups 1 and 49 in 7 × 7 case with 225 basis functions in one group.In this example for HM-iQR in matrix Z m minimum and maximum values are 1.33 • 10 −9 and 3.85 • 10 −7 , respectively.For "pure" iQR minimum and maximum values are 6.47 • 10 −12 and 1.31 • 10 −8 , respectively.We noticed also, that in this case r is 14 and 36, for HM-iQR and "pure" iQR, respectively.Those values of r are below average results presented in Tab. 1, and it is expected, because groups 1 and 49 are extremely far away from each other.
Comparison of the times (Tab. 2 and in Tab. 3) shows, that HM-iQR method may be, depending on an example, almost 30% faster.This is visible and valuable especially in calculation of RCS not only for single frequency, but for a range of frequencies (in this paper for 21 frequency points).All calculations in HM-iQR and iQR were done with tolerance 10 −10 .For case 7 × 7 the proper result was achieved slower than in others examples.If the tolerance was changed to 9.8 −10 (in both methods), the HM-iQR was about 30% faster and gave correct result.It means, that algorithm is sensitive to tolerance defined in conjugate gradient method.
During numerical tests, we have also found that when increasing the distance between patches beyond a certain limit, the algorithm loses its beneficial properties.This may be attributed to the fact that for far-lying patches the coupling mechanism is dominated by surface waves phenomena, so extracting the FMM part does not reduce r in the iQR, while adding unnecessary overhead into the computation time.This issue should however be studied in more detail in future investigations.

Conclusions
The paper presents an enhancement into MoM computations in layered media problems concerning electrically large objects.The Hybrid Multipole -incomplete QR algorithm has been applied to determine RCS of a finite array of rectangular microstrip patches printed on a grounded dielectric slab.The described HM-iQR algorithm allowed to get accurate and rigorous solution to the given problem and required less computation time than the pure iQR method.Although the speed-up is not spectacular, the method allows to get some additional insight into correlation between physical phenomena and numerical methods.
3/2 , and, what is more important, to indicate what physical mechanisms are relevant, when applying iQR to multilayered problems.Throughout the paper the e jωt time convention is used.