Efficient Compression of Far Field Matrices in Multipole Algorithms based on Spherical Harmonics and Radiating Modes

This paper proposes a compression of far field matrices in the fast multipole method and its multilevel extension for electromagnetic problems. The compression is based on a spherical harmonic representation of radiation patterns in conjunction with a radiating mode expression of the surface current. The method is applied to study near field effects and the far field of an antenna placed on a ship surface. Furthermore, the electromagnetic scattering of an electrically large plate is investigated. It is demonstrated, that the proposed technique leads to a significant memory saving, making multipole algorithms even more efficient without compromising the accuracy.


Introduction
The method of moments (MoM) [1] is well established for solving surface integral equations (SIEs) in the range of electromagnetic compatibility.Acceleration techniques like the fast multipole method (FMM) or the multilevel fast multipole algorithm (MLFMA) are often applied when analyzing electrically large objects to reduce both computation time and memory [2].In order to gain a time-and memory-efficient FMM implementation, various measures have been introduced in the past [3].This contribution proposes a compression of the far field matrices combining a spherical harmonic representation of the radiation pattern [4] with a radiating mode expression of the currents in a FMM group [5].
The paper is organized as follows: In Section 2 the fundamentals of MoM and FMM are briefly outlined.Section 3 describes the theoretical background of the proposed compression technique and in Section 4 two numerical examples are considered to investigate the impact of the developed approach on the required memory and the accuracy.

Formulation of MoM and FMM
The MoM can be applied to solve various types of SIEs.In this paper, it is focused on the electric field integral equation (EFIE) for perfect electric conducting (PEC) bodies.Discretizing the EFIE with MoM leads to a system of linear equations Here, Z ∈ C N ×N is the dense system matrix, I ∈ C N ×1 includes the unknown current amplitudes and V ∈ C N ×1 represents the excitation of the system.N is the number of unknowns.Applying FMM to MoM results in a blockwise approximated system matrix.The matrix block Z r,s couples a group s of basis functions and a group r of test functions.It is approximated as [5] with α α α r,s and R θ,φ r being the translation matrix and the disaggregation matrices, respectively.The columns of the aggregation matrices are radiation patterns of the N s basis functions in group s evaluated at K discrete points k = k 1 ... k K on the unit sphere (k-space formulation).For the i-th basis function of group s the radiation pattern is given by where d i is the distance between source point r and the center of group s.Ī denotes the unit dyadic and k is the normalized k-vector, pointing in radial direction on the unit sphere [2].

Compression of Far Field Matrices
In the following a compression technique for the aggregation matrices is described, which consist of two parts: The radiation pattern given by Equation ( 4) is represented by spherical harmonics and only radiating current modes are considered for the aggregation matrices given by Eq. (3).

Spherical Harmonic Representation of the Aggregation Matrices
The number of sampling points K is determined by the spectral content of the translation operator and the radiation patterns [2].This leads to an oversampling of radiation patterns in the k-space formulation [4].Thus, a more appropriate basis for the representation of f si k) should be applied.Since each radiation pattern is quasi-bandlimited (bandwidth P ) [2], it can be expressed by spherical harmonics Y m l as where c l,m si ∈ C 3×1 are weighting coefficients for the x-, yand z-component of f si (k).Following Eq. ( 5), the aggregation matrices of group s read with C s ∈ C M sph ×Ns .The discretized spherical harmonics are globally stored in Y ∈ C K×M sph and have to be computed only once.The number of spherical harmonics is M sph = (P + 1) 2 according to Eq. ( 5).This representation leads to a memory reduction, since only coefficients c l,m si have to be stored instead of an explicit k-space pattern [4].Note that F x,y,z s are expressed in Cartesian coordinates to avoid Gibb's phenomenon caused by discontinuities at the poles of F θ,φ s in spherical coordinates.The order of spherical harmonics P required for an approximation accuracy of 10 −α can be determined by [2]: where d denotes the half diameter of a group and k is the wavenumber.A box size of 0.5 λ as depicted in Figure 1 leads to d = √ 3/2 • 0.5 λ.Choosing an accuracy of 10 −4 results in P = 8.305 ≈ 8 for example.
In [4] a radiation pattern is computed in k-space and then transferred to spherical harmonics.Instead of doing that, we propose to determine the coefficients c l,m si directly from a spherical wave expansion for plane waves [6] e jk•di ≈ 4π where j l is the spherical Bessel function of order l and Ȳ denotes the conjugate of Y .Substituting Eq. (8) into Eq.( 4) leads to the spherical harmonic representation of Eq. ( 5).Hence, the corresponding weighting coefficients c l,m si can be computed via In the implementation, this integral is numerically evaluated by a quadrature rule.

