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

Clément Henry2, Davide Consoli1, Alexandre Dély1, Lyes Rahmouni1, Adrien
Merlini2, Francesco P. Andriulli1  1Department of Electronics and
Telecommunications, Politecnico di Torino, Turin, Italy,
<EMAIL_ADDRESS>2Microwave Department, IMT Atlantique, Brest,
France<EMAIL_ADDRESS>
###### 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_{\mathrm{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}^{\mathrm{inc}}$ and the out-of-plane component of the incident magnetic
field $h_{z}^{\mathrm{inc}}$, respectively, are

$\displaystyle(\mathrm{i}k)^{-1}\mathcal{N}j_{t}$
$\displaystyle=-\eta^{-1}e_{t}^{\mathrm{inc}}$ (1)
$\displaystyle\left(\mathcal{I}/2-\mathcal{D}\right)j_{t}$
$\displaystyle=-h_{z}^{\mathrm{inc}}$ (2)

with

$\displaystyle\left(\mathcal{N}j_{t}\right)(\bm{r})$
$\displaystyle\coloneqq-\frac{\partial}{\partial
n}\int_{\Gamma_{2}}\frac{\partial}{\partial
n^{\prime}}g_{\mathrm{2}}(\bm{r},\bm{r}^{\prime})j_{t}(\bm{r}^{\prime})\,\mathrm{d}\bm{r}^{\prime}$
(3) $\displaystyle\left(\mathcal{D}j_{t}\right)(\bm{r})$
$\displaystyle\coloneqq\int_{\Gamma_{2}}\frac{\partial}{\partial
n^{\prime}}g_{\mathrm{2}}(\bm{r},\bm{r}^{\prime})j_{t}(\bm{r}^{\prime})\,\mathrm{d}\bm{r}^{\prime}\,,$
(4)

and
$g_{\mathrm{2}}(\bm{r},\bm{r}^{\prime})\coloneqq\mathrm{i}/4H_{0}^{1}(k|\bm{r}-\bm{r}^{\prime}|)$.
Leveraging the Calderón formula
$\mathcal{S}\mathcal{N}=\mathcal{I}/4+\mathcal{C}_{2}$, where $\mathcal{I}$ is
the identity, $\mathcal{C}_{2}$ is a compact operator, and
$\left(\mathcal{S}j_{t}\right)(\bm{r})\coloneqq\int_{\Gamma_{2}}g_{\mathrm{2}}(\bm{r},\bm{r}^{\prime})j_{t}(\bm{r}^{\prime})\,\mathrm{d}\bm{r}^{\prime}$,
(1) can be turned into a second kind integral equation.

After expanding the unknown current density with piecewise 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}[\mathsfbfit{j}_{2}]_{i}\varphi_{i}$ and after
Galerkin testing, the discretized forms of the normalized Calderón TE-EFIE and
normalized TE-MFIE are

$\displaystyle\mathbf{Z}_{2}\mathsfbfit{j}_{2}=\mathsfbfit{v}_{e,2}\,,$ (5)
$\displaystyle(\mathbf{I}/2-\mathbf{G}_{\varphi\varphi}^{-\frac{1}{2}}\mathbf{D}\mathbf{G}_{\varphi\varphi}^{-\frac{1}{2}})\mathsfbfit{j}_{2}=\mathsfbfit{v}_{h,2}\,,$
(6)

where
$\mathbf{Z}_{2}\coloneqq(\mathrm{i}k)^{-1}\mathbf{G}_{\varphi\varphi}^{-\frac{1}{2}}\mathbf{S}\mathbf{G}_{\varphi\varphi}^{-1}\mathbf{N}\mathbf{G}_{\varphi\varphi}^{-\frac{1}{2}}$,
$\mathsfbfit{v}_{e,2}\coloneqq-\eta^{-1}\mathbf{G}_{\varphi\varphi}^{-\frac{1}{2}}\mathbf{S}\mathbf{G}_{\varphi\varphi}^{-1}\mathsfbfit{e}_{2}$,
$\mathsfbfit{v}_{h,2}\coloneqq-\mathbf{G}_{\varphi\varphi}^{-\frac{1}{2}}\mathsfbfit{h}_{2}$,
$[\mathbf{N}]_{ij}\coloneqq\langle\varphi_{i},\mathcal{N}(\varphi_{j})\rangle$,
$[\mathbf{S}]_{ij}\coloneqq\langle\varphi_{i},\mathcal{S}(\varphi_{j})\rangle$,
$[\mathbf{D}]_{ij}\coloneqq\langle\varphi_{i},\mathcal{D}(\varphi_{j})\rangle$,
$[\mathsfbfit{e}_{2}]_{i}\coloneqq\langle\varphi_{i},e_{t}^{\mathrm{inc}}\rangle$,
$[\mathsfbfit{h}_{2}]_{i}\coloneqq\langle\varphi_{i},h_{z}^{\mathrm{inc}}\rangle$,
$[\mathbf{G}_{\varphi\varphi}]_{ij}\coloneqq\langle\varphi_{i},\varphi_{j}\rangle$,
and $\langle a,b\rangle\coloneqq\int_{\Gamma_{2}}ab\,\mathrm{d}s$.

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{\bm{n}}$. The
electric current density $\bm{j}$ induced on the object by impinging electric
and magnetic fields $\bm{e}^{\mathrm{inc}}$ and $\bm{h}^{\mathrm{inc}}$ is
obtained by solving the 3D EFIE or MFIE

