# On the Computation of Square Roots and Inverse Square Roots of Gram Matrices for Surface Integral Equations in Electromagnetics

Rui Chen, Adrien Merlini, Francesco P Andriulli

► To cite this version:

Rui Chen, Adrien Merlini, Francesco P Andriulli. On the Computation of Square Roots and Inverse Square Roots of Gram Matrices for Surface Integral Equations in Electromagnetics. IEEE Transactions on Antennas and Propagation, inPress, 10.1109/tap.2025.3568083. hal-05087045

HAL Id: hal-05087045

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

Submitted on 27 May 2025

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
Communication

On the Computation of Square Roots and Inverse Square Roots of Gram
Matrices for Surface Integral Equations in Electromagnetics

Rui Chen, Adrien Merlini, and Francesco P. Andriulli

**Abstract**—Surface integral equations (SIEs) are widely used for model-
ing electromagnetic scattering problems. However, after their discretiza-
tion via the boundary element method, the spectra and eigenvectors of
resulting matrices are not usually representative of those of underlying
surface integral operators, which can be problematic for methods that
rely heavily on spectral properties. To address this issue, the spectrum of
the integral operators can be recovered, while preserving symmetry, using
square roots and inverse square roots of Gram matrices. In this work,
several algorithms are delineated for the computation of square roots and
inverse square roots of several relevant Gram matrices. The algorithms
we detail rely on properly chosen polynomial expansions of the scalar
square root and inverse square root functions and the theory of matrix
functions. Tables containing different sets of expansion coefficients are
provided along with comparative numerical experiments that evidence
advantages and disadvantages of the different approaches. In addition,
the use of the proposed techniques as a tool for the recovery of the
spectrum of integral operators is illustrated in the case of a spherical
geometry for which the analytic spectrum is known.

***Index Terms***—Gram matrix, matrix function, polynomial expansion, spectral analysis, surface integral equation.

I. INTRODUCTION

The boundary element method (BEM) is widely used for analyz-
ing electromagnetic scattering from perfectly electrically conducting
(PEC) objects. For such scatterers, BEM is typically used to dis-
cretize surface integral equations (SIEs), such as the electric field
integral equation (EFIE), the magnetic field integral equation, and
their linear combinations. Different from differential equation (DE)-
based methods, BEM does not require the use of artificial absorbing
boundary conditions for truncating the physical domain into a finite
computational domain. Moreover, BEM requires only the surface of
the scatterer to be discretized, using boundary elements, which results

This work has received support from the National Natural Science Foun-
dation of China (NSFC) under Grant 62201264 and Grant 62331016, from
the Fundamental Research Funds for the Central Universities under Grant
30924010207, from the Fund Program for the Scientific Activities of Selected
Returned Overseas Professionals in Shanxi Province under Grant 20240063,
from the Foundation of Anhui Key Laboratory of Industrial Energy-Saving
and Safety under Grant KFKT202306, from the European Innovation Council
(EIC) through the European Union's Horizon Europe research Programme
under Grant 101046748 (Project CEREBRO), by the European Union –
Next Generation EU within the PNRR project "Multiscale modeling and
Engineering Applications" of the Italian National Center for HPC, Big Data
and Quantum Computing (Spoke 6) – PNRR M4C2, Investimento 1.4 –
Avviso n. 3138 del 16/12/2021 - CN00000013 National Centre for HPC,
Big Data and Quantum Computing (HPC) - CUP E13C22000990001, and
has received a French government support granted to the Labex CominLabs
excellence laboratory and managed by the National Research Agency in the
"Investing for the Future" program under reference ANR-10-LABX-07-01
(project "CYCLE"). (Corresponding author: Rui Chen.)

Rui Chen was with the Department of Electronics and Telecommunications, Politecnico di Torino, 10129 Turin, Italy. He is now with the School of Microelectronics, Nanjing University of Science and Technology, Nanjing 210094, China and also with North Automatic Control Technology Institute, Taiyuan 030006, China (e-mail: rui.chen@njust.edu.cn).

Adrien Merlini is with the Microwave Department, IMT Atlantique, 29238 Brest, France (e-mail: adrien.merlini@imt-atlantique.fr).

Francesco P. Andriulli is with the Department of Electronics and Telecommunications, Politecnico di Torino, 10129 Turin, Italy (e-mail: francesco.andriulli@polito.it).

in comparatively fewer unknowns compared to the DE-based methods
[1].

To numerically solve SIEs, surface integral operators are dis-
cretized by expanding the surface current with basis functions and
testing the resulting equation, which yields a matrix system. Then,
direct solvers or iterative solvers coupled with fast matrix-vector
product algorithms are used to solve the system for the current
expansion coefficients.

However, the spectrum and eigenvectors of the resulting boundary
element matrices are not usually representative of those of the
underlying surface integral operators. It is partly due to the fact
that, in general, the boundary elements do not form an orthonormal
basis of the spaces they span and the operators' eigenvalues can
be recovered through a generalized eigenvalue problem [2]. While
it is not necessarily problematic for standard solution strategies,
some numerical methods rely on delicate spectral manipulations that
necessitate the correspondence between the spectra of the operator to
be discretized and that of the corresponding boundary element matrix
[3]. To solve this issue, one could use orthonormal basis functions,
but they are often impractical.

Other approaches leveraging the Gram matrices [4] of the bases
used in the discretization of the integral operators have been proposed.
For instance, when the symmetry of the original problem does not
need to be maintained, the multiplication by the inverse of the
Gram matrices [5] can be used to reconstruct the operator spectrum.
Instead, when the symmetry must be preserved, an *a posteriori*
orthonormalization of the bases using square roots and inverse square
roots of these Gram matrices—as delineated for instance in [3]
and references therein—can be used. To compute the square roots
and inverse square roots of the Gram matrices, the singular value
decomposition (SVD) [4] could be used. However, because of its
prohibitive computational cost, SVD quickly becomes impractical for
large problems. Alternatively, several algorithms have been presented
to compute these matrices—such as the Schur method [6], the Newton
method [7], or on the matrix sign function [8] to name a few—and
their theory is well-established. Nonetheless, the application of these
schemes to the variety of the Gram matrices appearing in the integral
equations-relevant scenarios is not trivial. A specific set of approaches
must be selected among many and the associated parameters need to
be optimized to avoid instability issues [7] and ensure a computational
complexity compatible with traditional fast solution and iterative
eigenproblem methods [9].

In this work, several strategies inherited from scalar polynomial
expansions approaches are delineated and tailored for the evaluation
of the square roots and inverse square roots of several of the most
relevant Gram matrices for SIEs in electromagnetics. To numerically
compute these matrices suitably chosen polynomial expansions of
the scalar square root and inverse square root functions are used
in conjunction with the theory of matrix functions [10]. Note that
the techniques delineated here could also be used in applications
beyond electromagnetics in which Gram matrices are useful, for
instance, in quantum mechanics and quantum theory [11]. Tables
containing expansion coefficients are provided along with numerical
experiments that illustrate the merits and the shortcomings of the different approaches along with a thorough performance analysis. To demonstrate the applicability of the proposed techniques, numerical examples showing how the spectrum of the integral operators can be recovered from a standard boundary element matrix are presented for the scattering from a PEC sphere for which the analytic spectrum is known. A preliminary version of this manuscript with limited implementation details was presented in [12].

## II. BACKGROUND

Let $S$ denote the surface of a PEC scatterer residing in an unbounded homogeneous background medium with permittivity $\epsilon$ and permeability $\mu$. The PEC scatterer is excited by an electromagnetic wave with the incident electric field $\mathbf{E}^{\text{inc}}(\mathbf{r})$. Upon excitation, the surface current $\mathbf{J}(\mathbf{r})$ induced on $S$ generates the scattered electric field $\mathbf{E}^{\text{sca}}(\mathbf{r})$ in the background medium. The total electric field $\mathbf{E}^{\text{tot}}(\mathbf{r}) = \mathbf{E}^{\text{inc}}(\mathbf{r}) + \mathbf{E}^{\text{sca}}(\mathbf{r})$ satisfies the boundary condition $\hat{\mathbf{n}}(\mathbf{r}) \times \mathbf{E}^{\text{tot}}(\mathbf{r}) = 0$ for $\mathbf{r} \in S$ from which EFIE is expressed as