Radiating Current Modes
To further reduce memory we exploit that only a few current modes contribute significantly to the radiation pattern of a group [5].These modes are determined by using a low-rank approximation of the coefficient matrices C x,y,z s according to: with A s ∈ C M sph ×qs and B s ∈ C qs×Ns .To achieve such a compression, a truncated singular value decomposition (TSVD) [7] is frequently used.For a matrix approximation with accuracy , only q s singular values larger than the accuracy threshold have to be considered.A mathematical description of TSVD can be found in [7], [8].The singular values for the coefficients matrices C These matrices can be interpreted as follows: B s transfers the considered basis functions into q s radiating modes and A s contains the weighting coefficients for the spherical harmonics of the far field pattern associated with these radiating modes.Fig. 2 illustrates the radiating current modes corresponding to the three largest singular values of matrix C x 1 .These modes are defined by the first three rows of B x 1 .

Memory and Computational Complexity
The number of elements needed to be stored for the aggregation matrices of one group with N s basis functions is shown in Table 1.Here, three compression techniques are compared with a conventional approach: A TSVD of the matrices F θ,φ as presented in [5], a compression with spherical harmonics [4] and the proposed technique.Note that the memory required to hold Y is neglected, because it does not depend on the total number of unknowns.The factor 2 takes into account that θ and φ components have to be stored for the conventional approach and for the compression with TSVD.Factor 3 reflects the use of Cartesian   1, the proposed method is superior to a pure spherical harmonic representation if the condition is fulfilled.Table 2 illustrates the elements to be stored for the considered example structure depicted in Fig. 1.In this example a TSVD threshold of = 10 −4 is applied, which leads to q s = 26 radiating modes.q s is usually different for each matrix C x , C y and C z as shown in Fig. 1.For the sake of simplicity q s is assumed to be the same for all matrices.Furthermore P = 8 has been chosen to achieve an accuracy of 10 −4 in accordance to Eq. ( 7); this leads to M sph = 81.The group under consideration includes N s = 95 basis functions and the number of points discretizing the unit sphere is K = 338.
Table 2: Number of elements for the aggregation matrices associated with the example structure in Fig. 1.

Method Elements Factor Conventional
64 220 -TSVD [5] 22 516 2.85 Spherical Harmonics [4] 23 093 2.78 Proposed Method 13 728 4.68 Fig. 3 illustrates the impact of various parameters P (order of spherical harmonics) and (TSVD threshold) on the compression factor and the relative errror of matrix F θ for the considered example.The compression factor is defined as ratio between memory for conventional approach and memory for proposed method.It is shown that the compression factor behaves inversely proportional to the accuracy, hence a trade-off between memory saving and accuracy must be made.Moreover one observes that for each pair of parameters the accuracy is bounded by either P or .Since compressing the aggregation matrices will influence the computational effort, a short complexity analysis is given below.The main steps which affect the computation time are: A TSVD of the coefficient matrices C x,y,z s for each group s and the multiplication of a vector with the aggregation matrices during the iterative solution procedure of the entire problem.In particular follows: 1.The complexity of a TSVD for one group is O(N 3 s ) if N s ≈ q s [8].This step has to be executed for M groups.N s and q s do not depend on the total number of basis functions N , but the number of groups M is proportional to N for a fixed group size.Thus the complexity for this stage is O(N ).