$\displaystyle\bm{\mathcal{T}}\bm{j}\coloneqq\mathrm{i}k\bm{\mathcal{T}}_{\mathrm{A}}\bm{j}+(\mathrm{i}k)^{-1}\bm{\mathcal{T}}_{\mathrm{\Phi}}\bm{j}=-\hat{\bm{n}}\times\bm{e}^{\mathrm{inc}}\,,$
(7)
$\displaystyle\bm{\mathcal{M}}\bm{j}\coloneqq\left(\bm{\mathcal{I}}/2+\bm{\mathcal{K}}\right)\bm{j}=\hat{\bm{n}}\times\bm{h}^{\mathrm{inc}}\,,$
(8)

where
$\left(\bm{\mathcal{T}}_{\mathrm{A}}\bm{j}\right)(\bm{r}):=\hat{\bm{n}}(\bm{r})\times\int_{\Gamma_{3}}g_{3}(\bm{r},\bm{r}^{\prime})\bm{j}(\bm{r}^{\prime})\,\mathrm{d}\bm{r}^{\prime}$,
$\left(\bm{\mathcal{T}}_{\mathrm{\Phi}}\bm{j}\right)(\bm{r}):=-\hat{\bm{n}}(\bm{r})\times\nabla_{\bm{r}}\int_{\Gamma_{3}}g_{3}(\bm{r},\bm{r}^{\prime})\nabla_{\bm{r}^{\prime}}\cdot\bm{j}(\bm{r}^{\prime})\,\mathrm{d}\bm{r}^{\prime}$,
$\left(\bm{\mathcal{K}}\bm{j}\right)(\bm{r}):=-\hat{\bm{n}}(\bm{r})\times\int_{\Gamma_{3}}\nabla_{r}g_{3}(\bm{r},\bm{r}^{\prime})\times\bm{j}(\bm{r}^{\prime})\,\mathrm{d}\bm{r}^{\prime}$,
and
$g_{3}(\bm{r},\bm{r}^{\prime})\coloneqq\exp{(\mathrm{i}k|\bm{r}-\bm{r}^{\prime}|)}/|\bm{r}-\bm{r}^{\prime}|$
is the free-space Green’s function. Leveraging the identity
$\bm{\mathcal{T}}^{2}=-\bm{\mathcal{I}}/4+\bm{\mathcal{C}}_{3}$, a second kind
integral equation can be obtained from the EFIE, yielding the standard
Calderón EFIE [2]

$-\bm{\mathcal{T}}^{2}\bm{j}=\bm{\mathcal{T}}\left(\hat{\bm{n}}\times\bm{e}^{\mathrm{inc}}\right)\,.$
(9)

