# On the Fast Direct Solution of a Preconditioned Electromagnetic Integral Equation

Davide Consoli, Clément Henry, Alexandre Dély, Lyes Rahmouni, John Erik Ortiz Guzman, Tiffany L. Chhim, Simon B. Adrian, Adrien Merlini, Francesco P. Andriulli

► To cite this version:

Davide Consoli, Clément Henry, Alexandre Dély, Lyes Rahmouni, John Erik Ortiz Guzman, et al. On the Fast Direct Solution of a Preconditioned Electromagnetic Integral Equation. International Conference on Electromagnetics in Advanced Applications (ICEAA 2022), Sep 2022, Le Cap, South Africa. 10.1109/ICEAA49419.2022.9899870. hal-03762873

HAL Id: hal-03762873

https://hal.science/hal-03762873v1

Submitted on 24 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# On the Fast Direct Solution of a Preconditioned Electromagnetic Integral Equation

Davide Consoli¹, Clément Henry², Alexandre Dély¹, Lyes Rahmouni¹, John Erik Ortiz Guzman³, Tiffany L. Chhim¹, Simon B. Adrian⁴, Adrien Merlini², and Francesco P. Andriulli¹

¹Politecnico di Torino, Turin, Italy

²IMT Atlantique, Brest, France

³Universidad de Nariño, Pasto, Colombia

⁴Universität Rostock, Rostock, Germany

**Abstract**—This work presents a fast direct solver strategy for electromagnetic integral equations in the high-frequency regime. The new scheme relies on a suitably preconditioned combined field formulation and results in a single skeleton form plus identity equation. This is obtained after a regularization of the elliptic spectrum through the extraction of a suitably chosen equivalent circulant problem. The inverse of the system matrix is then obtained by leveraging the Woodbury matrix identity, the low-rank representation of the extracted part of the operator, and fast circulant algebra yielding a scheme with a favorable complexity and suitable for the solution of multiple right-hand sides. Theoretical considerations are accompanied by numerical results both of which are confirming and showing the practical relevance of the newly developed scheme.

**Index Terms**—Integral equations, preconditioning, fast direct solvers.

## I. INTRODUCTION

Integral equation strategies are effective for solving scattering and radiation problems over a wide range of frequencies since they do not require absorbing boundary conditions, they are free from numerical dispersion, and they only require the discretization of the surfaces of the scatterers and the radiators. Unfortunately, matrices arising from their discretization are dense and have dimensions that grow for increasing frequencies. Naive solutions would scale cubically with the number of unknowns and easily jeopardize the advantages of integral equation approaches. Fast matrix-vector multiplication algorithms [1] can substantially decrease this computational burden, but they require iterative solutions, which renders them less attractive when many right-hand sides (RHSs) must be considered. An alternative is the use of fast direct solvers that are building, in reduced complexity, directly the inverse of the system matrix and are efficient for solving problems with multiple RHSs. There are several effective strategies for direct solutions that often rely on hierarchical decompositions both for problems in the low- and in the high-frequency regimes (see [2]–[4] and references therein).

This work has been funded in part by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (ERC project 321, grant No.724846), in part by the Italian Ministry of University and Research within the Program PRIN2017, EMVISIONING, under Grant 2017HZJXSZ, CUP:E64I190025300, and in part by the ANR Labex CominLabs under the project “CYCLE”.

In this work, we leverage a different approach. We first precondition the combined field integral equation (CFIE) in a suitable way [5]. We then exploit the particular structure of the preconditioned spectrum to separate its principal part contribution and compress its remainder. This is made practical by extracting an equivalent circulant problem, which has the advantage of automatically extracting the principal part contribution and the remainder is compressed in favorable complexity. This results in a non-hierarchical scheme in which a single skeleton form is required and the solution for several RHSs can be obtained efficiently. Theoretical considerations are matched with numerical results confirming and showing the practical relevance of the new scheme.

## II. NOTATION AND BACKGROUND

Consider the scattering of a TE time-harmonic incident electromagnetic field ($E^{\text{inc}}$, $H^{\text{inc}}$) of angular frequency $\omega$ impinging on a perfectly electrically conducting (PEC) body, modeled by a 2-dimensional convex contour $\gamma$ that resides in a medium of permittivity $\epsilon$, permeability $\mu$, impedance $\eta = (\mu/\epsilon)^{1/2}$, and corresponding wavenumber $k = \omega\sqrt{\mu\epsilon}$. The combined field integral equation (CFIE) that relates the induced tangential current density $j$ on $\gamma$ to the incident electric and magnetic fields reads

$$ \frac{1}{ik} \mathcal{N}^k(j_t) + \left(\frac{\mathcal{I}}{2} - \mathcal{D}^k\right) (j_t) = -\frac{1}{\eta} E_t^{\text{inc}} - H_z^{\text{inc}} \quad (1) $$

with

