# Fast Direct Solvers for Integral Equations at Low-Frequency Based on Operator Filtering

Clément Henry, Davide Consoli, Alexandre Dély, Lyes Rahmouni, Adrien

Merlini, Francesco Andriulli

► To cite this version:

Clément Henry, Davide Consoli, Alexandre Dély, Lyes Rahmouni, Adrien Merlini, et al. .. Fast Di-rect Solvers for Integral Equations at Low-Frequency Based on Operator Filtering. 17th European Conference on Antennas and Propagation (EuCAP 2023), Mar 2023, Florence, Italy. 10.23919/Eu-CAP57121.2023.10133635 . hal-04164122

HAL Id: hal-04164122

https://hal.science/hal-04164122

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 Attribution 4.0 International License
# Fast Direct Solvers for Integral Equations at Low-Frequency Based on Operator Filtering

Clément Henry†, Davide Consoli*, Alexandre Dély*, Lyes Rahmouni*, Adrien Merlini†, Francesco P. Andriulli*

*Department of Electronics and Telecommunications, Politecnico di Torino, Turin, Italy, francesco.andriulli@polito.it

†Microwave Department, IMT Atlantique, Brest, France, clement.henry@imt-atlantique.fr

**Abstract**—This paper focuses on fast direct solvers for integral equations in the low-to-moderate-frequency regime obtained by leveraging preconditioned first kind or second kind operators regularized with Laplacian filters. The spectral errors arising from boundary element discretizations are properly handled by filtering that, in addition, allows for the use of low-rank representations for the compact perturbations of all operators involved. Numerical results show the effectiveness of the approaches and their effectiveness in the direct solution of integral equations.

**Index Terms**—Operator filtering, Calderón preconditioning, fast direct solvers, integral equations.

## I. INTRODUCTION

Several inverse strategies require the solutions of multiple forward problems that, with subsequent optimizations, converge to the desired result. In this context, direct solutions (as opposed to iterative ones) are preferable because they are favorably performing when dealing with multiple right-hand-side problems.

Among forward solution strategies, the boundary element method (BEM) is a powerful technique to model time-harmonic electromagnetic scattering from arbitrary geometries. Its advantages are, among others, that only the scatterer's boundary needs to be discretized and that it automatically enforces the radiation conditions. A popular choice as a forward solution equation is the Electric Field Integral Equation (EFIE). Although this formulation is widespread, its numerical solution comes with its own set of challenges: the linear systems that stem from its discretization are dense and ill-conditioned, with a condition number that grows with the inverses of the frequency and the average mesh edge length [1].

On the one hand, fast algorithms such as the fast multipole method are commonly used to obtain a solution in linear or quasi-linear time and memory complexities. On the other hand, the conditioning problem is usually tackled by Calderón preconditioning [2], multilevel schemes, algebraic approaches, etc. In fact, the classical approach is to separate the problem of acceleration from the one of regularization in order to fix them independently.

This work will present a different paradigm, proposing a new family of techniques that combine these two aspects by taking advantage of the common structure present in fast solvers and preconditioners for the matrices of discretized boundary integral operators. In particular, by an operation of regularization we obtain a formulation from which a fast solver

naturally emerges. This effect, however, cannot be obtained in a straightforward manner because of the spectral pollution of standard integral operators. A key enabling factor will thus be a recently introduced tool: the quasi-Helmholtz Laplacian Filter that will be able to fix the spectral behavior of the operators involved allowing for their fast direct inversion with a focus on problems from low to moderate frequencies.

Numerical results will corroborate the theory, showing the practical relevance of all newly proposed techniques.

## II. NOTATION AND BACKGROUND

Consider a smooth, closed curve $\Gamma_2 \subset \mathbb{R}^2$ modeling a PEC scatterer in a medium of permeability $\mu$, permittivity $\epsilon$, wavenumber $k = \omega\sqrt{\mu\epsilon}$, and impedance $\eta = \sqrt{\mu/\epsilon}$. The TE electric and magnetic field integral equations (TE-EFIE and TE-MFIE) that relate the tangential component of the current $j_t$ with the tangential component of the incident electric field $e_t^{\text{inc}}$ and the out-of-plane component of the incident magnetic field $h_z^{\text{inc}}$, respectively, are

$$(ik)^{-1} \mathcal{N} j_t = -\eta^{-1} e_t^{\text{inc}} \quad (1)$$