To discretize (7) following a Petrov-Galerkin approach, we first expand the
unknown current density as
$\bm{j}\approx\sum_{i=1}^{N}[\mathsfbfit{j}_{3}]_{i}\bm{f}_{i}$, where
$\\{\bm{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{\bm{n}}\times\bm{f}_{i}\\}$, Buffa-Christiansen (BC) functions
$\\{\bm{\tilde{f}}_{i}\\}$ [2], rotated BC functions
$\\{\hat{\bm{n}}\times\bm{\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}_{\bm{u}\bm{v}}]_{ij}\coloneqq\langle\bm{u}_{i},\bm{v}_{j}\rangle$,
and
$[\mathbf{\tilde{G}}_{\bm{u}\bm{v}}]_{ij}\coloneqq\langle\hat{\bm{n}}\times\bm{u}_{i},\bm{v}_{j}\rangle$.
Finally, upon adequate Petrov-Galerkin testing, the discrete normalized
Calderón EFIE (9) and normalized MFIE (8) are

$\displaystyle\mathbf{Z}_{3}\mathsfbfit{j}_{3}=\mathsfbfit{v}_{e,3}\,,$ (10)
$\displaystyle\left(\mathbf{I}/2+\mathbf{\tilde{K}}\right)\mathsfbfit{j}_{3}=\mathsfbfit{v}_{h,3}\,,$
(11)

where

$\displaystyle\mathbf{Z}_{3}=-\mathbf{G}_{ff}^{\frac{1}{2}}\mathbf{\tilde{G}}_{\tilde{f}f}^{-1}\mathbf{T_{BC}}\mathbf{\tilde{G}}_{f\tilde{f}}^{-1}\mathbf{T}\mathbf{G}_{ff}^{-\frac{1}{2}}\,,$
(12)
$\displaystyle\mathbf{\tilde{K}}=\mathbf{G}_{ff}^{\frac{1}{2}}\mathbf{\tilde{G}}_{\tilde{f}f}^{-1}\mathbf{K}\mathbf{G}_{ff}^{-\frac{1}{2}}\,,$
(13)
$\displaystyle\mathsfbfit{v}_{e,3}=\mathbf{G}_{ff}^{\frac{1}{2}}\mathbf{\tilde{G}}_{\tilde{f}f}^{-1}\mathbf{T_{BC}}\mathbf{\tilde{G}}_{f\tilde{f}}^{-1}\mathsfbfit{e}_{3}\,,$
(14)
$\displaystyle\mathsfbfit{v}_{h,3}=\mathbf{G}_{ff}^{\frac{1}{2}}\mathbf{\tilde{G}}_{\tilde{f}f}^{-1}\mathsfbfit{h}_{3}\,,$
(15)

$[\mathbf{T}]_{ij}\coloneqq\langle\hat{\bm{n}}\times\bm{f}_{i},\bm{\mathcal{T}}(\bm{f}_{j})\rangle$,
$[\mathbf{K}]_{ij}\coloneqq\langle\hat{\bm{n}}\times\bm{\tilde{f}}_{i},\bm{\mathcal{K}}(\bm{f}_{j})\rangle$,
$[\mathbf{T_{BC}}]_{ij}\coloneqq\langle\hat{\bm{n}}\times\bm{\tilde{f}}_{i},\bm{\mathcal{T}}(\bm{\tilde{f}}_{j})\rangle$,
$[\mathsfbfit{e}_{3}]_{i}\coloneqq\langle\hat{\bm{n}}\times\bm{f}_{i},\hat{\bm{n}}\times\bm{e}^{\mathrm{inc}}\rangle$,
$[\mathsfbfit{h}_{3}]_{i}\coloneqq\langle\hat{\bm{n}}\times\bm{\tilde{f}}_{i},\hat{\bm{n}}\times\bm{h}^{\mathrm{inc}}\rangle$,
and
$\langle\bm{a},\bm{b}\rangle\coloneqq\int_{\Gamma_{3}}\bm{a}\cdot\bm{b}\,\mathrm{d}s$.

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

$[\mathbf{\Sigma}]_{ij}\coloneqq\begin{dcases}\phantom{-}1&\text{if triangle
$j$ is $c_{i}^{+}$,}\\\ -1&\text{if triangle $j$ is $c_{i}^{-}$,}\\\
\phantom{-}0&\text{otherwise,}\end{dcases}$ (16)

and

$[\mathbf{\Lambda}]_{ij}\coloneqq\begin{dcases}\phantom{-}1&\text{if vertex
$j$ is $v_{i}^{+}$,}\\\ -1&\text{if vertex $j$ is $v_{i}^{-}$,}\\\
\phantom{-}0&\text{otherwise.}\end{dcases}$ (17)

In addition to these standard definitions, we also introduce their
counterparts acting on orthonormalized bases
$\mathbf{\tilde{\Sigma}}=\mathbf{G}_{ff}^{-1/2}\mathbf{\Sigma}\mathbf{G}_{pp}^{1/2}$
and
$\mathbf{\tilde{\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]

$\displaystyle\mathbf{P}^{\Sigma}$
$\displaystyle\coloneqq\mathbf{\Sigma}\left(\mathbf{\Sigma}^{T}\mathbf{\Sigma}\right)^{+}\mathbf{\Sigma}^{T}\,,$
(18) $\displaystyle\mathbb{P}^{\Lambda}$
$\displaystyle\coloneqq\mathbf{\Lambda}\left(\mathbf{\Lambda}^{T}\mathbf{\Lambda}\right)^{+}\mathbf{\Lambda}^{T}\,,$
(19)

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

$\displaystyle\mathbf{P}^{\tilde{\Sigma}}$
$\displaystyle\coloneqq\mathbf{\tilde{\Sigma}}\left(\mathbf{\tilde{\Sigma}}^{T}\mathbf{\tilde{\Sigma}}\right)^{+}\mathbf{\tilde{\Sigma}}^{T}\,,$
(20) $\displaystyle\mathbb{P}^{\tilde{\Lambda}}$
$\displaystyle\coloneqq\mathbf{\tilde{\Lambda}}\left(\mathbf{\tilde{\Lambda}}^{T}\mathbf{\tilde{\Lambda}}\right)^{+}\mathbf{\tilde{\Lambda}}^{T}\,.$
(21)

$\bm{e}_{i}$$\bm{v}_{i}^{-}$$\bm{v}_{i}^{+}$$\bm{r}_{i}^{-}$$\bm{r}_{i}^{+}$$c_{i}^{-}$$c_{i}^{+}$
Figure 1: Illustration of the conventions used in defining the RWG functions
$\\{\bm{f}_{i}\\}$: the support of $\bm{f}_{i}$ is the union of the triangles
$c_{i}^{+}$ and $c_{i}^{-}$, that have the RWG-defining, oriented edge
$\bm{e}_{i}$ in common and are complemented by $\bm{r}_{i}^{+}$ and
$\bm{r}_{i}^{-}$, respectively.

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}\coloneqq\mathbf{Z}_{2,3}-\mathbf{I}/4\,,$ (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}}^{T}$
(23)

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