$$ (\mathcal{N}^k j_t)(r) := -\frac{\partial}{\partial n} \int_{\gamma} \frac{\partial}{\partial n'} g_k(\mathbf{r}, \mathbf{r}') j_t(r') dr', \quad (2) $$

$$ (\mathcal{D}^k j_t)(r) := \int_{\gamma} \frac{\partial}{\partial n'} g_k(\mathbf{r}, \mathbf{r}') j_t(r') dr', \quad (3) $$

and $g_k(\mathbf{r}, \mathbf{r}') = i/4H_0^{(1)}(k\|\mathbf{r} - \mathbf{r}'\|)$. After boundary element discretization with piecewise linear functions $\lambda_i(\mathbf{r})$, (1) is discretized as

$$ \frac{1}{ik} \mathbf{N}^k \mathbf{j} + \left( \frac{\mathbf{G}}{2} - \mathbf{D}^k \right) \mathbf{j} = -\frac{1}{\eta} \mathbf{e} - \mathbf{h} \quad (4) $$

with $[\mathbf{e}]_i = \int_{\gamma} \lambda_i(\mathbf{r}) E_t^{\text{inc}}(\mathbf{r}) d\mathbf{r}, [\mathbf{h}]_i = \int_{\gamma} \lambda_i(\mathbf{r}) H_z^{\text{inc}}(\mathbf{r}) d\mathbf{r}, [\mathbf{N}^k]_{ij} = \langle \lambda_i, \mathcal{N}^k \lambda_j \rangle, [\mathbf{D}^k]_{ij} = \langle \lambda_i, \mathcal{D}^k \lambda_j \rangle, [\mathbf{G}]_{ij} = \langle \lambda_i, \lambda_j \rangle$, with $\langle f, g \rangle = \int_{\gamma} f(\mathbf{r}) g(\mathbf{r}) d\mathbf{r}$.
<figure><img src="image_2.png" /><figcaption>Fig. 1: Spectrum of the equation of this work (continuous counterpart of (5)) obtained for a cylinder.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2: Continuous spectrum of the 2D TE-MFIE obtained for a cylinder.</figcaption></figure>

### III. A SUITABLE INTEGRAL EQUATION AND ITS SPECTRAL ANALYSIS

A well conditioned equation can be obtained from (1) by leveraging a modified version of the preconditioner proposed in [5]

$$ \left( \mathbf{S}^{\tilde{k}} \mathbf{G}^{-1} \mathbf{N}^{k} + \left( \frac{\mathbf{G}}{2} + \mathbf{D}^{\tilde{k}} \right) \mathbf{G}^{-1} \left( \frac{\mathbf{G}}{2} - \mathbf{D}^{k} \right) \right) \mathbf{j} = - \frac{i k}{\eta} \mathbf{S}^{\tilde{k}} \mathbf{G}^{-1} \mathbf{e} - \left( \frac{\mathbf{G}}{2} + \mathbf{D}^{\tilde{k}} \right) \mathbf{G}^{-1} \mathbf{h}, \quad (5) $$

where $\tilde{k} := k + 0.4ik^{1/3}a^{-2/3}$ following [6], [7] with a suitable average of the radii of curvature along $\gamma$. When comparing the spectrum of the above equation (Fig. 1) with that of the standard formulations such as the EFIE, the MFIE (Fig. 2), or even the Calderón-preconditioned CFIE, one can see that the spectral content of (5) clusters, and is maximal, around the surface resonant point attained when the spatial frequency equals $ka$. Note that for ease of reading the singular values are ordered by the Fourier content of their singular vectors. This suggests that selecting the singular vectors corresponding to this region of maximal spectral strength and filtering out the deviation from the halved identity of the others would lead to a compression of the electromagnetic operator in the high frequency regime—in which the discretization parameter is kept at a fixed proportion of the wavelength. In particular, this would correspond to searching the solution in the space of corrections to the RHS in (5) of the form $s \mapsto M(s)e^{iks}$, with $M(s)$ slowly varying.

### IV. HANDLING THE ERROR ON THE ELLIPTIC SPECTRUM

Although eigenvector extraction (after proper shifting and symmetrization) is an option for obtaining a direct solver out of (5), at higher precision this may become less practical and skeleton form obtaining algorithms can be used instead. However, the discretization error in the elliptic part of the spectrum can interfere with this approach. The discretization will cause a constant ($O(1)$) error in the elliptic spectrum for first kind preconditioned operators (see Fig. 3). There are

<figure><img src="image_4.png" /><figcaption>Fig. 3: Spectrum of the matrix C for increasing frequency, with a compression parameter $\epsilon = 0.015$. The green markers indicate which singular values/singular vectors are selected by the compression scheme.</figcaption></figure>

several strategies to overcome this problem. Here we opted for subtracting an equivalent circular problem from the operators defined on $\gamma$. This will, in practice, extract the second kind part of the original operators, leaving only a compact operator for which the relative error (with respect to the identity) will be decreasing.

### V. THE DIRECT SOLVER

To exploit these considerations to the fullest, we will study the deviation of the electromagnetic operators defined on $\gamma$ from their equivalents defined on the circle $\gamma_c$ that has the same perimeter as $\gamma$. These auxiliary matrices, denoted by the subscript $c$, are obtained with the same uniform discretization in the curvilinear abscissa as the main matrices defined on $\gamma$—both of which are discretized with piecewise linear basis functions in the curvilinear abscissa. Note that under these conditions we have $\mathbf{G} = \mathbf{G}_c$, among other critical properties. Then by defining

$$ \mathbf{C} := \left( \left( \frac{\mathbf{G}}{2} + \mathbf{D}^{\tilde{k}} \right) \mathbf{G}^{-1} \left( \frac{\mathbf{G}}{2} - \mathbf{D}^{k} \right) + \mathbf{S}^{\tilde{k}} \mathbf{G}^{-1} \mathbf{N}^{k} \right) \quad (6) $$
and its circle counterpart $C_c$, we compute the skeleton form of the matrix

$$C - C_c := UV^T, \quad (7)$$

which can be obtained efficiently based on the expected low rank of the difference and on the fact that all operations involving the matrices and operators on $\gamma_c$ are done rapidly via fast Fourier transform (FFT). Finally, the solution of the original problem can be obtained efficiently, for any number of RHSs, by inverting

$$C = C_c (I + (C_c)^{-1} UV^T), \quad (8)$$

where $I$ is the identity matrix, using the Woodbury matrix identity [8]

$$C^{-1} = C_c^{-1} - C_c^{-1}U (I + V^T C_c^{-1} U)^{-1} V^T C_c^{-1} \quad (9)$$

that can be computed efficiently given the low-rank (i.e., the inner dimension) nature of the skeleton form and circulant algebra.

## VI. NUMERICAL RESULTS

The direct solution strategy of this work has been tested on an ellipse with semi-major axis 2 m and semi-minor axis 1 m. A close inspection of the spectrum of $C$ (Fig. 3) reveals that the spectral corruption caused by the discretization error poses a real challenge for compression purposes, since they would artificially increase the rank of the compression if the looked-for accuracy was too high (the limit case is indicated by the red lines in the figure). On the contrary, the elliptic spectrum of $(C-C_c)$ (Fig. 4) does not remain constant, but it decreases for the high spatial frequencies, which makes it an ideal candidate for applying the proposed compression scheme. Fig. 5 shows the behavior of the inner rank of the skeleton in (7) that shows a growth not larger than $k^{1/3}$ (which is a bound that can be theoretically predicted for the circular case). This provides in this case an upper bound for the computational complexity not larger than $O(N^{4/3})$ where in practice lower complexities can be obtained given that for a constant truncation error on the matrix the accuracy of the solution increases with frequency, as can be seen in Fig. 5. Thus, the effectiveness of the scheme is ensured in the high-frequency regime.

## REFERENCES

[1] W. C. Chew, J.-M. Jin, E. Michielssen, and J. M. Song, Eds., *Fast and Efficient Algorithms in Computational Electromagnetics*. Artech House, 2001.

[2] R. J. Adams, Y. Xu, X. Xu, J.-s. Choi, S. D. Gedney, and F. X. Canning, “Modular Fast Direct Electromagnetic Analysis Using Local-Global Solution Modes,” *IEEE Transactions on Antennas and Propagation*, vol. 56, no. 8, pp. 2427–2441, Aug. 2008.

[3] H. Guo, Y. Liu, J. Hu, and E. Michielssen, “A Butterfly-Based Direct Integral-Equation Solver Using Hierarchical LU Factorization for Analyzing Scattering From Electrically Large Conducting Objects,” *IEEE Transactions on Antennas and Propagation*, vol. 65, no. 9, pp. 4742–4750, Sep. 2017.

[4] A. Sharshевsky, Y. Brick, and A. Boag, “Direct Solution of Scattering Problems Using Generalized Source Integral Equations,” *IEEE Transactions on Antennas and Propagation*, vol. 68, no. 7, pp. 5512–5523, Jul. 2020.

<figure><img src="image_5.png" /><figcaption>Fig. 4: Spectrum of the $C-C_c$ for increasing frequencies.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5: Accuracy of the solution and rank of the compressed matrix as a function of $k$, with a compression parameter $\epsilon = 0.015$.</figcaption></figure>

[5] F. P. Andriulli, I. Bogaert, and K. Cools, “On the high frequency behavior and stabilization of a preconditioned and resonance-free formulation,” in *2015 International Conference on Electromagnetics in Advanced Applications (ICEAA)*, Sep. 2015, pp. 1321–1324.

[6] M. Darbas, “Generalized combined field integral equations for the iterative solution of the three-dimensional Maxwell equations,” *Applied Mathematics Letters*, vol. 19, no. 8, pp. 834–839, Aug. 2006.

[7] Y. Boubendir and C. Turc, “Well-conditioned boundary integral equation formulations for the solution of high-frequency electromagnetic scattering problems,” *Computers & Mathematics with Applications*, vol. 67, no. 10, pp. 1772–1805, Jun. 2014.

[8] H. V. Henderson and S. R. Searle, “On Deriving the Inverse of a Sum of Matrices,” *SIAM Review*, vol. 23, no. 1, pp. 53–60, Jan. 1981.