$$(\mathcal{I}/2 - \mathcal{D}) j_t = -h_z^{\text{inc}} \quad (2)$$

with

$$(\mathcal{N} j_t)(\mathbf{r}) := -\frac{\partial}{\partial n} \int_{\Gamma_2} \frac{\partial}{\partial n'} g_2(\mathbf{r}, \mathbf{r}') j_t(\mathbf{r}') \, d\mathbf{r}' \quad (3)$$

$$(\mathcal{D} j_t)(\mathbf{r}) := \int_{\Gamma_2} \frac{\partial}{\partial n'} g_2(\mathbf{r}, \mathbf{r}') j_t(\mathbf{r}') \, d\mathbf{r}', \quad (4)$$

and $g_2(\mathbf{r}, \mathbf{r}') := i/4H_0^1(k|\mathbf{r} - \mathbf{r}'|)$. Leveraging the Calderón formula $\mathcal{N} = \mathcal{I}/4 + C_2$, where $\mathcal{I}$ is the identity, $C_2$ is a compact operator (i.e. its spectrum accumulates at the origin), and $(S j_t)(\mathbf{r}) := \int_{\Gamma_2} g_2(\mathbf{r}, \mathbf{r}') j_t(\mathbf{r}') \, d\mathbf{r}'$, (1) can be turned into a second kind integral equation.

After expanding the unknown current density with piece-wise linear Lagrange interpolants $\{\varphi_i\}$, defined on a mesh of $\Gamma_2$ made out of segments of average length $h$, as $j_t \approx \sum_{i=1}^N [\dot{j}_2]_i \varphi_i$ and after Galerkin testing, the discretized forms of the normalized Calderón TE-EFIE and normalized TE-MFIE are

$$\mathbf{Z}_2 \dot{j}_2 = \mathbf{v}_{e,2}, \quad (5)$$

$$(\mathcal{I}/2 - \mathbf{G}_{\varphi\varphi}^{-\frac{1}{2}} \mathbf{D} \mathbf{G}_{\varphi\varphi}^{-\frac{1}{2}}) \dot{j}_2 = \mathbf{v}_{h,2}, \quad (6)$$

where $\mathbf{Z}_2 := (ik)^{-1} \mathbf{G}_{\varphi\varphi}^{-\frac{1}{2}} \mathbf{S} \mathbf{G}_{\varphi\varphi}^{-1} \mathbf{N} \mathbf{G}_{\varphi\varphi}^{-\frac{1}{2}}$, $\mathbf{v}_{e,2} := -\eta^{-1} \mathbf{G}_{\varphi\varphi}^{-\frac{1}{2}} \mathbf{S} \mathbf{G}_{\varphi\varphi}^{-1} \mathbf{e}_2$, $\mathbf{v}_{h,2} := -\mathbf{G}_{\varphi\varphi}^{-\frac{1}{2}} \mathbf{h}_2$, $[\mathbf{N}]_{ij} := \langle \varphi_i, \mathcal{N}(\varphi_j) \rangle$,
$[\mathbf{S}]_{ij} := \langle \varphi_i, \mathcal{S}(\varphi_j) \rangle$, $[\mathbf{D}]_{ij} := \langle \varphi_i, \mathcal{D}(\varphi_j) \rangle$, $[\mathbf{e}_2]_i := \langle \varphi_i, e_t^{\text{inc}} \rangle$, $[\mathbf{h}_2]_i := \langle \varphi_i, h_z^{\text{inc}} \rangle$, $[\mathbf{G}_{\varphi\varphi}]_{ij} := \langle \varphi_i, \varphi_j \rangle$, and $\langle a, b \rangle := \int_{\Gamma_2} ab \, ds$.

In the following, we also consider the 3-dimensional counterpart of the problem described above. Consider a closed, simply-connected domain $\Omega \subset \mathbb{R}^3$ delimited by smooth boundary $\Gamma_3$ and associated outgoing unit surface normal vector field $\hat{n}$. The electric current density $j$ induced on the object by impinging electric and magnetic fields $e^{\text{inc}}$ and $h^{\text{inc}}$ is obtained by solving the 3D EFIE or MFIE

$$ \mathcal{T}j := ik\mathcal{T}_A j + (ik)^{-1}\mathcal{T}_{\Phi}j = -\hat{n} \times e^{\text{inc}}, \quad (7) $$

$$ \mathcal{M}j := (\mathcal{I}/2 + \mathcal{K})j = \hat{n} \times h^{\text{inc}}, \quad (8) $$

where $(\mathcal{T}_A j)(\mathbf{r}) := \hat{n}(\mathbf{r}) \times \int_{\Gamma_3} g_3(\mathbf{r}, \mathbf{r}') j(\mathbf{r}') d\mathbf{r}',$ $(\mathcal{T}_{\Phi}j)(\mathbf{r}) := -\hat{n}(\mathbf{r}) \times \nabla_r \int_{\Gamma_3} g_3(\mathbf{r}, \mathbf{r}') \nabla_{\mathbf{r}'} \cdot j(\mathbf{r}') d\mathbf{r}',$ $(\mathcal{K}j)(\mathbf{r}) := -\hat{n}(\mathbf{r}) \times \int_{\Gamma_3} \nabla_r g_3(\mathbf{r}, \mathbf{r}') \times j(\mathbf{r}') d\mathbf{r}',$ and $g_3(\mathbf{r}, \mathbf{r}') := \exp(ik|\mathbf{r} - \mathbf{r}'|)/|\mathbf{r} - \mathbf{r}'|$ is the free-space Green's function. Leveraging the identity $\mathcal{T}^2 = -\mathcal{I}/4 + \mathcal{C}_3$, a second kind integral equation can be obtained from the EFIE, yielding the standard Calderón EFIE [2]

$$ -\mathcal{T}^2 j = \mathcal{T} (\hat{n} \times e^{\text{inc}}). \quad (9) $$

To discretize (7) following a Petrov-Galerkin approach, we first expand the unknown current density as $j \approx \sum_{i=1}^N [\tilde{f}_i]_i \hat{f}_i$, where $\{\tilde{f}_i\}$ is the set of Rao-Wilton-Glisson (RWG) basis functions defined on a triangular discretization of $\Gamma_3$ in which the average edge length of the triangles is denoted $h$. We also define on this discretized geometry the sets of rotated RWG functions $\{\hat{n} \times \tilde{f}_i\}$, Buffa-Christiansen (BC) functions $\{\tilde{f}_i\}$ [2], rotated BC functions $\{\hat{n} \times \tilde{f}_i\}$, pyramid functions $\{\lambda_i\}$, and patch functions $\{p_i\}$. We also define the (potentially mixed) Gram matrices for these bases as $[\mathbf{G}_{uv}]_{ij} := \langle \mathbf{u}_i, \mathbf{v}_j \rangle$, and $[\mathbf{\tilde{G}}_{uv}]_{ij} := \langle \hat{n} \times \mathbf{u}_i, \mathbf{v}_j \rangle$. Finally, upon adequate Petrov-Galerkin testing, the discrete normalized Calderón EFIE (9) and normalized MFIE (8) are

$$ \mathbf{Z}_3 j_3 = \mathbf{v}_{e,3}, \quad (10) $$

$$ (\mathbf{I}/2 + \tilde{\mathbf{K}}) j_3 = \mathbf{v}_{h,3}, \quad (11) $$

where

$$ \mathbf{Z}_3 = -\mathbf{G}_{ff}^{\frac{1}{2}} \tilde{\mathbf{G}}_{\tilde{f}f}^{-1} \mathbf{T}_{BC} \tilde{\mathbf{G}}_{\tilde{f}f}^{-1} \mathbf{T} \mathbf{G}_{ff}^{-\frac{1}{2}}, \quad (12) $$

$$ \tilde{\mathbf{K}} = \mathbf{G}_{ff}^{\frac{1}{2}} \tilde{\mathbf{G}}_{\tilde{f}f}^{-1} \mathbf{K} \mathbf{G}_{ff}^{-\frac{1}{2}}, \quad (13) $$

$$ \mathbf{v}_{e,3} = \mathbf{G}_{ff}^{\frac{1}{2}} \tilde{\mathbf{G}}_{\tilde{f}f}^{-1} \mathbf{T}_{BC} \tilde{\mathbf{G}}_{\tilde{f}f}^{-1} \mathbf{e}_3, \quad (14) $$

$$ \mathbf{v}_{h,3} = \mathbf{G}_{ff}^{\frac{1}{2}} \tilde{\mathbf{G}}_{\tilde{f}f}^{-1} \mathbf{h}_3, \quad (15) $$

$[\mathbf{T}]_{ij} := \langle \hat{n} \times \mathbf{f}_i, \mathcal{T}(\mathbf{f}_j) \rangle$, $[\mathbf{K}]_{ij} := \langle \hat{n} \times \tilde{\mathbf{f}}_i, \mathcal{K}(\mathbf{f}_j) \rangle$, $[\mathbf{T}_{BC}]_{ij} := \langle \hat{n} \times \tilde{\mathbf{f}}_i, \mathcal{T}(\tilde{\mathbf{f}}_j) \rangle$, $[\mathbf{e}_3]_i := \langle \hat{n} \times \mathbf{f}_i, \hat{n} \times e^{\text{inc}} \rangle$, $[\mathbf{h}_3]_i := \langle \hat{n} \times \tilde{\mathbf{f}}_i, \hat{n} \times h^{\text{inc}} \rangle$, and $\langle a, b \rangle := \int_{\Gamma_3} a \cdot b \, ds$.

To allow for complete analyses of the matrices of the 3-dimensional problem, we introduce the loop-to-RWG and star-to-RWG transformation matrices $\mathbf{\Lambda} \in \mathbb{R}^{N \times N_l}$ and $\mathbf{\Sigma} \in \mathbb{R}^{N \times N_s}$ that are mappings between the solenoidal and RWG subspaces, and between the non-solenoidal and RWG subspaces, respectively. Here $N_s$ is the number of triangles in the discretization of $\Gamma_3$ and $N_l$ the number of vertices. These mappings can be defined using the notation illustrated in Fig. 1 as

<figure><img src="image_2.png" /><figcaption>Fig. 1. Illustration of the conventions used in defining the RWG functions $\{\tilde{f}_i\}$: the support of $\tilde{f}_i$ is the union of the triangles $c_i^+$ and $c_i^-$, that have the RWG-defining, oriented edge $e_i$ in common and are complemented by $r_i^+$ and $r_i^-$, respectively.</figcaption></figure>

$$ [\mathbf{\Sigma}]_{ij} := \begin{cases} 1 & \text{if triangle } j \text{ is } c_i^+, \\ -1 & \text{if triangle } j \text{ is } c_i^-, \\ 0 & \text{otherwise,} \end{cases} \quad (16) $$

and

$$ [\mathbf{\Lambda}]_{ij} := \begin{cases} 1 & \text{if vertex } j \text{ is } v_i^+, \\ -1 & \text{if vertex } j \text{ is } v_i^-, \\ 0 & \text{otherwise.} \end{cases} \quad (17) $$

In addition to these standard definitions, we also introduce their counterparts acting on orthonormalized bases $\tilde{\mathbf{\Sigma}} = \mathbf{G}_{ff}^{-1/2} \mathbf{\Sigma} \mathbf{G}_{pp}^{1/2}$ and $\tilde{\mathbf{\Lambda}} = \mathbf{G}_{ff}^{1/2} \mathbf{\Lambda} \mathbf{G}_{\lambda\lambda}^{-1/2}$. These transformation matrices can also be used to form the quasi-Helmhotlz projectors [2]

$$ \mathbf{P}^{\mathbf{\Sigma}} := \mathbf{\Sigma} (\mathbf{\Sigma}^T \mathbf{\Sigma})^+ \mathbf{\Sigma}^T, \quad (18) $$

$$ \mathbf{P}^{\mathbf{A}} := \mathbf{\Lambda} (\mathbf{\Lambda}^T \mathbf{\Lambda})^+ \mathbf{\Lambda}^T, \quad (19) $$

and their counterparts based $\tilde{\mathbf{\Sigma}}$ and $\tilde{\mathbf{\Lambda}}$

$$ \mathbf{P}^{\tilde{\mathbf{\Sigma}}} := \tilde{\mathbf{\Sigma}} (\tilde{\mathbf{\Sigma}}^T \tilde{\mathbf{\Sigma}})^+ \tilde{\mathbf{\Sigma}}^T, \quad (20) $$

$$ \mathbf{P}^{\tilde{\mathbf{A}}} := \tilde{\mathbf{\Lambda}} (\tilde{\mathbf{\Lambda}}^T \tilde{\mathbf{\Lambda}})^+ \tilde{\mathbf{\Lambda}}^T. \quad (21) $$

From the boundary element matrices $\mathbf{Z}_{2,3}$ of the normalized Calderón TE-EFIE (2D) and normalized Calderón EFIE (3D), we define the matrices $\mathbf{C}_{2,3}$

$$ \mathbf{C}_{2,3} := \mathbf{Z}_{2,3} - \mathbf{I}/4, \quad (22) $$

where $\mathbf{I}$ is the identity matrix. In the following, to simplify the notation, we will omit the subscripts 2 and 3 when the developments apply to both 2-dimensional and 3-dimensional problems.

### III. LAPLACIAN FILTERS

In this section we leverage on operator filtering results presented in [3]. The following notation will be used for the
singular value decomposition (SVD) [4] of a symmetric matrix
$\mathbf{X} \in \mathbb{R}^{N_x \times N_x}$

$$
\mathbf{X} = \mathbf{V}_{\mathbf{X}} \mathbf{S}_{\mathbf{X}} \mathbf{V}_{\mathbf{X}}^{\mathrm{T}} \quad (23)
$$

where $\sigma_{\mathbf{X},i} = [\mathbf{S}_{\mathbf{X}}]_{ii}$ are the singular values of $\mathbf{X}$. We define
the diagonal matrix $\mathbf{S}_{\mathbf{X},n}$ as

$$
[\mathbf{S}_{\mathbf{X},n}]_{ii} = \begin{cases} \sigma_{\mathbf{X},i} & \text{if } i > N_x - n, \\ 0 & \text{otherwise,} \end{cases} \quad (24)
$$

and the filtered matrix $\mathbf{X}_n$ as

$$
\mathbf{X}_n := \mathbf{V}_\mathbf{X} \mathbf{S}_{\mathbf{X},n} \mathbf{V}_\mathbf{X}^T . \quad (25)
$$

A. Two-Dimensional Case

Consider the discretized two dimensional variational Lapla-
cian **L** where

$$
[\mathbf{L}]_{ij} = \langle \nabla \varphi_i, \nabla \varphi_j \rangle, \quad (26)
$$

and its orthonormalized counterpart $\tilde{\mathbf{L}} := \mathbf{G}_{\varphi\varphi}^{-1/2} \mathbf{L} \mathbf{G}_{\varphi\varphi}^{-1/2}$.
Consistently with the general definition (25), we obtain the
filtered Laplacian

$$
\mathbf{L}_n := \mathbf{V}_L \mathbf{S}_{L,n} \mathbf{V}_L^T \quad (27)
$$

from which we can define the associated Laplacian filter

$$
\mathbf{P}_n^L := \mathbf{L}_n^+ \mathbf{L}_n = \mathbf{V}_L \mathbf{S}_{L,n}^+ \mathbf{S}_{L,n} \mathbf{V}_L^T . \quad (28)
$$

The reader should note that the SVD is used in the definitions
and will never be necessary in practice, as the FFT and related
algorithms will be used to keep the quasi-linear complexity
for the sparse matrices. In the definitions above, the Laplacian
matrix **L** should be replaced by **L̃** when dealing with non-
homogeneously discretized structures as is assumed in the
following.

B. Three-Dimensional Case

Following [3], we define the primal filters as

$$
\mathbf{P}_n^\Sigma = \mathbf{\Sigma} (\mathbf{\Sigma}^T \mathbf{\Sigma})_n^+ \mathbf{\Sigma}^T, \quad (29)
$$

$$
\mathbf{P}_n^{AH} = \mathbf{\Lambda} (\mathbf{\Lambda}^T \mathbf{\Lambda})_n^+ \mathbf{\Lambda}^T + \mathbf{I} - \mathbf{P}^\Sigma - \mathbb{P}^A, \quad (30)
$$

and the dual ones as

$$
\mathbb{P}_n^A = \mathbf{\Lambda} (\mathbf{\Lambda}^T \mathbf{\Lambda})_n^+ \mathbf{\Lambda}^T, \quad (31)
$$

$$
\mathbb{P}_n^{\Sigma H} = \mathbf{\Sigma} (\mathbf{\Sigma}^T \mathbf{\Sigma})_n^+ \mathbf{\Sigma}^T + \mathbf{I} - \mathbb{P}^A - \mathbf{P}^\Sigma. \quad (32)
$$

The above loop and star matrices $\mathbf{\Lambda}$ and $\mathbf{\Sigma}$ should be re-
placed by $\tilde{\mathbf{\Lambda}}$ and $\tilde{\mathbf{\Sigma}}$ when dealing with non-homogeneously
discretized structures as is assumed in the following.

IV. FILTERS BASED FAST DIRECT SOLVER

The filters defined in the previous section can be used to
obtain a fast direct solver, in the low-frequency regime, for
the two-dimensional and three-dimensional cases.

A. Two-dimensional case

Considering first the TE-EFIE formulation

$$
\mathbf{Z}_2 \mathbf{j}_2 = \mathbf{v}_{e,2}, \quad (33)
$$

the compact part of which is represented by the matrix $\mathbf{C}_2$,
the filtered equation we propose reads

$$
(\mathbf{I}/4 + \mathbf{P}_n^L \mathbf{C}_2) \mathbf{j}_2 = \mathbf{v}_{e,2}. \quad (34)
$$

Dually, for the TE-MFIE, the compact part is just the operator
$\mathcal{D}$, so that the filtered TM-MFIE reads

$$
(\mathbf{I}/2 + \mathbf{P}_n^L \mathbf{D}) \mathbf{j}_2 = \mathbf{v}_{h,2}. \quad (35)
$$

By combining the two previous formulations, after local-
izing the preconditioner as is standard in Calderón-
preconditioning and grouping the compact filtered parts, we
obtain the filtered Yukawa-Calderón TE-CFIE

$$
\left( \frac{(1 + 2\alpha)}{4} \mathbf{I} + \mathbf{P}_n^L (\mathbf{C}_2 + \alpha \mathbf{D}) \right) \mathbf{j}_2 = \mathbf{v}_{e,2} + \alpha \mathbf{v}_{h,2}. \quad (36)
$$

B. Three-dimensional case

Similarly, in the 3D case, the filtered EFIE we propose is

$$
(\mathbf{I}/4 + (\mathbf{P}_n^{\mathrm{AH}} + \mathbf{P}_n^{\Sigma})\mathbf{C}_3) \mathbf{j}_3 = \mathbf{v}_{e,3}, \quad (37)
$$

and, because its compact part is just the operator $\tilde{\mathbf{K}}$, the filtered
MFIE reads

$$
\left( \mathbf{I} / 2 + (\mathbf{P}_n^{\mathrm{AH}} + \mathbf{P}_n^{\Sigma}) \tilde{\mathbf{K}} \right) \mathbf{j}_3 = \mathbf{v}_{h,3}. \quad (38)
$$

By combining the previous two formulations, after local-
izing the preconditioner as is standard in Calderón precondi-
tioning, and grouping the compact filtered parts, we obtain the
filtered Yukawa-Calderón CFIE

$$
\left( \frac{(1 + 2\alpha)}{4} \mathbf{I} + (\mathbf{P}_n^{\mathrm{AH}} + \mathbf{P}_n^{\Sigma}) (\mathbf{C}_3 + \alpha \tilde{\mathbf{K}}) \right) \mathbf{j}_3 \\
= \mathbf{v}_{e,3} + \alpha \mathbf{v}_{h,3}. \quad (39)
$$

C. Inversion of the structured matrices

Each of the above equations is in the form $\beta \mathbf{I} + \mathbf{C}$, where $\beta$ is a constant and $\mathbf{C}$ is a low rank matrix that can be expressed in the form $\mathbf{C} =: \mathbf{U}\mathbf{V}^T$ by using a skeleton-detecting algorithm such as the adaptive cross approximation (ACA) [5], [6]. The peculiar single-skeleton structure we obtain in this way allows for a linear-in-complexity direct inverse in the form

$$
(\beta \mathbf{I} + \mathbf{C})^{-1} = (\beta \mathbf{I} + \mathbf{U}\mathbf{V}^T)^{-1} = \\
\mathbf{I}\beta^{-1} - \beta^{-2}\mathbf{U}(\mathbf{I}_r + \beta^{-1}\mathbf{V}^T\mathbf{U})^{-1}\mathbf{V}^T, \quad (40)
$$

where $\mathbf{I}_r$ is the identity matrix that has as dimension the rank
of $\mathbf{C}$, using the Woodbury matrix identity [7].
<table><caption>TABLE I
COMPRESSION EFFICIENCY AND SOLUTION ACCURACY OBTAINED WITH THE PERTURBED CIRCLE WITH A COMPRESSION TOLERANCE OF $\epsilon = 10^{-3}$ AND FILTERING POINT OF 200.</caption><thead><tr><th></th><th>1004</th><th>2008</th><th>4016</th><th>8032</th></tr></thead><tbody><tr><td>Number of unknowns</td><td>1004</td><td>2008</td><td>4016</td><td>8032</td></tr><tr><td>Solution relative error</td><td>0.23 × 10<sup>-3</sup></td><td>0.24 × 10<sup>-3</sup></td><td>0.24 × 10<sup>-3</sup></td><td>0.24 × 10<sup>-3</sup></td></tr><tr><td>Memory usage without compression (MBytes)</td><td>16</td><td>64</td><td>258</td><td>1032</td></tr><tr><td>Memory usage skeleton (MBytes)</td><td>2.37</td><td>1.54</td><td>3.08</td><td>6.16</td></tr></tbody></table>

## V. NUMERICAL RESULTS

We first evaluated the fast direct solution scheme on an ellipse characterized by a semi-major axis of 1.42 m and a semi-minor axis of 1.32 m and illuminated by a dipole source radiating at $k = 0.4 \text{ rad m}^{-1}$. Figure 2 shows the projections of the compact part (C), the filtered compact part ($\tilde{C} = P_n^L C$), the compressed and filtered compact part ($UV^T$), and the dipole right hand side (RHS) $v$ on the singular vectors of the filter $P_N^L$. The projection of C on this basis makes apparent the spectral deviation induced by the discretization error (i.e. growth of the projected magnitude for increasing singular value index) which prevents its proper compression using standard techniques. The projection of $\tilde{C}$ is free from this shortcoming, thanks to the filtering which allows for the compression of $\tilde{C}$ at any tolerance above the last singular value filtered out (see Fig. 2). Moreover, the projected RHS in Fig. 2 is band-limited with respect to the filter, which ensures that the filtering does not deteriorate the precision of the solution.

To demonstrate the accuracy of the solution and the compression efficiency, we show in Fig. 3 the error on the solution and the inner rank of the skeleton as a function of the inverse edge length of the discretized ellipse. For increasing refinement, the error with respect to the non compressed operator decreases and the rank saturates to a constant, which validates the effectiveness of the scheme.

Finally, the ellipse has been substituted by a perturbed circle with radius $r(\theta) = 2 + 0.2 \sin(8\theta)$. As shown in Table I, the storage requirement for the skeleton, from a given refinement, grows linearly with the number of unknowns and, at the same time, the error on the solution remains below the desired tolerance.

## ACKNOWLEDGMENT

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".

## REFERENCES

[1] S. B. Adrian, A. Dély, D. Consoli, A. Merlini, and F. P. Andriulli, "Electromagnetic Integral Equations: Insights in Conditioning and Pre-conditioning," *IEEE Open Journal of Antennas and Propagation*, vol. 2, pp. 1143–1174, 2021.

[2] F. P. Andriulli, K. Cools, I. Bogaert, and E. Michielssen, "On a Well-Conditioned Electric Field Integral Operator for Multiply Connected Geometries," *IEEE Transactions on Antennas and Propagation*, vol. 61, no. 4, pp. 2077–2087, Apr. 2013.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Projections of the compact part (C), the filtered compact part ($\tilde{C}$), and the dipole RHS on the singular vectors of the filter $P_N^L$. The yellow markers indicate which singular values/singular vectors are kept after compression.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Accuracy of the solution and inner rank of the skeleton as a function of the inverse edge length, with a compression tolerance of $\epsilon = 6 \times 10^{-6}$ and a filtering point of 21.</figcaption></figure>

[3] A. Merlini, C. Henry, D. Consoli, L. Rahmouni, and F. P. Andriulli, "Laplacian Filters for Integral Equations: Further Developments and Fast Algorithms," in *2022 IEEE Antennas and Propagation Society International Symposium (APS URSI)*, Jul. 2022.

[4] G. H. Golub and C. F. Van Loan, *Matrix Computations*. JHU Press, 2012, vol. 3.

[5] M. Bebendorf, *Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems*, ser. Lecture Notes in Computational Science and Engineering. Berlin: Springer, 2008, no. 63.

[6] E. Liberty, F. Woolfe, P.-G. Martinsson, V. Rokhlin, and M. Tygert, "Randomized algorithms for the low-rank approximation of matrices," *Proceedings of the National Academy of Sciences*, vol. 104, no. 51, pp. 20167–20172, 2007.

[7] 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.