$\left[\mathbf{S}_{\mathbf{X},n}\right]_{ii}=\begin{dcases}\sigma_{\mathbf{X},i}&\text{if
}i>N_{x}-n\,,\\\ \phantom{-}0&\text{otherwise,}\end{dcases}$ (24)

and the filtered matrix $\mathbf{X}_{n}$ as

$\mathbf{X}_{n}\coloneqq\mathbf{V}_{\mathbf{X}}\mathbf{S}_{\mathbf{X},n}\mathbf{V}_{\mathbf{X}}^{T}\,.$
(25)

###  III-A Two-Dimensional Case

Consider the discretized two dimensional variational Laplacian $\mathbf{L}$
where

$[\mathbf{L}]_{ij}=\langle\nabla\varphi_{i},\nabla\varphi_{j}\rangle\,,$ (26)

and its orthonormalized counterpart
$\mathbf{\tilde{L}}\coloneqq\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}\coloneqq\mathbf{V}_{\mathbf{L}}\mathbf{S}_{\mathbf{L},n}\mathbf{V}_{\mathbf{L}}^{T}$
(27)

from which we can define the associated Laplacian filter

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

The reader should note that the SVD is used in the definitions and will never
be necessary in practice, as the FFT will be used to keep the quasi-linear
complexity. In the definitions above, the Laplacian matrix $\mathbf{L}$ should
be replaced by $\mathbf{\tilde{L}}$ when dealing with non-homogeneously
discretized structures as is assumed in the following.

###  III-B Three-Dimensional Case

Following [3], we define the primal filters as

$\displaystyle\mathbf{P}_{n}^{\Sigma}$
$\displaystyle=\mathbf{\Sigma}\left(\mathbf{\Sigma}^{T}\mathbf{\Sigma}\right)^{+}_{n}\mathbf{\Sigma}^{T}\,,$
(29) $\displaystyle\mathbf{P}_{n}^{\Lambda H}$
$\displaystyle=\mathbf{\Lambda}\left(\mathbf{\Lambda}^{T}\mathbf{\Lambda}\right)^{+}_{n}\mathbf{\Lambda}^{T}+\mathbf{I}-\mathbf{P}^{\Sigma}-\mathbb{P}^{\Lambda}\,,$
(30)

and the dual ones as

$\displaystyle\mathbb{P}_{n}^{\Lambda}$
$\displaystyle=\mathbf{\Lambda}\left(\mathbf{\Lambda}^{T}\mathbf{\Lambda}\right)^{+}_{n}\mathbf{\Lambda}^{T}\,,$
(31) $\displaystyle\mathbb{P}_{n}^{\Sigma H}$
$\displaystyle=\mathbf{\Sigma}\left(\mathbf{\Sigma}^{T}\mathbf{\Sigma}\right)^{+}_{n}\mathbf{\Sigma}^{T}+\mathbf{I}-\mathbb{P}^{\Lambda}-\mathbf{P}^{\Sigma}\,.$
(32)

The above loop and star matrices $\mathbf{\Lambda}$ and $\mathbf{\Sigma}$
should be replaced 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.

###  IV-A Two-dimensional case

Considering first the TE-EFIE formulation

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

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