$$ -\hat{\mathbf{n}}(\mathbf{r}) \times \mathbf{E}^{\text{inc}}(\mathbf{r}) = \mathcal{T}[\mathbf{J}(\mathbf{r})] = \hat{\mathbf{n}}(\mathbf{r}) \times \left\{ -j\omega\mu \int_S g(\mathbf{r}, \mathbf{r}') \mathbf{J}(\mathbf{r}') ds' - \frac{j}{\omega\epsilon} \nabla \int_S g(\mathbf{r}, \mathbf{r}') \nabla_s' \mathbf{J}(\mathbf{r}') ds' \right\}. \quad (1) $$

Here, $\hat{\mathbf{n}}(\mathbf{r})$ is the outward pointing unit normal vector, $g(\mathbf{r}, \mathbf{r}') = e^{-j\kappa R}/(4\pi R)$ is the Green function, $k = \omega\sqrt{\mu\epsilon}$ is the wavenumber, and $R = |\mathbf{r} - \mathbf{r}'|$ represents the distance between the field and source points. Note that, the time-harmonic factor $e^{j\omega t}$ is assumed and suppressed throughout this work.

To numerically solve (1) for $\mathbf{J}(\mathbf{r})$, $S$ is meshed using triangular patches. Then, $\mathbf{J}(\mathbf{r})$ is expanded using Rao-Wilton-Glisson (RWG) basis functions [13] as $\mathbf{J}(\mathbf{r}) = \sum_{n=1}^{N_E} \{\mathbf{j}_n\}_n f_n(\mathbf{r})$. Here, $N_E$ denotes the number of edges of the mesh, $\mathbf{j}$ is the vector containing the $N_E$ expansion coefficients to be solved for, and the $n^{\text{th}}$ RWG basis function $f_n(\mathbf{r})$ is defined as $f_n(\mathbf{r}) = (\mathbf{r} - \mathbf{r}_n^+)/(2A_n^+)$ if $\mathbf{r} \in t_n^+$ and $f_n(\mathbf{r}) = (\mathbf{r}_n^- - \mathbf{r})/(2A_n^-)$ if $\mathbf{r} \in t_n^-$, where $t_n^+$ and $t_n^-$ are two triangular patches sharing the $n^{\text{th}}$ edge, $\mathbf{r}_n^+$ and $\mathbf{r}_n^-$ are the vertices of $t_n^+$ and $t_n^-$ that do not belong to the $n^{\text{th}}$ edge, and $A_n^+$ and $A_n^-$ are the areas of $t_n^+$ and $t_n^-$, respectively.

Next, rotated RWG functions $\hat{\mathbf{n}}(\mathbf{r}) \times f_m(\mathbf{r})$, $m = 1, 2, ..., N_E$, are used to test the equation yielding a fully discretized matrix equation

$$ \mathbf{T} \mathbf{j} = \mathbf{v}. \quad (2) $$

Here, the entries of the boundary element matrix $\mathbf{T}$ has dimension $N_E \times N_E$ and the right-hand side vector $\mathbf{v}$ has $N_E$ elements. Their elements are $\{\mathbf{T}\}_{mn} = \langle \hat{\mathbf{n}}(\mathbf{r}) \times f_m(\mathbf{r}), \mathcal{T}[\mathbf{f}_n](\mathbf{r}) \rangle$ and $\{\mathbf{v}\}_m = -\langle \hat{\mathbf{n}}(\mathbf{r}) \times f_m(\mathbf{r}), \hat{\mathbf{n}}(\mathbf{r}) \times \mathbf{E}^{\text{inc}}(\mathbf{r}) \rangle$, respectively, where $\langle \cdot, \cdot \rangle$ denotes the appropriate bilinear form and, which for all cases studied in the manuscript, can be computed as $\langle \mathbf{a}(\mathbf{r}), \mathbf{b}(\mathbf{r}) \rangle = \int_{S_a} \mathbf{a}(\mathbf{r}) \cdot \mathbf{b}(\mathbf{r}) ds$ for vector functions, $\langle \mathbf{a}(\mathbf{r}), \mathbf{b}(\mathbf{r}) \rangle = \int_{S_b} \mathbf{a}(\mathbf{r})\mathbf{b}(\mathbf{r}) ds$ for scalar functions, and $S_a$ is the support of $\mathbf{a}(\mathbf{r})$ or $\mathbf{a}(\mathbf{r})$.

For certain applications, such as preconditioning [14], the electric field integral operator (EFIO) $\mathcal{T}$ must also be discretized using Buffa-Christiansen (BC) basis functions $g_n(\mathbf{r})$, $n = 1, 2, ..., N_E$ [15]. The BC basis functions are dual functions with respect to the RWG basis functions and are defined on barycentrically refined triangular patches. The explicit definition of the BC basis function $g_n(\mathbf{r})$ is omitted here and can be found in [15]. We denote the Gram matrices associated with the RWG and BC functions as $\mathbf{G}_{\mathbf{r},f}$ and $\mathbf{G}_{\mathbf{g},g}$. Their entries are $\{\mathbf{G}_{\mathbf{r},f}\}_{mn} = \langle f_m(\mathbf{r}), f_n(\mathbf{r}) \rangle$ and $\{\mathbf{G}_{\mathbf{g},g}\}_{mn} = \langle g_m(\mathbf{r}), g_n(\mathbf{r}) \rangle$, respectively.

Besides the above mentioned vector basis functions, used for the discretization of EFIE, scalar basis functions can also be used for

the discretization of SIEs. Here, two commonly used sets of scalar basis functions are considered the pyramid basis functions $\lambda_n(\mathbf{r})$, $n = 1, 2, ..., N_V$ and their dual functions $\tilde{\lambda}_n(\mathbf{r})$, $n = 1, 2, ..., N_P$, where $N_V$ and $N_P$ denote the number of vertices and triangular patches of the mesh, respectively. The dual pyramid basis functions $\tilde{\lambda}_n(\mathbf{r})$ are defined on the barycentrically refined triangular patches. The $n^{\text{th}}$ pyramid basis function $\lambda_n(\mathbf{r})$ is defined to be equal to 1 at the $n^{\text{th}}$ vertex, 0 at the other vertices, and linear on the triangular patches sharing the $n^{\text{th}}$ vertex [16]. The explicit definition of $\tilde{\lambda}_n(\mathbf{r})$ is omitted here and can be found in [15]. Similarly, the Gram matrices associated with the pyramid and dual pyramid functions are $\mathbf{G}_{\lambda,\lambda}$ with dimension $N_V \times N_V$ and $\mathbf{G}_{\tilde{\lambda},\tilde{\lambda}}$ with dimension $N_P \times N_P$, whose entries are $\{\mathbf{G}_{\lambda,\lambda}\}_{mn} = \langle \lambda_m(\mathbf{r}), \lambda_n(\mathbf{r}) \rangle$ and $\{\mathbf{G}_{\tilde{\lambda},\tilde{\lambda}}\}_{mn} = \langle \tilde{\lambda}_m(\mathbf{r}), \tilde{\lambda}_n(\mathbf{r}) \rangle$, respectively.

The Gram matrices introduced above ($\mathbf{G}_{\mathbf{f},f}$, $\mathbf{G}_{\mathbf{g},g}$, $\mathbf{G}_{\lambda,\lambda}$, and $\mathbf{G}_{\tilde{\lambda},\tilde{\lambda}}$) are symmetric positive definite (SPD)—since we consider only Gram matrices of bases [4]—and can be used for the normalizations of the EFIE boundary element matrix $\mathbf{T}$, the loop-to-RWG transformation matrix $\mathbf{\Lambda}$, and the star-to-RWG transformation matrix $\mathbf{\Sigma}$, where the explicit definitions of $\mathbf{\Lambda}$ and $\mathbf{\Sigma}$ can be found in [17]. As mentioned in [3], using the matrices multiplied by the inverse square roots of the corresponding Gram matrices $\hat{\mathbf{T}}$, $\hat{\mathbf{\Lambda}}$, and $\hat{\mathbf{\Sigma}}$ instead of their standard counterparts $\mathbf{T}$, $\mathbf{\Lambda}$, and $\mathbf{\Sigma}$ allows for more efficient spectral manipulations for problems involving non-uniformly discretized geometries. These new matrices are obtained by multiplying the standard matrices with the inverse square roots of the corresponding Gram matrices, e.g., $\hat{\mathbf{T}} = \sqrt{\mathbf{G}_{\mathbf{f},f}^{-1}}\mathbf{T}\sqrt{\mathbf{G}_{\mathbf{f},f}^{-1}}$. In Section III, strategies for the computation of the square roots and inverse square roots of the real, sparse, and SPD Gram matrices introduced above are presented.

## III. COMPUTATION OF THE SQUARE ROOTS AND INVERSE SQUARE ROOTS OF THE GRAM MATRICES

To compute the square roots and inverse square roots of the Gram matrices of interest, a simple-to-define but expensive-to-compute approach is leveraging SVD [4]. SVD of a real SPD matrix $\mathbf{G}$ of dimension $N \times N$ is

$$ \mathbf{G} = \mathbf{U}\mathbf{S}\mathbf{U}^T, \quad (3) $$

where $\mathbf{S}$ is a diagonal matrix with non-zero elements that are the $N$ singular values of $\mathbf{G}$, $\mathbf{U}$ is a unitary matrix. Here and in the following $\mathbf{G}$ stands for any of the matrices $\mathbf{G}_{\mathbf{f},f}$, $\mathbf{G}_{\mathbf{g},g}$, $\mathbf{G}_{\lambda,\lambda}$, and $\mathbf{G}_{\tilde{\lambda},\tilde{\lambda}}$. Given that the Gram matrices under consideration are also invertible, their square roots and inverse square roots can be obtained as

$$ \sqrt{\mathbf{G}} = \mathbf{U}\bar{\mathbf{S}}\mathbf{U}^T \quad (4) $$

$$ \sqrt{\mathbf{G}}^{-1} = \mathbf{U}\bar{\mathbf{S}}^{-1}\mathbf{U}^T, \quad (5) $$

where $\bar{\mathbf{S}}$ is a diagonal matrix of which the diagonal entries are the square roots of that of $\mathbf{S}$. However, because of the computational cost of the SVD algorithm, this method quickly becomes unpractical for large problems. Note that, since the matrices under consideration are SPD, their respective square roots are unique [18].

In this work, the square roots and inverse square roots of the Gram matrices are expressed as

$$ \sqrt{\mathbf{G}} = \sqrt{\|\mathbf{G}\|_2} \sqrt{\mathbf{G}/\|\mathbf{G}\|_2}, \quad (6) $$

$$ \sqrt{\mathbf{G}}^{-1} = \sqrt{\|\mathbf{G}\|_2^{-1}} \sqrt{\mathbf{G}/\|\mathbf{G}\|_2^{-1}}, \quad (7) $$

where $\|\mathbf{G}\|_2$ denotes the $L_2$-norm of $\mathbf{G}$. Even though this normalization by $\|\mathbf{G}\|_2$ is not a necessary step of our proposed techniques, it is helpful for our analysis and the value of $\|\mathbf{G}\|_2$ can be computed fast using power methods [19] for sparse, SPD, and well-conditioned
matrices. Then, because the Gram matrices under consideration are sparse, $\sqrt{\mathbf{G}}/\|\mathbf{G}\|_2$ and $\sqrt{\mathbf{G}}/\|\mathbf{G}\|_2^{-1}$ are numerically approximated using polynomial expansions-based algorithms combined with the theory of matrix functions [10] as presented in Section III-A, -B, and -C. Finally, $\sqrt{\mathbf{G}}$ and $\sqrt{\mathbf{G}}^{-1}$ are computed by multiplying $\sqrt{\|\mathbf{G}\|_2}$ and $\sqrt{\|\mathbf{G}\|_2^{-1}}$ with the numerical approximations of $\sqrt{\mathbf{G}}/\|\mathbf{G}\|_2$ and $\sqrt{\mathbf{G}}/\|\mathbf{G}\|_2^{-1}$, respectively. This approach is made possible by the fact that all the Gram matrices under consideration are sparse, meaning that they can be computed, stored, and multiplied to a vector in linear complexity, which is crucial to maintain compatibility fast solvers for the overall electromagnetic problem.

To condense the treatments, given a real, sparse, and SPD matrix $\mathbf{X}$ of dimension $N \times N$, we denote $f(\mathbf{X})$ as a matrix function that corresponds to either $\mathbf{X} \mapsto \sqrt{\mathbf{X}}$ or $\mathbf{X} \mapsto \sqrt{\mathbf{X}}^{-1}$. The numerical expression of $f(\mathbf{X})$ can be obtained by generalizing the numerical (polynomial or rational) approximation of its scalar function counterpart $f(x)$, which corresponds to either $x \mapsto \sqrt{x}$ or $x \mapsto \sqrt{x}^{-1}$, $x \in \mathbb{R}^{+}$, to matrix arguments [10]. Then, the numerical approximations of $\sqrt{\mathbf{G}}/\|\mathbf{G}\|_2$ and $\sqrt{\mathbf{G}}/\|\mathbf{G}\|_2^{-1}$ are obtained by using the expansions of $\sqrt{\mathbf{X}}$ and $\sqrt{\mathbf{X}}^{-1}$ with $\mathbf{X} = \mathbf{G}/\|\mathbf{G}\|_2$. We denote by $\delta$ the relative error between $f(\mathbf{X})$—computed, for instance, by SVD—and its approximation $f^{\text{num}}(\mathbf{X})$: $\delta = \|f(\mathbf{X}) - f^{\text{num}}(\mathbf{X})\|_2/\|f(\mathbf{X})\|_2$.

### A. Taylor Series Expansion

As the first algorithmic approach, the Taylor series expansion (TSE) [20] is used for the polynomial approximations of $x \mapsto \sqrt{x}$ and $x \mapsto \sqrt{x}^{-1}$. Using TSEs, $x \mapsto \sqrt{x}$ and $x \mapsto \sqrt{x}^{-1}$ can be numerically expanded around $x = 1$ as

$$x \mapsto \sqrt{x} \approx \sum_{n=0}^{N_T} c_n^{\text{TSE}}(x-1)^n \quad \text{and} \quad x \mapsto \sqrt{x}^{-1} \approx \sum_{n=0}^{N_T} \tilde{c}_n^{\text{TSE}}(x-1)^n. \qquad (8)$$

Here, $N_T$ is the order of the TSE, and $c_n^{\text{TSE}}$ and $\tilde{c}_n^{\text{TSE}}$ are the $n^{\text{th}}$ expansion coefficients of TSE for the approximations of $\sqrt{x}$ and $\sqrt{x}^{-1}$, respectively. These coefficients can be computed as

$$c_n^{\text{TSE}} = \frac{\sqrt{x}^{(n)}}{n!} \quad \text{and} \quad \tilde{c}_n^{\text{TSE}} = \frac{(1/\sqrt{x})^{(n)}}{n!}. \qquad (9)$$

The approximations of their matrix function counterparts $\mathbf{X} \mapsto \sqrt{\mathbf{X}}$ and $\mathbf{X} \mapsto \sqrt{\mathbf{X}}^{-1}$ are obtained by generalizing (8) to $\mathbf{X}$ as

$$\mathbf{X} \mapsto \sqrt{\mathbf{X}} \approx \sum_{n=0}^{N_T} c_n^{\text{TSE}}(\mathbf{X}-\mathbf{I})^n, \quad \mathbf{X} \mapsto \sqrt{\mathbf{X}}^{-1} \approx \sum_{n=0}^{N_T} \tilde{c}_n^{\text{TSE}}(\mathbf{X}-\mathbf{I})^n, \qquad (10)$$

where $\mathbf{I}$ is the identity matrix with the same dimension as $\mathbf{X}$.

The values of the first ten expansion coefficients of the TSEs $\{c_n^{\text{TSE}}\}_{n=0}^9$ and $\{\tilde{c}_n^{\text{TSE}}\}_{n=0}^9$ are $\{1, \frac{1}{2}, -\frac{1}{8}, \frac{1}{16}, -\frac{5}{128}, \frac{7}{256}, -\frac{21}{1024}, \frac{33}{2048}, -\frac{429}{32768}, \frac{715}{65536}\}$ and $\{1, -\frac{1}{8}, \frac{3}{8}, -\frac{5}{16}, \frac{35}{128}, -\frac{63}{256}, \frac{231}{1024}, -\frac{429}{2048}, \frac{6435}{32768}, -\frac{12155}{65536}\}$, respectively.

### B. Chebyshev Polynomial Expansion

As the second algorithmic approach, the truncated Chebyshev polynomial series expansion (CPE) [21] is used for the polynomial approximations of the square root function and its reciprocal on the interval $x \in [n_0, 1]$ as

$$x \mapsto \sqrt{x} \approx \sum_{n=0}^{N_C} c_n^{\text{CPE}} T_n(x) \quad \text{and} \quad x \mapsto \sqrt{x}^{-1} \approx \sum_{n=0}^{N_C} \tilde{c}_n^{\text{CPE}} T_n(x). \qquad (11)$$

Here, the primed summation indicates that the first term is halved, $n_0$ is a positive real number, $N_C$ is the expansion order of CPE,

$T_n(x)$ is the shifted Chebyshev polynomial of the first kind of order $n$ defined on the interval $x \in [n_0, 1]$ as

$$T_n(x) = \begin{cases} 1 & \text{if } n = 0 \\ \frac{2x-(n_0+1)}{1-n_0} & \text{if } n = 1 \\ 2 \left( \frac{2x-(n_0+1)}{1-n_0} \right) T_{n-1}(x) - T_{n-2}(x) & \text{otherwise} \end{cases}, \qquad (12)$$

and the Chebyshev series coefficients $c_n^{\text{CPE}}$ and $\tilde{c}_n^{\text{CPE}}$ for the approximations of $\sqrt{x}$ and $\sqrt{x}^{-1}$ can be computed as [21]

$$c_n^{\text{CPE}} = \frac{2}{\pi} \int_{n_0}^{1} \frac{\sqrt{x} T_n(x)}{\sqrt{-x^2 + (n_0 + 1)x - n_0}} dx \qquad (13)$$

$$\tilde{c}_n^{\text{CPE}} = \frac{2}{\pi} \int_{n_0}^{1} \frac{\sqrt{x}^{-1} T_n(x)}{\sqrt{-x^2 + (n_0 + 1)x - n_0}} dx, \qquad (14)$$

which can be rewritten with special functions as

$$c_n^{\text{CPE}} = 2\sqrt{n_0} \frac{{}_3F_2\left(\begin{matrix} -1/2 & 1/2 & 1 \\ 1-n & 1+n & \end{matrix}; 1-n_0^{-1}\right)}{\Gamma(1-n)\Gamma(1+n)} \qquad (15)$$

$$\tilde{c}_n^{\text{CPE}} = \frac{2}{\sqrt{n_0}} \frac{{}_3F_2\left(\begin{matrix} 1/2 & 1/2 & 1 \\ 1-n & 1+n & \end{matrix}; 1-n_0^{-1}\right)}{\Gamma(1-n)\Gamma(1+n)}, \qquad (16)$$

where ${}_3F_2$ is the generalized hypergeometric function [22] and $\Gamma$ is the Euler gamma function.

Similarly, after the series expansions of $x \mapsto \sqrt{x}$ and $x \mapsto \sqrt{x}^{-1}$ in (11) are derived, their matrix function counterparts can be obtained by generalizing (11) to the matrix arguments as

$$\mathbf{X} \mapsto \sqrt{\mathbf{X}} \approx \sum_{n=0}^{N_C} \tilde{c}_n^{\text{CPE}} T_n(\mathbf{X}) \quad \mathbf{X} \mapsto \sqrt{\mathbf{X}}^{-1} \approx \sum_{n=0}^{N_C} \tilde{c}_n^{\text{CPE}} T_n(\mathbf{X}), \qquad (17)$$

where direct evaluation of the sum in (17) is not necessarily optimal and other algorithms can be employed [23], and $T_n(\mathbf{X})$ is defined as

$$T_n(\mathbf{X}) = \begin{cases} \mathbf{I} & \text{if } n = 0 \\ \frac{2\mathbf{X}-(n_0+1)\mathbf{I}}{1-n_0} & \text{if } n = 1 \\ \frac{2[2\mathbf{X}-(n_0+1)\mathbf{I}]}{1-n_0} T_{n-1}(\mathbf{X}) - T_{n-2}(\mathbf{X}) & \text{otherwise} \end{cases}. \qquad (18)$$

Note that, due to the normalization and the properties of the Gram matrices we study, the spectrum of $\mathbf{X}$—to which the polynomial approximations are applied—lives in the interval $[n_0, 1]$, where $n_0$ is the inverse of the condition number of $\mathbf{X}$. When using (17), the coefficients $c_n^{\text{CPE}}$ and $\tilde{c}_n^{\text{CPE}}$ must be computed for each matrix, since they depend on the condition number of the matrix argument $n_0^{-1}$. If this turns out to be inconvenient, tabulated values of coefficients can be used for the matrices that have a value of $n_0$ living in known ranges. Usage of such tabulated values will typically yield higher errors than the computation with the matrix-dependent coefficients. In Section IV, we will provide numerical results for these two approaches designated as “CPE-1” for the matrix specific approach and “CPE-2” for the tabulated approach, respectively. For convenience, the values of $c_n^{\text{CPE}}$ and $\tilde{c}_n^{\text{CPE}}$ of CPE-2 are provided in Table I and II for given bounds on the matrices conditioning (different $n_0$).

### C. Padé Approximant Expansion

The last algorithmic approach we explore is the Padé approximant expansion (PAE) [24] for $x \mapsto \sqrt{x}$ and $x \mapsto \sqrt{x}^{-1}$ at $x = 1$ as

$$x \mapsto \sqrt{x} \approx \frac{\sum_{n=0}^{N_A} c_n^{\text{PAE}} x^n}{\sum_{n=0}^{N_A} c_n^{\text{PAE}} x^{N_A-n}}, \quad x \mapsto \sqrt{x}^{-1} \approx \frac{\sum_{n=0}^{N_A} c_n^{\text{PAE}} x^{N_A-n}}{\sum_{n=0}^{N_A} c_n^{\text{PAE}} x^n}. \qquad (19)$$
<table><caption>TABLE I
FIRST TWENTY CHEBYSHEV SERIES COEFFICIENTS FOR THE SQUARE ROOT FUNCTION FOR DIFFERENT $n_0$.</caption><thead><tr><th>n<sub>0</sub></th><th>c<sub>n</sub><sup>CPE</sup></th><th>c<sub>0</sub><sup>CPE</sup></th><th>c<sub>1</sub><sup>CPE</sup></th><th>c<sub>2</sub><sup>CPE</sup></th><th>c<sub>3</sub><sup>CPE</sup></th><th>c<sub>4</sub><sup>CPE</sup></th><th>c<sub>5</sub><sup>CPE</sup></th><th>c<sub>6</sub><sup>CPE</sup></th><th>c<sub>7</sub><sup>CPE</sup></th><th>c<sub>8</sub><sup>CPE</sup></th><th>c<sub>9</sub><sup>CPE</sup></th></tr><tr><th></th><th></th><th>c<sub>0</sub><sup>CPE</sup></th><th>c<sub>1</sub><sup>CPE</sup></th><th>c<sub>2</sub><sup>CPE</sup></th><th>c<sub>3</sub><sup>CPE</sup></th><th>c<sub>4</sub><sup>CPE</sup></th><th>c<sub>5</sub><sup>CPE</sup></th><th>c<sub>6</sub><sup>CPE</sup></th><th>c<sub>7</sub><sup>CPE</sup></th><th>c<sub>8</sub><sup>CPE</sup></th><th>c<sub>9</sub><sup>CPE</sup></th></tr><tr><th></th><th></th><th></th><th>c<sub>10</sub><sup>CPE</sup></th><th>c<sub>11</sub><sup>CPE</sup></th><th>c<sub>12</sub><sup>CPE</sup></th><th>c<sub>13</sub><sup>CPE</sup></th><th>c<sub>14</sub><sup>CPE</sup></th><th>c<sub>15</sub><sup>CPE</sup></th><th>c<sub>16</sub><sup>CPE</sup></th><th>c<sub>17</sub><sup>CPE</sup></th><th>c<sub>18</sub><sup>CPE</sup></th><th>c<sub>19</sub><sup>CPE</sup></th></tr></thead><tbody><tr><td>n<sub>0</sub> ≥ 1 × 10<sup>-1</sup></td><td>50720584<br>36057897</td><td>58192354<br>176379781</td><td>—<br>1188</td><td>—<br>409</td><td>599673<br>145</td><td>—<br>101639721</td><td>—<br>25723455</td><td>—<br>15</td><td>58340<br>103329721</td><td>—<br>15</td><td>7566<br>30591204</td><td>—<br>2764<br>79582277</td></tr><tr><td>n<sub>0</sub> ≥ 5 × 10<sup>-2</sup></td><td>141760564<br>104989387</td><td>35002745<br>95238932</td><td>—<br>4645</td><td>—<br>4539</td><td>—<br>154958201</td><td>—<br>131648930</td><td>—<br>567</td><td>—<br>61162961</td><td>—<br>225</td><td>—<br>74184637</td><td>—<br>91</td><td>—<br>89427045</td><td>—<br>43265<br>239826181</td></tr><tr><td>n<sub>0</sub> ≥ 1 × 10<sup>-2</sup></td><td>94293998<br>72892523</td><td>11654085<br>28548136</td><td>—<br>43604</td><td>—<br>29207</td><td>—<br>66081562</td><td>—<br>162168061</td><td>—<br>110193549</td><td>—<br>16229961</td><td>—<br>139399337</td><td>—<br>106738132</td><td>—<br>18</td><td>—<br>129639190</td><td>—<br>77<br>52900755</td></tr><tr><td>n<sub>0</sub> ≥ 5 × 10<sup>-3</sup></td><td>111282682<br>86634283</td><td>19778972<br>47632141</td><td>—<br>77118</td><td>—<br>61735</td><td>—<br>78300108</td><td>—<br>43125</td><td>—<br>71884333</td><td>—<br>93187</td><td>—<br>96161515</td><td>—<br>1853</td><td>—<br>8242610</td><td>—<br>102602<br>52715219</td></tr><tr><td>n<sub>0</sub> ≥ 1 × 10<sup>-3</sup></td><td>67448647<br>52859297</td><td>18465778<br>43750239</td><td>—<br>272353</td><td>—<br>312</td><td>—<br>225334648</td><td>—<br>3836927</td><td>—<br>1863218</td><td>—<br>1438469</td><td>—<br>1071366</td><td>—<br>1846182</td><td>—<br>1551939</td><td>—<br>400555<br>127615892</td></tr></tbody></table>

<table><caption>TABLE II
FIRST TWENTY CHEBYSHEV SERIES COEFFICIENTS FOR THE INVERSE SQUARE ROOT FUNCTION FOR DIFFERENT $n_0$.</caption><thead><tr><th>n<sub>0</sub></th><th>c<sub>n</sub><sup>CPE</sup></th><th>c<sub>0</sub><sup>CPE</sup></th><th>c<sub>1</sub><sup>CPE</sup></th><th>c<sub>2</sub><sup>CPE</sup></th><th>c<sub>3</sub><sup>CPE</sup></th><th>c<sub>4</sub><sup>CPE</sup></th><th>c<sub>5</sub><sup>CPE</sup></th><th>c<sub>6</sub><sup>CPE</sup></th><th>c<sub>7</sub><sup>CPE</sup></th><th>c<sub>8</sub><sup>CPE</sup></th><th>c<sub>9</sub><sup>CPE</sup></th></tr><tr><th></th><th></th><th>c<sub>10</sub><sup>CPE</sup></th><th>c<sub>11</sub><sup>CPE</sup></th><th>c<sub>12</sub><sup>CPE</sup></th><th>c<sub>13</sub><sup>CPE</sup></th><th>c<sub>14</sub><sup>CPE</sup></th><th>c<sub>15</sub><sup>CPE</sup></th><th>c<sub>16</sub><sup>CPE</sup></th><th>c<sub>17</sub><sup>CPE</sup></th><th>c<sub>18</sub><sup>CPE</sup></th><th>c<sub>19</sub><sup>CPE</sup></th></tr></thead><tbody><tr><td>n<sub>0</sub> ≥ 1 × 10<sup>-1</sup></td><td>374048017<br>113951175</td><td>88499941<br>99875048</td><td>—<br>17881</td><td>15797788<br>41555977</td><td>—<br>23207</td><td>—<br>154174661</td><td>—<br>1821274</td><td>1646349<br>105728006</td><td>—<br>87247730</td><td>187127<br>51011883</td><td>—<br>257569<br>142972536</td></tr><tr><td>n<sub>0</sub> ≥ 5 × 10<sup>-2</sup></td><td>249779459<br>67453088</td><td>350792999<br>280590374</td><td>—<br>491723</td><td>132037463<br>217146298</td><td>—<br>360646</td><td>—<br>17309909</td><td>—<br>11726340</td><td>9883564<br>162320651</td><td>—<br>9377515</td><td>1014741<br>47325860</td><td>—<br>16780163<br>1303748560</td></tr><tr><td>n<sub>0</sub> ≥ 1 × 10<sup>-2</sup></td><td>237852977<br>50542861</td><td>39264271<br>17952244</td><td>—<br>4967215</td><td>195265677<br>3834427</td><td>—<br>490771</td><td>—<br>10341695</td><td>—<br>81979913</td><td>100102199<br>16961213</td><td>—<br>93897642</td><td>19174519<br>80051251</td><td>—<br>10291450<br>55779169</td></tr><tr><td>n<sub>0</sub> ≥ 5 × 10<sup>-3</sup></td><td>564146905<br>109693414</td><td>271445188<br>103894561</td><td>—<br>2511703</td><td>203729517<br>112913626</td><td>—<br>67778415</td><td>—<br>84437949</td><td>—<br>49111051</td><td>248249039<br>1048651</td><td>—<br>133947555</td><td>200509969<br>19174519</td><td>—<br>1699069<br>56928095</td></tr><tr><td>n<sub>0</sub> ≥ 1 × 10<sup>-3</sup></td><td>67588883<br>248678426</td><td>—<br>9946203</td><td>—<br>16184137</td><td>—<br>21196048</td><td>—<br>76778415</td><td>—<br>120119399</td><td>—<br>17846601</td><td>301832257<br>1020511615</td><td>—<br>237360333</td><td>871395101<br>117123</td><td>—<br>182008350<br>62114496</td></tr><tr><td>n<sub>0</sub> ≥ 1 × 10<sup>-3</sup></td><td>630048624<br>102215551</td><td>—<br>94808470</td><td>—<br>95599410</td><td>—<br>184341026</td><td>—<br>55931467</td><td>—<br>80408470</td><td>—<br>494742647</td><td>—<br>1224960907</td><td>—<br>130317757</td><td>109757935<br>122418745</td><td>—<br>36868111<br>35375072</td></tr><tr><td></td><td>65195223<br>69543947</td><td>—<br>112964291</td><td>—<br>201622718</td><td>—<br>184059003</td><td>—<br>55931467</td><td>—<br>68082220</td><td>—<br>262965811</td><td>—<br>122566312</td><td>—<br>99758659</td><td>102912463<br>53675438</td><td>—<br>34220975<br>84624381</td></tr></tbody></table>

Here, $N_A$ is the expansion order of PAE, and $c_n^{\text{PAE}}$ is the $n^{\text{th}}$ expansion coefficient of PAE for the approximations of $\sqrt{x}$ and $\sqrt{x}^{-1}$, which depends on $N_A$ and can be obtained after solving the linear equation systems as presented in [25].

The numerical expressions of their matrix function counterparts are obtained by generalizing (19) to the matrix arguments as

$$ \mathbf{X} \mapsto \begin{cases} \sqrt{\mathbf{X}} \approx (\sum_{n=0}^{N_A} c_n^{\text{PAE}} \mathbf{X}^{N_A-n})^{-1} (\sum_{n=0}^{N_A} c_n^{\text{PAE}} \mathbf{X}^n) \\ \sqrt{\mathbf{X}}^{-1} \approx (\sum_{n=0}^{N_A} c_n^{\text{PAE}} \mathbf{X}^n)^{-1} (\sum_{n=0}^{N_A} c_n^{\text{PAE}} \mathbf{X}^{N_A-n}) \end{cases}. \quad (20) $$

The detailed values of the expansion coefficients $c_n^{\text{PAE}}$ of PAE for $N_A = 0, 1, \dots, 9$ are provided in Table III.

### D. General Comments

Before the presentation of numerical experiments in Section IV, we would like to attract the readers' attention on the following points.

For the reader's convenience, Table IV lists the order $N_T$, $N_C$, and $N_A$ at which the Taylor, Chebyshev, Padé series should be truncated for the numerical expansion of $\sqrt{\mathbf{X}}$ and $\sqrt{\mathbf{X}}^{-1}$ to reach a target relative error $\delta$ for $\mathbf{X}$ with different condition number (different $n_0$); these truncation orders have been obtained through extensive numerical experiments.

## IV. NUMERICAL RESULTS

In this section, some numerical experiments are presented to investigate the advantages and disadvantages of TSE, CPE, and PAE for the

<table><caption>TABLE III
FIRST TEN PADÉ SERIES COEFFICIENTS FOR DIFFERENT $N_A$.</caption><thead><tr><th>c<sub>n</sub><sup>CPE</sup></th><th>N<sub>A</sub></th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th></tr></thead><tbody><tr><td>c<sub>0</sub><sup>PAE</sup></td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>c<sub>1</sub><sup>PAE</sup></td><td>—</td><td>3</td><td>10</td><td>21</td><td>36</td><td>55</td><td>78</td><td>105</td><td>136</td><td>171</td><td>171</td></tr><tr><td>c<sub>2</sub><sup>PAE</sup></td><td>—</td><td>—</td><td>5</td><td>35</td><td>126</td><td>330</td><td>715</td><td>1365</td><td>2380</td><td>3876</td><td>3876</td></tr><tr><td>c<sub>3</sub><sup>PAE</sup></td><td>—</td><td>—</td><td>—</td><td>7</td><td>84</td><td>462</td><td>1716</td><td>5005</td><td>12376</td><td>27132</td><td>27132</td></tr><tr><td>c<sub>4</sub><sup>PAE</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>9</td><td>165</td><td>1287</td><td>6435</td><td>24310</td><td>75582</td><td>75582</td></tr><tr><td>c<sub>5</sub><sup>PAE</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>11</td><td>286</td><td>3003</td><td>19448</td><td>92378</td><td>92378</td></tr><tr><td>c<sub>6</sub><sup>PAE</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>13</td><td>455</td><td>6188</td><td>50388</td><td>50388</td></tr><tr><td>c<sub>7</sub><sup>PAE</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>15</td><td>680</td><td>11628</td><td>11628</td></tr><tr><td>c<sub>8</sub><sup>PAE</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>17</td><td>969</td><td>969</td></tr><tr><td>c<sub>9</sub><sup>PAE</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>19</td><td>19</td></tr></tbody></table>

computation of $\sqrt{\mathbf{G}}$ and $\sqrt{\mathbf{G}}^{-1}$ for $\mathbf{G} \in \{\mathbf{G}_{\text{f,f}}, \mathbf{G}_{\text{g,g}}, \mathbf{G}_{\lambda,\lambda}, \mathbf{G}_{\bar{\lambda},\bar{\lambda}}\}$ with thorough performance analysis and demonstrate the applicability of the proposed techniques to spectral analysis.

### A. Accuracy of the Computation of $\sqrt{\mathbf{G}}$ and $\sqrt{\mathbf{G}}^{-1}$

As the first numerical experiment, the accuracy of the computation of $\sqrt{\mathbf{G}}$ and $\sqrt{\mathbf{G}}^{-1}$ using the proposed techniques is investigated. A PEC sphere of radius 0.5 m centered at the origin and residing in free space is considered. The sphere surface is non-uniformly discretized into three meshes designated as “mesh-1”, “mesh-2”, and “mesh-3” with $\{N_E, N_P, N_V\} = \{1230, 820, 412\}$, $\{N_E, N_P, N_V\} =$
<table><caption>TABLE IV
TRUNCATION ORDER $N_T$, $N_C$, AND $N_A$ FOR THE EXPANSION OF $\sqrt{\mathbf{X}}$ (LEFT) AND $\sqrt{\mathbf{X}}^{-1}$ (RIGHT) TO REACH TARGET $\delta$ FOR $\mathbf{X}$ WITH DIFFERENT $n_0$.</caption><thead><tr><th>n<sub>0</sub></th><th></th><th>10<sup>-2</sup></th><th>10<sup>-3</sup></th><th>10<sup>-4</sup></th><th>10<sup>-5</sup></th><th>10<sup>-6</sup></th></tr></thead><tbody><tr><td rowspan="3">n<sub>0</sub> ≥ 1 × 10<sup>-1</sup></td><td>N<sub>T</sub></td><td>12</td><td>25</td><td>41</td><td>59</td><td>—</td></tr><tr><td>N<sub>C</sub></td><td>3</td><td>5</td><td>8</td><td>11</td><td>14</td></tr><tr><td>N<sub>A</sub></td><td>3</td><td>5</td><td>7</td><td>8</td><td>10</td></tr><tr><td rowspan="3">n<sub>0</sub> ≥ 5 × 10<sup>-2</sup></td><td>N<sub>T</sub></td><td>21</td><td>47</td><td>—</td><td>—</td><td>—</td></tr><tr><td>N<sub>C</sub></td><td>4</td><td>7</td><td>11</td><td>15</td><td>19</td></tr><tr><td>N<sub>A</sub></td><td>4</td><td>7</td><td>9</td><td>12</td><td>14</td></tr><tr><td rowspan="3">n<sub>0</sub> ≥ 1 × 10<sup>-2</sup></td><td>N<sub>T</sub></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>N<sub>C</sub></td><td>6</td><td>13</td><td>21</td><td>31</td><td>40</td></tr><tr><td>N<sub>A</sub></td><td>8</td><td>13</td><td>19</td><td>—</td><td>—</td></tr><tr><td rowspan="3">n<sub>0</sub> ≥ 5 × 10<sup>-3</sup></td><td>N<sub>T</sub></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>N<sub>C</sub></td><td>8</td><td>17</td><td>29</td><td>41</td><td>—</td></tr><tr><td>N<sub>A</sub></td><td>10</td><td>18</td><td>—</td><td>—</td><td>—</td></tr><tr><td rowspan="3">n<sub>0</sub> ≥ 1 × 10<sup>-3</sup></td><td>N<sub>T</sub></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>N<sub>C</sub></td><td>12</td><td>30</td><td>—</td><td>—</td><td>—</td></tr><tr><td>N<sub>A</sub></td><td>16</td><td>—</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>

<table><thead><tr><th>n<sub>0</sub></th><th></th><th>10<sup>-2</sup></th><th>10<sup>-3</sup></th><th>10<sup>-4</sup></th><th>10<sup>-5</sup></th><th>10<sup>-6</sup></th></tr></thead><tbody><tr><td rowspan="3">n<sub>0</sub> ≥ 1 × 10<sup>-1</sup></td><td>N<sub>T</sub></td><td>31</td><td>51</td><td>72</td><td>—</td><td>—</td></tr><tr><td>N<sub>C</sub></td><td>5</td><td>9</td><td>12</td><td>15</td><td>19</td></tr><tr><td>N<sub>A</sub></td><td>4</td><td>6</td><td>8</td><td>9</td><td>11</td></tr><tr><td rowspan="3">n<sub>0</sub> ≥ 5 × 10<sup>-2</sup></td><td>N<sub>T</sub></td><td>64</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>N<sub>C</sub></td><td>8</td><td>13</td><td>17</td><td>22</td><td>27</td></tr><tr><td>N<sub>A</sub></td><td>6</td><td>8</td><td>11</td><td>13</td><td>16</td></tr><tr><td rowspan="3">n<sub>0</sub> ≥ 1 × 10<sup>-2</sup></td><td>N<sub>T</sub></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>N<sub>C</sub></td><td>18</td><td>28</td><td>39</td><td>—</td><td>—</td></tr><tr><td>N<sub>A</sub></td><td>13</td><td>19</td><td>—</td><td>—</td><td>—</td></tr><tr><td rowspan="3">n<sub>0</sub> ≥ 5 × 10<sup>-3</sup></td><td>N<sub>T</sub></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>N<sub>C</sub></td><td>25</td><td>40</td><td>—</td><td>—</td><td>—</td></tr><tr><td>N<sub>A</sub></td><td>19</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td rowspan="3">n<sub>0</sub> ≥ 1 × 10<sup>-3</sup></td><td>N<sub>T</sub></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>N<sub>C</sub></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>N<sub>A</sub></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Fig. 1. Relative error $\delta$ of the computation of (a) $\sqrt{G_{f,f}}$, (b) $\sqrt{G_{g,g}}$, (c) $\sqrt{G_{\lambda,\bar{\lambda}}}$, and (d) $\sqrt{G_{\bar{\lambda},\bar{\lambda}}}$ using TSE, CPE-1, CPE-2, and PAE for three mesh sets against the expansion order.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Relative error $\delta$ of the computation of (a) $\sqrt{G_{f,f}}^{-1}$, (b) $\sqrt{G_{g,g}}^{-1}$, (c) $\sqrt{G_{\lambda,\bar{\lambda}}^{-1}}$, and (d) $\sqrt{G_{\bar{\lambda},\bar{\lambda}}^{-1}}$ using TSE, CPE-1, CPE-2, and PAE for three mesh sets against the expansion order.</figcaption></figure>

$\{2193, 1462, 733\}$, and $\{N_E, N_P, N_V\} = \{4449, 2966, 1485\}$, respectively. These sets of meshes are selected so that the resulting Gram matrices are increasingly ill-conditioned for mesh-1, mesh-2, and mesh-3 with the condition number of $\{G_{f,f}, G_{g,g}, G_{\lambda,\lambda}, G_{\bar{\lambda},\bar{\lambda}}\}$ equal to $\{5.16, 5.74, 6.72, 9.11\}$, $\{20.45, 14.49, 29.12, 24.58\}$, and $\{615.61, 344.95, 361.91, 269.31\}$, respectively.

Fig. 1(a)–(d) and Fig. 2(a)–(d) show the relative error $\delta$ of the computation of $\sqrt{G_{f,f}}$, $\sqrt{G_{g,g}}$, $\sqrt{G_{\lambda,\lambda}}$, $\sqrt{G_{\bar{\lambda},\bar{\lambda}}}$ and $\sqrt{G_{f,f}}^{-1}$, $\sqrt{G_{g,g}}^{-1}$, $\sqrt{G_{\lambda,\lambda}}^{-1}$, $\sqrt{G_{\bar{\lambda},\bar{\lambda}}^{-1}}$ for three sets of meshes against the expansion order, respectively, where $\delta$ is the relative error between $\sqrt{G}$ or $\sqrt{G}^{-1}$ computed using SVD and its numerical computation as described in Section III. As expected, the accuracy of the computation of $\sqrt{G}$ and $\sqrt{G}^{-1}$ using the proposed techniques for all the sets of meshes increases with the expansion order. Moreover, for a fixed expansion order and a given algorithm, the accuracy of the numerical computation of $\sqrt{G}$ and $\sqrt{G}^{-1}$ improves with mesh uniformity. This is due to the fact that, as the mesh quality improves, $G$ exhibits a lower condition number approaching 1 and

the numerical computation of $\sqrt{G}$ and $\sqrt{G}^{-1}$ improves too (see Section III-D).

Moreover, for mesh-1 (resulting in lower condition number for $G$), the PAE converges faster than the TSE and the CPE. In this case, one could choose the PAE that needs fewer expansion order for a given computational accuracy.

In addition, the figures show that we could use CPE-2 as a good approximation of CPE-1 by directly using the bounds of the known ranges in which $n_0$ lives and the corresponding coefficients in Table I and II for CPE for convenience. As one would expect, the error behaviour of CPE-1 and CPE-2 converges towards one another when the actual value of $n_0$ for a given set of meshes and a given Gram matrix corresponds to the lower bound of the range of $n_0$ (used to select the CPE-2 tabulated coefficients). Additionally, the Chebyshev series truncation order $N_C$ of CPE-2 for the computation of $\sqrt{G}$ and $\sqrt{G}^{-1}$ to reach target $\delta$ for $G$ with different condition number in Fig. 1(a)–(d) and Fig. 2(a)–(d) generally matches with that reported in Table IV. Also, the accuracy of the computation of $\sqrt{G}^{-1}$ is
<figure><img src="image_4.png" /><figcaption>Fig. 3. Comparison of the singular values of $\mathbf{T}$ and $\tilde{\mathbf{T}}$ with the spectrum of $\mathcal{T}$ including a zoom around singular value index 3380 (a) for $k = 0.1 \text{ rad m}^{-1}$ and (b) for $k = 2.744 \text{ rad m}^{-1}$.</figcaption></figure>

generally lower than that of $\sqrt{\mathbf{G}}$ using the same expansion order, the same algorithm, and the same set of meshes for the same $\mathbf{G}$ in Fig. 1(a)–(d) and Fig. 2(a)–(d), which is in accordance with the reports in Table IV.

## B. Application to Spectral Analysis

As the second numerical experiment, the use of the proposed techniques as a tool for the study, manipulation, and accurate recovery of the spectrum of the relevant integral operators is showcased. The surface of a PEC sphere of radius 1 m is uniformly discretized with $\{N_E, N_P, N_V\} = \{5070, 3380, 1692\}$ and the corresponding boundary element matrix $\mathbf{T}$ is obtained as delineated in Section II. The spectral analysis is demonstrated for a spurious resonance of the EFIE operator and for a non-resonant wavenumber.

Fig. 3(a) and (b) compare the singular values of the standard EFIE matrix $\mathbf{T}$ and its normalized counterpart $\tilde{\mathbf{T}}$ with the analytic spectrum of the surface integral operator $\mathcal{T}$ obtained from a spherical harmonics analysis of the problem [17]. The non-resonant and resonant wavenumbers [26] are respectively $k = 0.1 \text{ rad m}^{-1}$ and $k = 2.744 \text{ rad m}^{-1}$ and the condition numbers of the corresponding $\mathbf{T}$ matrices are $1.7 \times 10^5$ and $2.9 \times 10^3$. The results are obtained with the PAE for $\sqrt{\mathbf{G}_{\mathbf{f},\mathbf{f}^{-1}}}$ with $N_A = 9$. Clearly, the singular values of $\tilde{\mathbf{T}}$ have a good agreement with the spectrum of $\mathcal{T}$ while those of $\mathbf{T}$ do not, even for the resonant scenario.

When assessing the quality of the spectrum obtained from the BEM matrices of the EM operators, one should keep in mind that several sources of errors will cause deviation of the matrix spectrum from the continuous operator's spectrum. They include the error in the computation of the elements of the BEM matrices—for instance in the computation of near singular integrals [27]—, the error caused by the truncation of the expansions of the reciprocal square roots of the Gram matrices, and the error introduced by the spectral properties of the discretization bases [28]. In Fig. 3, these errors degrade the precision of the singular values associated to highly oscillating singular vectors while the rest of the spectrum shows a good agreement with spherical harmonics.

## V. CONCLUSION

In this work, three polynomial-expansion based algorithms are delineated for the computation of the square roots and inverse square roots of some of the most relevant Gram matrices for SIEs. Properly chosen polynomial expansions of the square root and inverse square root of the scalar functions and the theory of matrix functions are used. Tables containing the different sets of the expansion coefficients are provided along with the comparative numerical experiments that illustrate the advantages and disadvantages of the different strategies. Moreover, the application of the proposed techniques to spectral analysis is demonstrated, which could be applied for most of the standard SIEs in electromagnetics and acoustics.

## REFERENCES

[1] J.-M. Jin, *Theory and Computation of Electromagnetic Fields*. 2nd ed. New York, NY, USA: Wiley, 2015.

[2] D. Boffi, “Compatible discretizations for eigenvalue problems,” in *Compatible Spatial Discretizations*, IMA Volumes in Mathematics and its Applications, vol. 142, pp. 121–142, Berlin: Springer, 2006.

[3] A. Merlini, C. Henry, D. Consoli, L. Rahmouni, A. Dély and F. P. Andriulli, “Laplacian filtered loop-star decompositions and quasi-Helmholtz filters: Definitions, analysis, and efficient algorithms,” *IEEE Trans. Antennas Propag.*, vol. 71, no. 12, pp. 9289–9302, Dec. 2023.

[4] R. A. Horn and C. R. Johnson, *Matrix Analysis*. Cambridge University Press, 2012.

[5] A. F. Peterson, C. F. Smith, and R. Mittra, “Eigenvalues of the moment-method matrix and their effect on the convergence of the conjugate gradient algorithm,” *IEEE Trans. Antennas Propag.*, vol. 36, no. 8, pp. 1177–1179, Aug. 1988.

[6] A. Bjorck and S. Hammarling, “A Schur method for the square root of a matrix,” *Linear Algebra Appl.*, vol. 52–53, pp. 127–140, Jul. 1983.

[7] N. J. Higham, “Newton’s method for the matrix square root,” *Math. Comp.*, vol. 46, no. 174, pp. 537–549, Apr. 1986.

[8] E. D. Denman and A. N. Beavers, “The matrix sign function and computations in systems,” *Applied Mathematics and Computation*, vol. 2, no. 1, pp. 63–94, Jan. 1976.

[9] G. W. Stewart, “A Krylov-Schur algorithm for large eigenproblems,” *SIAM J. Matrix Anal. Appl.*, vol. 23, issue 3, pp. 601–614, 2001.

[10] N. J. Higham, *Functions of Matrices: Theory and Computation*. Philadelphia, PA, USA: SIAM, 2008.

[11] Y. Sun, S. Luo, and X. Lei, “Gram matrices of mixed-state ensembles,” *Int. J. Theor. Phys.*, vol. 60, pp. 3211–3224, 2021.

[12] R. Chen, A. Merlini, and F. P. Andriulli, “Calculation of inverse square root of RWG Gram matrix based on Chebyshev polynomial expansion,” in *Proc. Int. Appl. Comput. Electromagn. Soc. Symp. (ACES-China)*, Xi’an, China, Aug. 2024, pp. 1–2.

[13] S. Rao, D. Wilton, and A. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” *IEEE Trans. Antennas Propag.*, vol. AP-30, no. 3, pp. 409–418, May 1982.

[14] F. P. Andriulli, K. Cools, H. Bagci, F. Olyslager, A. Buffa, S. Christiansen, and E. Michielssen, “A multiplicative Calderon preconditioner for the electric field integral equation,” *IEEE Trans. Antennas Propag.*, vol. 56, no. 8, pp. 2398–2412, Aug. 2008.

[15] A. Buffa and S. Christiansen, “A dual finite element complex on the barycentric refinement,” *Math. Comput.*, vol. 76, pp. 1743–1769, 2007.

[16] S. Y. Chen, W. C. Chew, J. M. Song, and J.-S. Zhao, “Analysis of low frequency scattering from penetrable scatterers,” *IEEE Trans. Geosci. Remote Sens.*, vol. 39, no. 4, pp. 726–735, Apr. 2001.

[17] S. B. Adrian, A. Dély, D. Consoli, A. Merlini, and F. P. Andriulli, “Electromagnetic integral equations: Insights in conditioning and preconditioning,” *IEEE Open J. Antennas Propag.*, vol. 2, pp. 1143–1174, Oct. 2021.

[18] C. R. Johnson, K. Okubo, and R. Reams, “Uniqueness of matrix square roots and an application,” *Linear Algebra and its Applications*, vol. 323, no. 1–3, pp. 51–60, Jan. 2001.

[19] G. H. Golub, and C. F. Van Loan, *Matrix Computations*. Johns Hopkins University Press, 2013.

[20] B. Taylor, *Methodus Incrementorum Directa et Inversa*. London, UK: Impensis Gulielmi Innyns, 1717.

[21] J. C. Mason and C. D. Handscomb, *Chebyshev Polynomials*. London, UK: Chapman & Hall, 2002.

[22] F. W. J. Olver, *NIST Handbook of Mathematical Functions*. Cambridge, UK: Cambridge University Press, 2010.

[23] C. W. Clenshaw, “A note on the summation of Chebyshev series,” *Mathematics of Computation*, vol. 9, no. 51, pp. 118–120, 1955.

[24] G. A. Baker and P. Graves-Morris, *Padé Approximants*. Cambridge, UK: Cambridge University Press, 1996.

[25] Y. Song, N. Sebe and W. Wang, “Fast differentiable matrix square root and inverse square root,” *IEEE Trans. Pattern Anal. Mach. Intell.*, vol. 45, no. 6, pp. 7367–7380, Jun. 2023.

[26] R. F. Harrington, “Spherical wave functions,” in *Time-Harmonic Electromagnetic Fields*. Piscataway, NJ, USA: Wiley-IEEE Press, 2001, ch. 6, pp. 264–316.

[27] M. M. Botha, “A family of augmented Duffy transformations for near-singularity cancellation quadrature,” *IEEE Trans. Antennas Propag.*, vol. 61, no. 6, pp. 3123–3134, Jun. 2013.

[28] Karl F. Warnick, *Numerical Analysis for Electromagnetic Integral Equations*. Artech House, 2008.