The multiplication with the compressed aggregation
matrix in the iterative solution procedure is accelerated, since the elements of these matrices have been reduced by the proposed technique.The complexity for this stage is O(N ).
It can be summarized that, although the first point will slow down the set-up phase and the second point will speed up the solution procedure, the complexity of multipole algorithms -i.e. O(N logN ) for MLFMA -remains unchanged when using the proposed method.the associated coefficient matrices C x s are drawn in Fig. 5.For all three geometries the magnitude of singular values is falling.The slope decreases with increasing geometrical complexity.It can be concluded, that the number of current modes q s becomes larger with more complex structures.This corresponds to investigations carried out in [5].
To obtain an accuracy of 10 −4 for the specified examples, q s = 26, q s = 43 and q s = 77 modes are required for structure a), b) and c), respectively, see Fig. 5.
Next, the proposed technique is investigated for varying mesh densities.Fig. 6 illustrates the compression rate as a function of mesh density starting with 5 elements per wavelength and ending with 40 elements per wavelength.Spherical Harmonics [4] Upper Bound In all cases, the compression rate grows with increasing mesh density.This behavior can be explained by the number of radiating current modes q s , depends the shape only as indicated by Fig. 5. Accordingly, the compression rate is large for less complex geometries.As q s is a fixed number for a specified geometry, an upper bound for the compression factor can be derived: Here, the upper bound is 8.67, 5.24 and 2.93 considering geometry a), b) and c), respectively for K = 338.These bounds are approached with increasing mesh density as shown in Fig. 6.In case of examples a) and b), the proposed method is superior to a pure spherical harmonic compression for mesh densities higher than 8 elements per wavelength.For example c) no significant memory saving compared to a pure spherical harmonic representation is at-tained, because Eq. ( 12) is not fulfilled.Thus, the method is expected to perform best for objects containing flat structure parts with high mesh densities.

Numerical Results
Results of the proposed compression technique are presented in the following.The technique was applied to a MoM Code [9] with FMM and MLFMA acceleration, respectively.All computations have been performed in double precision on a PC with a AMD Opteron 6140 2.6 GHz Processor.The transpose-free quasi-minimal residual (TFQMR) algorithm [10] was used to solve the equation system iteratively.The first example to be considered is a ship model with a total length of 120 m as depicted in Fig. 7 and Fig. 8.The discretized surface involves 126 962 Rao-Wilton-Glisson (RWG) basis functions [11].The ship is excited by a monopole antenna at a frequency of 72 MHz.
The compression was carried out using the techniques described by [4], [5] and the proposed method for different parameters P (order of spherical harmonics) and (accuracy of TSVD).In approach [4] only spherical harmonics are used with parameter P ; in approach [5] only TSVD is used with parameter .Each method has been applied to the FMM with a box size of 0.5 λ and MLFMA with a box size of 0.25 λ on the finest level.In all cases Galerkin testing with R = F H was used.Here, matrix R is not computed explicitly.Table 3 illustrates the memory reduction for the 289.5 2.17 aggregation matrices.The proposed technique has a better performance than the techniques suggested in [4] and [5] because it combines both aspects, i. e. spherical harmonics and radiating modes.It has to be emphasized, that the approach under consideration requires slightly more memory than the spherical harmonic representation, choosing P = 4 and = 10 −3 .This is because the requirement of Eq. ( 12) has been violated for some groups of basis functions.In this first example the total memory for all required FMM matrices is reduced by 27% from 4.24 GB to 3.09 GB.The total memory for MLFMA is decreased from 1.30 GB to 0.82 GB, which corresponds to a 37% reduction.Fig. 7 shows the electric field distribution along an observation path 1 m above the deck.Only small deviations are caused by the proposed technique.In this example we aimed at 10 −3 accuracy with P = 8, = 10 −3 and 10 −2 accuracy with P = 5, = 10 −2 following the analysis given in Section 3 .The maximum relative error of the computed field is below 10 −2 and 10 −3 , respectively, which confirms an error controllability for the proposed method.Furthermore the radiation pattern of the antenna is given in Fig. 8, where the deviation to a conventional approach is hardly noticeable.The time for compressing the far field  The second example is an electrically large PEC plate excited by a linearly polarized plane wave as illustrated in Fig. 9. Two different cases are investigated: The square plate with an edge length of a = 32 λ and a = 64 λ, respectively.In both cases the plate is discretized with 1 085 544 basis functions.Computational results for a conventional approach, the proposed technique with P = 4 and = 0.01 as well as a pure spherical harmonic representation with P = 4 are given in Table 4.The memory  required for the aggregation matrices with the introduced technique is 0.63 GB in case of a = 32 λ and 1.07 GB in case of a = 64 λ compared to 5.24 GB for the conventional approach.For a = 64 λ, the total memory of all MLFMA matrices is decreased from 7.89 GB to 3.72 GB, which relates to 53%.The compression time was 78 s.Solving this problem with a residual tolerance res = 3 • 10 −3 took 6.9 h on a single CPU.Table 4 points out, that the compression rate grows with increasing mesh density, as predicted in Section 3.For a = 32 λ the average number of basis functions in one group was 66.3, for a = 64 λ the average was 16.6 basis functions per group.To validate the accuracy of our approach, MLFMA is compared to an analytical estimation for the plate with a = 64 λ.Based on assumptions known from physical optics, the surface current density is approximated by J ≈ 2 • n × H inc .From this the bistatic radar cross section (RCS) can analytically be computed using Green's function for the far field [12].Since the plate is large in terms of wavelength, accurate results are expected