$\left(\mathbf{I}/4+\mathbf{P}^{L}_{n}\mathbf{C}_{2}\right)\mathsfbfit{j}_{2}=\mathsfbfit{v}_{e,2}\,.$
(34)

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

$\left(\mathbf{I}/2+\mathbf{P}^{L}_{n}\mathbf{D}\right)\mathsfbfit{j}_{2}=\mathsfbfit{v}_{h,2}\,.$
(35)

By combining the two previous formulations, after localizing 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}^{L}_{n}\left(\mathbf{C}_{2}+\alpha\mathbf{D}\right)\right)\mathsfbfit{j}_{2}=\mathsfbfit{v}_{e,2}+\alpha\mathsfbfit{v}_{h,2}\,.$
(36)

###  IV-B Three-dimensional case

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

$\left(\mathbf{I}/4+(\mathbf{P}_{n}^{\Lambda
H}+\mathbf{P}_{n}^{\Sigma})\mathbf{C}_{3}\right)\mathsfbfit{j}_{3}=\mathbf{v}_{e,3}\,,$
(37)

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

$\left(\mathbf{I}/2+(\mathbf{P}_{n}^{\Lambda
H}+\mathbf{P}_{n}^{\Sigma})\mathbf{\tilde{K}}\right)\mathsfbfit{j}_{3}=\mathsfbfit{v}_{h,3}\,.$
(38)

By combining the previous two formulations, after localizing the
preconditioner as is standard in Calderón preconditioning, 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}^{\Lambda
H}+\mathbf{P}_{n}^{\Sigma})\left(\mathbf{C}_{3}+\alpha\mathbf{\tilde{K}}\right)\right)\mathsfbfit{j}_{3}\\\
=\mathsfbfit{v}_{e,3}+\alpha\mathsfbfit{v}_{h,3}\,.$ (39)

###  IV-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}\eqqcolon\mathbf{U}\mathbf{V}^{T}$ by using a
skeleton-detecting algorithm [5, 6]. The peculiar single-skeleton structure we
obtain in this way allows for a linear-in-complexity direct inverse in the
form

$\left(\beta\mathbf{I}+\mathbf{C}\right)^{-1}=\left(\beta\mathbf{I}+\mathbf{U}\mathbf{V}^{T}\right)^{-1}=\\\
\mathbf{I}\beta^{-1}-\beta^{-2}\mathbf{U}\left(\mathbf{I}_{r}+\beta^{-1}\mathbf{V}^{T}\mathbf{U}\right)^{-1}\mathbf{V}^{T}\,,$
(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 I: Compression efficiency and solution accuracy obtained with the perturbed circle with a compression tolerance of $\epsilon=${10}^{-3}$$ and filtering point of $200$. Number of unknowns | 1004 | 2008 | 4016 | 8032
---|---|---|---|---
Solution relative error | $0.23\text{\times}{10}^{-3}$ | $0.24\text{\times}{10}^{-3}$ | $0.24\text{\times}{10}^{-3}$ | $0.24\text{\times}{10}^{-3}$
Memory usage without compression (MBytes) | $16$ | $64$ | $258$ | $1032$
Memory usage skeleton (MBytes) | $2.37$ | $1.54$ | $3.08$ | $6.16$

##  V Numerical Results

We first evaluated the fast direct solution scheme on an ellipse characterized
by a semi-major axis of $1.42\text{\,}\mathrm{m}$ and a semi-minor axis of
$1.32\text{\,}\mathrm{m}$ and illuminated by a dipole source radiating at
$k=$0.4\text{\,}\mathrm{rad}\text{\,}{\mathrm{m}}^{-1}$$. Figure 2 shows the
projections of the compact part ($\mathbf{C}$), the filtered compact part
($\mathbf{\tilde{C}}=\mathbf{P}^{L}_{n}\mathbf{C}$), the compressed and
filtered compact part ($\mathbf{U}\mathbf{V}^{T}$), and the dipole right hand
side (RHS) $\mathsfbfit{v}$ on the singular vectors of the filter
$\mathbf{P}^{L}_{N}$. The projection of $\mathbf{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
$\mathbf{\tilde{C}}$ is free from this shortcoming, thanks to the filtering
which allows for the compression of $\mathbf{\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.

Figure 2: Projections of the compact part ($\mathbf{C}$), the filtered compact
part ($\mathbf{\tilde{C}}$), and the dipole RHS on the singular vectors of the
filter $\mathbf{P}^{L}_{N}$. The yellow markers indicate which singular
values/singular vectors are kept after compression. Figure 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\text{\times}{10}^{-6}$$
and a filtering point of $21$.

## 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 Preconditioning,” _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.
  * [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. 20 167–20 172, 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.