Conclusions
An efficient compression of far field matrices in multipole methods is proposed.The compression technique has been applied to the conventional FMM and its multilevel extension MLFMA, leading to a significant memory reduction up to 53%.It turns out that the developed method is superior to techniques known from literature.In addition to this the proposed approach is shown to be very accurate.

Figure 1 :
Figure 1: Singular values of the coefficient matrices C x 1 , C y 1 and C z 1 for an example structure inside a group with edge length of 0.5 λ.The box includes a flat plate as an example for the content of boxes in the entire MLFMA.

Figure 2 :
Figure 2: Current modes in a group corresponding to the first, second and third singular value of matrix C x 1 for the example structure given in Fig 1.

Figure 3 :
Figure 3: Top: Relative error of the aggregation matrix F θ estimated in Frobenius norm [8].Bottom: Compression factor for aggregation matrix (compared to a noncompressed matrix).P : Maximum order of spherical harmonics.: Accuracy of compression with TSVD.Example: Flat plate included in a group with 0.5 λ edge length as depicted in Fig. 1.

Figure 5 :
Figure 5: Singular values of the coefficient matrix C x given in Eq. (10) for three different geometries a) square plate, b) dihedral corner and c) cube as depicted in Fig. 4.

Figure 6 :
Figure 6: Compression factor versus mesh density for three different geometries: a) square plate, b) dihedral corner and c) cube, see Fig. 4. Dashed line: upper bounds for the compression rate given by Eq. (13).

Table 3 :FMM
Memory for the aggregation matrices with the proposed technique using FMM (box size 0.5 λ) and MLFMA (box size 0.25 λ).Example: ship surface with 126 962 unknowns.

Figure 7 :
Figure 7: Top: Max.magnitude of E field along the indicated observation path (length 21 m, 1 m above the deck), computed with a conventional FMM approach and the proposed compression for aggregation matrices.P : Maximum order of spherical harmonics.: Accuracy of compression with TSVD.Bottom: Relative error of E field (referring to the conventional approach).

Figure 8 :
Figure 8: Radiation pattern in the upper hemisphere due to an monopole antenna placed on the ship surface.Comparison of proposed approach with conventional method.P : Maximum order of spherical harmonics.: Accuracy of compression with TSVD.

Table 4 :
Memory for the aggregation matrices with different compression techniques using MLFMA (box size: 0.25 λ).Example: square plate with 1 085 544 unknowns for the edge lengths a = 32 λ and a = 64 λ.MLFMA (0.25 λ) Mem.(GB) Factor wave around θ = 0 • .Fig. 10 compares MLFMA results with the analytical estimation.Here, the bistatic radar cross section for θ = 0 • ... 90 • and φ = 90 • is given.A very good agreement can be observed in the range between θ = 0 • and θ = 10 • .With increasing angle the deviations become more pronounced due to the optical approximation for the analytical computation.In order to have a more confident analysis, the results of a conventional MLFMA computation are given.A very good agreement between both MLFMA variants can be observed and only in the range θ = 80 • ... 90 • marginal deviations are noticeable.

Figure 10 :
Figure 10: Bistatic RCS of a square PEC plate (a = 64 λ) with 1 085 544 unknowns, computed with the a conventional approach, the proposed technique and an analytical estimation.P : Maximum order of spherical harmonics.: Accuracy of compression with TSVD.

Table 1 :
Number of elements for aggregation matrices of one group with different compression methods.
s K + N s Spherical Harmonics [4] 3M sph N s Proposed Method 3q s M sph + N s components for the other techniques.As implied by Table