# Reduced Order Modeling for Parameterized Electromagnetic Simulation Based on Tensor Decomposition

Xiao-Feng He, Liang Li, Stéphane Lanteri, Kun Li

**Abstract**—We present a data-driven surrogate modeling for parameterized electromagnetic simulation. This method extracts a set of reduced basis (RB) functions from full-order solutions through a two-step proper orthogonal decomposition (POD) method. A mapping from the time/parameter to the principal components of the projection coefficients, extracted by canonical polyadic decomposition (CPD), is approximated by a cubic spline interpolation (CSI) approach. The reduced-order model (ROM) is trained in the offline phase, while the RB solution of a new time/parameter value is recovered fast during the online phase. We evaluate the performance of the proposed method with numerical tests for the scattering of a plane wave by a 2-D multi-layer dielectric disk and a 3-D multi-layer dielectric sphere.

**Index Terms**—electromagnetic simulation, model order reduction, proper orthogonal decomposition, canonical polyadic decomposition, cubic spline interpolation

## I. INTRODUCTION

ELECTROMAGNETIC scattering problems are generally described as the parameterized time-domain Maxwell equations [16], where the parameters may represent geometric shapes, material properties, system configurations and so on. Solving the full-order model (FOM) repeatedly for varying parameter values by high-fidelity simulations is often required in many applications such as real-time system control, weather forecast and structure design, leading to great CPU time and memory demands. It is urgent to study the parametric surrogate modeling to accelerate the solution of the partial differential equations (PDEs) in an efficient and accurate way. Model order reduction (MOR) techniques have been proposed to solve this problem in the past decades. The main idea of MOR is to replace a high-dimensional FOM with a surrogate model, which has a lower dimensionality and acceptable accuracy loss [20].

Some techniques such as the fast multipole method (FMM) [29], subdomain multilevel approach (SMA) [30] and characteristic basis functions method (CBFM) [25] have been proposed to accelerate the solution of scattering problems. For example, CBFM has been developed for an efficient solution of Method of Moments matrix equations. Compared

to conventional RWG basis functions [28] bounded by the $\lambda/20$ domain discretization, this approach is based on high order basis functions over enlarged blocks, thus reducing the number of degrees of freedom (DoFs). A reduced basis (RB) method adopted in this work reduces the dimensionality of FOM in a completely different manner. In general, repeated simulations are usually not necessary because intrinsic similarities among the solutions of different parameters can be used to predict the solutions for new parameter values [11]. The RB method with an offline-online procedure, is one of the most widely used techniques for the surrogate modeling of time-dependent parameterized problems. In the offline phase, a reduced space representing the main dynamics of the FOM is constructed from full-order snapshots of different time/parameter values. There are several methods to extract the reduced basis functions, such as the proper orthogonal decomposition (POD) [14], [21], [22], greedy algorithms [18], [24] and an autoencoder [8], [10]. As the most widely-used technique, POD is used in this paper to generate the RB functions with a singular-value decomposition (SVD) [11]. As an approximation to the full-order solution, the RB solution, also named the reduced-order solution is recovered as a linear combination of the RB functions, where only the combination coefficients, also named projection coefficients will be calculated in the online phase [20]. It is advised that the methods mentioned above can be combined to accelerate the solution of parameterized electromagnetic scattering problems. For example, we could build a data-driven ROM based on snapshots generated by CBFM offline, then recover the RB solution for a new parameter online.

RB methods are divided into intrusive and non-intrusive (or data-driven) ones according to the methodology of calculating the projection coefficients. The former determines the projection coefficients by solving a low dimensional linear system, i.e., a projection of the full-order model onto the reduced space. And the Galerkin procedure is often used for the projection [7], [26]. However, the main drawback is that intrusive ROMs require access to the high dimensional model, and provide limited computational gain especially for complex nonlinear problems [1]. The non-intrusive ROMs approximate the projection coefficients over the parameter domain by using a database of reduced order information [5]. The high-fidelity solvers such as FEM, FDM, DGTD are only used to generate snapshots, leading to a full decoupling of the online stage and the high-fidelity scheme. Many excellent non-intrusive ROMs are proposed in the past. For instance, a feed-forward neural

The first two authors are with the School of Mathematical Sciences, University of Electronic Science and Technology of China (hexf@std.uestc.edu.cn, plum.liliang@gmail.com, plum_liliang@uestc.edu.cn). Stéphane Lanteri is with Inria ATLANTIS project-team (Stephane.Lanteri@inria.fr). Kun Li is with the School of Mathematics, Southwestern University of Finance and Economics (relikun@163.com, likun@swufe.edu.cn). The corresponding author is Liang Li.

Manuscript received March 25, 2023
network (FNN) has been designed for parameterized unsteady flows [15], [31], where a FNN is trained to approximate the projection coefficients. Interpolation and regression methods such as cubic spline interpolation (CSI) [3], [20], radial basis functions (RBF) [33] and Gaussian process regression (GPR) [11], [34] have also been adopted to approximate the projection coefficients. By combining an FNN-based regression model with a physics-informed neural network (PINN) [23], [27], a hybrid surrogate model is devised in [6], which requires less training data.

Data-driven reduced-order modeling based on RB method has been applied to electromagnetic scattering simulations [20], [34] recently. These related works use the POD technique to extract a set of time-/parameter-independent orthogonal basis functions and compute a series of projection coefficient matrices, then perform SVD for each projection matrix to extract the dominating temporal and parametric modes. Interpolation or regression models are built to approximate the temporal and parametric modes respectively. These methods are accurate enough for time-dependent parametric problems, but require to build a large number of interpolation or regression models (determined by the number of basis functions and the number of truncated singular values of each projection coefficient matrix), which increases the test time online. Motivated by the above problem, a ROM composed of two-step POD, canonical polyadic decomposition (CPD), and cubic spline interpolation (CSI) is developed for parameterized electromagnetic simulation. We employ a canonical polyadic decomposition (CPD) [17] to extract the principal components of the projection coefficients tensor, instead of repeated SVD for all projection coefficient matrices, which significantly reduces the number of interpolation or regression models, thus reducing the test time online with an acceptable loss of accuracy. The proposed method, named POD-CPD-CSI ROM is validated to be efficient and accurate for parameterized electromagnetic scattering problems.

The symbols and notations used in this paper are introduced as follows. Notations $t$ and $\mu$ denote time and parameter respectively. $S_u^j$ is the time trajectory matrix and $S_u$ is the snapshot matrix. $N_h$ is the geometric DoFs, $N_t$ and $N_p$ are the number of time and parameter samples respectively. $u_h(t, \mu)$ and $u_h^r(t, \mu)$ are the full-order solution and reduced-order solution respectively. $\alpha_u(t, \mu)$ denotes the projection coefficient vector and $V$ is the POD basis matrix. $P_u^l$ is the projection coefficient matrix and $G_u$ is the projection coefficient tensor.

## II. TWO-STEP POD METHOD

Assume that $\mathcal{P}_h^{tr} = \{\mu_1, \mu_2, \dots, \mu_{N_p}\}$ is a parameter sampling that we're interested in. For each $\mu_j \in \mathcal{P}_h^{tr}$, we have a time trajectory matrix

$$ \mathbf{S}_{\mathrm{u}}^{j} = [\mathbf{u}_{h}(t_{1}, \mu_{j}) | \mathbf{u}_{h}(t_{2}, \mu_{j}) | \cdots | \mathbf{u}_{h}(t_{N_{t}}, \mu_{j})] \in \mathbb{R}^{N_{h} \times N_{t}}, \quad (1) $$

with $\mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}$, which could be obtained either from high-fidelity solutions or from experiments, for a time sampling $\mathcal{T}_h^{tr} = \{t_1, t_2, \dots, t_{N_t}\}$. Here, $N_h$ denotes the number of

geometric DoFs. A snapshot matrix for all parameters in $\mathcal{P}_h^{tr}$ is then formed

$$ \mathbf{S}_{\mathrm{u}} = [\mathbf{S}_{\mathrm{u}}^{1} | \mathbf{S}_{\mathrm{u}}^{2} | \cdots | \mathbf{S}_{\mathrm{u}}^{N_{p}}] \in \mathbb{R}^{N_{h} \times N_{s}}, \quad (2) $$

with $N_s = N_t \cdot N_p$. To reduce the dimensionality of FOM, we seek a low-rank decomposition to the snapshot matrix and construct a reduced space $V_{u,rb} = \text{span}\{\mathbf{v}_{u,1}, \mathbf{v}_{u,2}, \dots, \mathbf{v}_{u,d_u}\}$ of dimension $d_u \ll \min\{N_h, N_s\}$, where $\{\mathbf{v}_{u,i}\}_{i=1}^{d_u}$ are time- and parameter-independent orthogonal RB functions [11], [34].

For the large snapshots matrix $S_u$, i.e. when $N_s$ is large, it will be very computationally expensive to directly perform POD on $S_u$. So, we adopt a two-step POD method [20] shown in Algorithm 1 to reduce the demands on computational cost. Details of the computational complexity and error analysis of this method can be found in [20]. The RB solution $\mathbf{u}_h^r(t, \mu)$, an approximation to the full-order solution $\mathbf{u}_h(t, \mu)$, is recovered as

$$ \mathbf{u}_h^r(t, \mu) = \mathbf{V}_u \alpha_u(t, \mu) = \sum_{i=1}^{d_u} \alpha_{u,i}(t, \mu) \mathbf{v}_{u,i}, \quad (3) $$

where $\mathbf{V}_u \in \mathbb{R}^{N_h \times d_u}$ is the POD basis matrix generated by the two-step POD and vector $\alpha_u(t, \mu) = [\alpha_{u,1}(t, \mu), \dots, \alpha_{u,d_u}(t, \mu)]^T$ is called the coefficient vector [34]. After calculating the POD basis functions, we need to build a mapping between the time/parameters and the projection coefficients.

### Algorithm 1: Two-step POD approach

input: $\mathbf{S}_u^1, \mathbf{S}_u^2, \dots, \mathbf{S}_u^{N_p}$ and truncation parameters $\rho_t$ and $\rho_\mu$

output: POD basis matrix $\mathbf{V}_u$

1. Calculate matrices $\mathbf{T}_u^j = \text{POD}(\mathbf{S}_u^j, \rho_t)$ for $j = 1, 2, \dots, N_p$;

2. Assemble matrix $\mathbf{T}_u = [\mathbf{T}_u^1 | \mathbf{T}_u^2 | \dots | \mathbf{T}_u^{N_p}]$;

3. Compute matrix $\mathbf{V}_u = \text{POD}(\mathbf{T}_u, \rho_\mu)$.

4. Function $\mathbf{V} = \text{POD}(\mathbf{A}, \rho)$

5. Compute the eigenvalues $\{\lambda_i\}_{i=1}^r$ and the eigenvectors $\{\mathbf{u}_i\}_{i=1}^r$ of matrix $\mathbf{A}^T\mathbf{A}$ with $r = \text{rank}(\mathbf{A})$;

6. Determine the number of basis functions $d$ based on $d = \arg\min\{\pi(d) : \pi(d) \ge 1 - \rho\}$ and $\pi(d) = \sum_{i=1}^d \lambda_i / \sum_{i=1}^r \lambda_i$;

7. Calculate the RB functions $\mathbf{v}_i = \frac{1}{\sqrt{\lambda_i}}\mathbf{A}\mathbf{u}_i$ for $i = 1, 2, \dots, d, d \ll r$;

8. Assemble matrix $\mathbf{V} = [\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_d]$.

9. end

## III. CP-DECOMPOSITION-BASED INTERPOLATION

Each vector in the snapshot matrix (2) can be approximated by a linear combination of the RB basis vectors (3) [34].
In such a way, we create training data and rewrite all the projection coefficients with a matrix form

$$
\begin{align}
\mathbf{P}_{\mathbf{u}}^{\ell} &= \begin{pmatrix}
\alpha_{\mathbf{u},\ell}(t_1, \boldsymbol{\mu}_1) & \alpha_{\mathbf{u},\ell}(t_1, \boldsymbol{\mu}_2) & \cdots & \alpha_{\mathbf{u},\ell}(t_1, \boldsymbol{\mu}_{N_p}) \\
\alpha_{\mathbf{u},\ell}(t_2, \boldsymbol{\mu}_1) & \alpha_{\mathbf{u},\ell}(t_2, \boldsymbol{\mu}_2) & \cdots & \alpha_{\mathbf{u},\ell}(t_2, \boldsymbol{\mu}_{N_p}) \\
\vdots & \vdots & \ddots & \vdots \\
\alpha_{\mathbf{u},\ell}(t_{N_t}, \boldsymbol{\mu}_1) & \alpha_{\mathbf{u},\ell}(t_{N_t}, \boldsymbol{\mu}_2) & \cdots & \alpha_{\mathbf{u},\ell}(t_{N_t}, \boldsymbol{\mu}_{N_p})
\end{pmatrix} \nonumber \\
&= [\mathbf{v}_{\mathbf{u},\ell}^{\mathrm{T}} \mathbf{S}_{\mathbf{u}}^{\mathrm{1}}; \mathbf{v}_{\mathbf{u},\ell}^{\mathrm{T}} \mathbf{S}_{\mathbf{u}}^{\mathrm{2}}; \cdots; \mathbf{v}_{\mathbf{u},\ell}^{\mathrm{T}} \mathbf{S}_{\mathbf{u}}^{\mathrm{N_p}}]^{\mathrm{T}} \in \mathbb{R}^{N_t \times N_p}, \tag{4}
\end{align}
$$

with $\ell = 1, 2, \dots, d_u$, $\mathbf{u} \in \{\mathbf{E}, \mathbf{H}\}$. Instead of dealing with matrix $\mathbf{P}_{\mathbf{u}}^{\ell}$ directly [20], [34], we concatenate the coefficient matrices to form a coefficient tensor $\mathcal{G}_{\mathbf{u}}$ of dimension $(N_t, N_p, d_{\mathbf{u}})$. The CP decomposition [17] illustrated in Fig. 1 factorizes $\mathcal{G}_{\mathbf{u}}$ into a sum of component rank-one tensors

$$
\mathcal{G}_{\mathbf{u}} \approx \hat{\mathcal{G}}_{\mathbf{u}} = \sum_{r=1}^{R_u} \varphi_{\mathbf{u},r} \circ \psi_{\mathbf{u},r} \circ \zeta_{\mathbf{u},r}, \quad (5)
$$

where $R_u$ is a positive integer and determines the error of
CP decomposition. The symbol $\circ$ represents the vector outer
product and $\varphi_{\mathbf{u},r} \in \mathbb{R}^{N_t}$, $\psi_{\mathbf{u},r} \in \mathbb{R}^{N_p}$, and $\zeta_{\mathbf{u},r} \in \mathbb{R}^{d_{\mathbf{u}}}$ for
$r = 1, \dots, R_u$ are vectors. Our goal is to compute a CP

<figure><img src="image_1.png" /><figcaption>Fig. 1. CP decomposition of a three-order tensor $\mathcal{G}_{\mathbf{u}}$.</figcaption></figure>

decomposition with $R_u$ components that best approximates
$\mathcal{G}_u$, i.e., to find

$$
\min_{\hat{\mathcal{G}}_u} \|\mathcal{G}_u - \hat{\mathcal{G}}_u\| \quad \text{with} \quad \hat{\mathcal{G}}_u = \sum_{r=1}^{R_u} \varphi_{\mathbf{u},r} \circ \psi_{\mathbf{u},r} \circ \zeta_{\mathbf{u},r}. \quad (6)
$$

We adopt an alternating least squares (ALS) method [4], [12]
restated in Algorithm 2 to find the CP decomposition in
this work, where $\mathcal{G}_{\mathbf{u}(n)}$ is the mode-$n$ unfolding of tensor
$\mathcal{G}$, and notation $\odot$ and $\oplus$ denote the Khatri-Rao product
and Hadamard product respectively. $\mathbf{A}^\dagger$ denotes the Moore-
Penrose pseudoinverse of $\mathbf{A}$. Readers are referred to [17] for
more details.

With the discrete data set, a cubic spline interpolation-based
model $\hat{\omega}_{\mathbf{u}}$ is trained to approximate the continuous modes as
$$
\begin{align}
t &\mapsto \hat{\varphi}_{\mathbf{u},r}(t) \quad \text{with} \quad \{(t_i, (\varphi_{\mathbf{u},r})_i), i = 1, 2, \dots, N_t\}, \nonumber \\
\boldsymbol{\mu} &\mapsto \hat{\psi}_{\mathbf{u},r}(\boldsymbol{\mu}) \quad \text{with} \quad \{(\boldsymbol{\mu}_j, (\psi_{\mathbf{u},r})_j), j = 1, 2, \dots, N_p\}, \tag{7}
\end{align}
$$

for $r = 1, 2, \dots, R_u$. Elementwise, we have

$$
(\mathcal{G}_{\mathbf{u}})_{ijk} \approx (\hat{\mathcal{G}}_{\mathbf{u}})_{ijk} \approx \sum_{r=1}^{R_u} \hat{\varphi}_{\mathbf{u},r}(t_i) \hat{\psi}_{\mathbf{u},r}(\boldsymbol{\mu}_j)(\zeta_{\mathbf{u},r})_k, \quad (8)
$$

with $1 \le i \le N_t$, $1 \le j \le N_p$, $1 \le k \le d_{\mathbf{u}}$. For a new
time/parameter value $(t^*, \boldsymbol{\mu}^*)$, the RB solution is recovered as

$$
\mathbf{u}_{h}^{r}(t^{*}, \boldsymbol{\mu}^{*}) = \mathbf{V}_{\mathbf{u}} \hat{\boldsymbol{\alpha}}_{\mathbf{u}}(t^{*}, \boldsymbol{\mu}^{*}) = \mathbf{V}_{\mathbf{u}} \sum_{r=1}^{R_u} \hat{\varphi}_{\mathbf{u},r}(t^{*}) \hat{\psi}_{\mathbf{u},r}(\boldsymbol{\mu}^{*}) \zeta_{\mathbf{u},r}. \quad (9)
$$

**Algorithm 2:** ALS method for CP decomposition of the projection coefficient tensor $\mathcal{G}_u$

input : projection coefficient tensor $\mathcal{G}_u \in \mathbb{R}^{N_t \times N_p \times d_{\mathbf{u}}}$ and a truncation parameter $R_u$.

output: factor matrices $\Phi_u = [\varphi_{u,1}, \dots, \varphi_{u,R_u}]$, $\Psi_u = [\psi_{u,1}, \dots, \psi_{u,R_u}]$ and $\Sigma_u = [\zeta_{u,1}, \dots, \zeta_{u,R_u}]$.

1. Initialize $\Phi_u \in \mathbb{R}^{N_t \times R_u}$, $\Psi_u \in \mathbb{R}^{N_p \times R_u}$ and $\Sigma_u \in \mathbb{R}^{d_{\mathbf{u}} \times R_u}$;

2 repeat

3. $\Phi_u = \mathcal{G}_{u(1)}(\Sigma_u \odot \Psi_u)(\Sigma_u^T \Sigma_u \oplus \Psi_u^T \Psi_u)^\dagger$

4. $\Psi_u = \mathcal{G}_{u(2)}(\Phi_u \odot \Sigma_u)(\Phi_u^T \Phi_u \oplus \Sigma_u^T \Sigma_u)^\dagger$

5. $\Sigma_u = \mathcal{G}_{u(3)}(\Psi_u \odot \Phi_u)(\Psi_u^T \Psi_u \oplus \Phi_u^T \Phi_u)^\dagger$

6. **until** fit cases to improve or maximum iterations exhausted;

The procedure of POD-CP ROM is decoupled into of-
line/online phases, which is shown in Algorithm3, Fig. 2
and Fig. 3. Compared to previous works [20], [34],fewer
interpolation or regression models are built in the offline stage,
thus reducing the test time to evaluate the RB solution online.
Specifically, for each EM field component, 2d<sub>u</sub> &sum;<sub>&ell;=1</sub><sup>d<sub>u</sub></sup> Q<sub>u,&ell;</sub>
interpolation-based models are built in [20], where Q<sub>u,&ell;</sub> de-
notes the number of singular values truncated of coefficient
matrix P<sub>u</sub><sup>&ell;</sup>, while only 2R<sub>u</sub> interpolation-based models are
built in POD-CP ROM.

The recovery error the POD-CPD-CSI ROM at $(t^*, \boldsymbol{\mu}^*)$ will
be written as

$$
\begin{align}
\epsilon &= \|\mathbf{u}_h(t^*, \boldsymbol{\mu}^*) - \mathbf{u}_h^r(t^*, \boldsymbol{\mu}^*)\|^2 \nonumber \\
&= \|\mathbf{u}_h(t^*, \boldsymbol{\mu}^*) - \mathbf{V}_u \hat{\boldsymbol{\alpha}}_{\mathbf{u}}(t^*, \boldsymbol{\mu}^*)\|^2 \nonumber \\
&= \|\mathbf{u}_h(t^*, \boldsymbol{\mu}^*) - \mathbf{V}_u \mathbf{V}_u^T \mathbf{u}_h(t^*, \boldsymbol{\mu}^*) + \mathbf{V}_u \mathbf{V}_u^T \mathbf{u}_h(t^*, \boldsymbol{\mu}^*) \nonumber \\
&\qquad - \mathbf{V}_u \hat{\boldsymbol{\alpha}}_{\mathbf{u}}(t^*, \boldsymbol{\mu}^*)\|^2 \nonumber \\
&\le \epsilon_{\text{POD}}^2 + \|\mathbf{V}_u\| \|\mathbf{V}_u^T \mathbf{u}_h(t^*, \boldsymbol{\mu}^*) - \hat{\boldsymbol{\alpha}}_{\mathbf{u}}(t^*, \boldsymbol{\mu}^*)\|^2 \nonumber \\
&= \epsilon_{\text{POD}}^2 + \sum_{i=1}^{N} (\alpha_{\mathbf{u},i}(t^*, \boldsymbol{\mu}^*) - \hat{\alpha}_{\mathbf{u},i}(t^*, \boldsymbol{\mu}^*))^2 \nonumber \\
&= \epsilon_{\text{POD}}^2 + \epsilon_{\text{coe}}^2, \tag{10}
\end{align}
$$

where $\epsilon_{\text{POD}}^2$ is the truncation error of the two-step POD. And
$\epsilon_{\text{coe}}^2$ denotes the error of predicting the projection coefficients,
which includes the error of CP decomposition and the CSI
error of $\hat{\varphi}_{\mathbf{u},r}(t^*)$ and $\hat{\psi}_{\mathbf{u},r}(\boldsymbol{\mu}^*)$.

IV. NUMERICAL RESULTS

In this section, we present four numerical tests to evalu-
ate the performance of the proposed POD-CPD-CSI ROM.
Moreover we compare the proposed ROM with the POD-ANN
ROM [31] and POD-CSI ROM [20] both in time-efficiency
and errors. We consider non-magnetic materials ($\nu_r = 1$) in
a vacuum for all numerical tests. The norm of a tensor $\mathcal{G}$ is
defined as

$$
\|\mathcal{G}\| = \sqrt{\sum_{i}^{N_t} \sum_{j}^{N_p} \sum_{k}^{d_{\mathbf{u}}} [(\mathcal{G})_{ijk}]^2}. \quad (11)
$$
<figure><img src="image_2.png" /><figcaption>Fig. 2. Offline stage of the proposed POD-CPD-CSI ROM.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 3. Online stage of the proposed POD-CPD-CSI ROM.</figcaption></figure>

Algorithm 3: POD-CPD-CSI ROM

1 **Function** $[V_u, \hat{\omega}_u]=\text{Offline}(\mathcal{P}, \mathcal{T})$

2 Prepare a training parameter sampling
$\mathcal{P}_h^{\text{tr}} = \{\mu_1, \mu_2, \dots, \mu_{N_p}\} \subset \mathcal{P}$;

3 Compute the high-fidelity solutions $\mathbf{u}_h(t_i, \boldsymbol{\mu}_j)$ in the time domain $\mathcal{T}_h \subset \mathcal{T}$ for each physical parameter $\boldsymbol{\mu}_j$, $j = 1, 2, \dots, N_p$;

4 Prepare a training time sampling
$\mathcal{T}_h^{\text{tr}} = \{t_1, t_2, \dots, t_{N_t}\} \subseteq \mathcal{T}_h \subset \mathcal{T}$ and assemble a snapshot matrix $\mathbf{S}_u$;

5 Compute the POD basis matrix $\mathbf{V}_u$ via the two-step POD method in Algorithm 1;

6 Perform CP decomposition of the projection coefficient tensor $\mathcal{G}_u$ via the ALS method in Algorithm 2;

7 Train the CSI-based model $\hat{\omega}_u$ based on 7.

8 end

9 **Function** $[\mathbf{u}_h^r(t^*, \boldsymbol{\mu}^*)] = \text{Online}(\mathbf{V}_u, \hat{\omega}_u, (t^*, \boldsymbol{\mu}^*))$

10 Calculate the projection coefficient $\hat{\alpha}_u(t^*, \boldsymbol{\mu}^*)$ through CSI-based model $\hat{\omega}_u$;

11 Evaluate the RB solution $\mathbf{u}_h^r(t^*, \boldsymbol{\mu}^*)$ via formula 9.

12 end

We determine the truncation parameter $R_u$ of the CP decomposition by the relative error indicator $e_{u,CP} = \|\mathcal{G}_u - \hat{\mathcal{G}}_u\| / \|\mathcal{G}_u\|$. In the online phase, for a new time value $t \in \mathcal{T}_h^{te}$ and physical parameter value $\boldsymbol{\mu} \in \mathcal{T}_h^{te}$, the relative projection error

$$e_{u,\text{Pro}}(t, \boldsymbol{\mu}) = \frac{\|\mathbf{u}_h(t, \boldsymbol{\mu}) - \mathbf{V}_u \mathbf{V}_u^T \mathbf{u}_h(t, \boldsymbol{\mu})\|_2}{\|\mathbf{u}_h(t, \boldsymbol{\mu})\|_2}, \quad (12)$$

provides a lower bound for the relative error of the POD-CPD-CSI ROM

$$e_{u,\text{POD-CPD-CSI}}(t, \boldsymbol{\mu}) = \frac{\|\mathbf{u}_h(t, \boldsymbol{\mu}) - \mathbf{V}_u \hat{\alpha}_u(t, \boldsymbol{\mu})\|_2}{\|\mathbf{u}_h(t, \boldsymbol{\mu})\|_2}. \quad (13)$$

Furthermore, the average relative errors on a test set $\mathcal{T}_h^{te} \times \mathcal{P}_h^{te}$ of size $N_{te}$ are calculated to measure the accuracy of the POD-CPD-CSI ROM, which are defined as

$$\bar{e}_{u,\text{POD-CPD-CSI}} = \frac{\sum_{(t, \boldsymbol{\mu}) \in \mathcal{T}_h^{te} \times \mathcal{P}_h^{te}} e_{u,\text{POD-CPD-CSI}}(t, \boldsymbol{\mu})}{N_{te}}, \quad (14)$$

$$\bar{e}_{u,\text{Pro}} = \frac{\sum_{(t, \boldsymbol{\mu}) \in \mathcal{T}_h^{te} \times \mathcal{P}_h^{te}} e_{u,\text{Pro}}(t, \boldsymbol{\mu})}{N_{te}}. \quad (15)$$

We use the function `cp_als` in the MATLAB library `Tensor Toolbox [2]` to decompose the coefficient tensor $\mathcal{G}_u$. The CSI-based approximation of the projection coefficients is built via the MATLAB functions `spline` and `griddedInterpolant` for single-layer cases and multi-layer cases respectively. We refer to [13], [32] for more details of the multivariate CSI method. The architecture of ANN used in the POD-ANN ROM is listed in Table VI. All simulations are run on a desktop with an Intel Core i9 10-core 2.8 GHz × 20 CPU and 64 GB RAM.
<table><caption>TABLE I
CONFIGURATION OF THE NUMERICAL TESTS.</caption><thead><tr><th rowspan="2">Case</th><th colspan="4">Relative permittivity μ<sub>i</sub> of each layer</th><th colspan="4">Radius of each layer (m)</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>disk</td><td>[1,5]</td><td>-</td><td>-</td><td>-</td><td>0.6</td><td>-</td><td>-</td><td>-</td></tr><tr><td>multi-layer disk</td><td>[5.0, 5.6]</td><td>[3.25, 3.75]</td><td>[2.0, 2.5]</td><td>[1.25,1.75]</td><td>0.15</td><td>0.3</td><td>0.45</td><td>0.6</td></tr><tr><td>sphere</td><td>[2,6]</td><td>-</td><td>-</td><td>-</td><td>1.0</td><td>-</td><td>-</td><td>-</td></tr><tr><td>multi-layer sphere</td><td>[4.0,4.5]</td><td>[8.0, 8.5]</td><td>[2.0, 2.5]</td><td>-</td><td>0.6</td><td>1.0</td><td>1.5</td><td>-</td></tr></tbody></table>

<table><caption>TABLE II
PARAMETERS SETTING IN THE TWO-STEP POD METHOD.</caption><thead><tr><th>case</th><th>ρ<sub>t</sub></th><th>ρ<sub>μ</sub></th><th>d<sub>H<sub>x</sub></sub></th><th>d<sub>H<sub>y</sub></sub></th><th>d<sub>H<sub>z</sub></sub></th><th>d<sub>E<sub>x</sub></sub></th><th>d<sub>E<sub>y</sub></sub></th><th>d<sub>E<sub>z</sub></sub></th></tr></thead><tbody><tr><td>disk</td><td>0.1</td><td>1e-5</td><td>20</td><td>35</td><td>-</td><td>-</td><td>-</td><td>22</td></tr><tr><td>multi-layer disk</td><td>0.1</td><td>1e-5</td><td>17</td><td>15</td><td>-</td><td>-</td><td>-</td><td>15</td></tr><tr><td>sphere</td><td>0.1</td><td>1e-6</td><td>45</td><td>44</td><td>39</td><td>43</td><td>37</td><td>36</td></tr><tr><td>multi-layer sphere</td><td>0.1</td><td>1e-6</td><td>48</td><td>47</td><td>43</td><td>45</td><td>44</td><td>40</td></tr></tbody></table>

<table><caption>TABLE III
PARAMETERS SETTING IN THE CP DECOMPOSITION OF POD-CPD-CSI ROM.</caption><thead><tr><th>case</th><th>R<sub>H<sub>x</sub></sub></th><th>R<sub>H<sub>y</sub></sub></th><th>R<sub>H<sub>z</sub></sub></th><th>R<sub>E<sub>x</sub></sub></th><th>R<sub>E<sub>y</sub></sub></th><th>R<sub>E<sub>z</sub></sub></th></tr></thead><tbody><tr><td>disk</td><td>50</td><td>50</td><td>-</td><td>-</td><td>-</td><td>50</td></tr><tr><td>multi-layer disk</td><td>60</td><td>60</td><td>-</td><td>-</td><td>-</td><td>60</td></tr><tr><td>sphere</td><td>60</td><td>60</td><td>60</td><td>60</td><td>60</td><td>60</td></tr><tr><td>multi-layer sphere</td><td>50</td><td>50</td><td>50</td><td>50</td><td>50</td><td>50</td></tr></tbody></table>

<table><caption>TABLE IV
NUMBER OF CSI MODELS BUILT IN THREE ROMS.</caption><thead><tr><th>case</th><th>POD-ANN</th><th>POD-CSI</th><th>POD-CPD-CSI</th></tr></thead><tbody><tr><td>disk</td><td>-</td><td>1330</td><td>300</td></tr><tr><td>multi-layer disk</td><td>-</td><td>819</td><td>360</td></tr><tr><td>sphere</td><td>-</td><td>2196</td><td>720</td></tr><tr><td>multi-layer sphere</td><td>-</td><td>2423</td><td>600</td></tr></tbody></table>

<table><caption>TABLE V
SIZE OF TRAINING AND TEST DATA SETS.</caption><thead><tr><th rowspan="2">Case</th><th colspan="2">Training data set</th><th colspan="2">Test data set</th></tr><tr><th>N<sub>p</sub></th><th>N<sub>t</sub></th><th>N<sub>p</sub></th><th>N<sub>t</sub></th></tr></thead><tbody><tr><td>disk</td><td>81</td><td>263</td><td>4</td><td>263</td></tr><tr><td>multi-layer disk</td><td>3<sup>4</sup></td><td>253</td><td>4</td><td>253</td></tr><tr><td>sphere</td><td>81</td><td>100</td><td>4</td><td>100</td></tr><tr><td>multi-layer sphere</td><td>5<sup>3</sup></td><td>100</td><td>4</td><td>100</td></tr></tbody></table>

<table><caption>TABLE VI
ARCHITECTURE OF THE ANN USED FOR THE POD-ANN ROM.</caption><thead><tr><th rowspan="2">Case</th><th colspan="6">Layer</th></tr><tr><th>i</th><th>h<sub>1</sub></th><th>h<sub>2</sub></th><th>h<sub>3</sub></th><th>h<sub>4</sub></th><th>h<sub>5</sub></th><th>o</th></tr></thead><tbody><tr><td>disk</td><td>2</td><td>20</td><td>40</td><td>60</td><td>-</td><td>-</td><td>77</td></tr><tr><td>multi-layer disk</td><td>5</td><td>10</td><td>20</td><td>30</td><td>-</td><td>-</td><td>47</td></tr><tr><td>sphere</td><td>2</td><td>50</td><td>55</td><td>60</td><td>65</td><td>70</td><td>242</td></tr><tr><td>multi-layer sphere</td><td>4</td><td>50</td><td>55</td><td>60</td><td>65</td><td>70</td><td>267</td></tr></tbody></table>

A. Scattering of a plane wave by a dielectric disk

First, we consider the scattering problem described by the 2D time-domain TM waves in the normalized form

$$ \left\{ \begin{aligned} & \nu_r \frac{\partial \mathbf{H}}{\partial t} + \text{curl } E = 0, \\ & \varepsilon_r \frac{\partial E}{\partial t} - \text{curl } \mathbf{H} = 0, \end{aligned} \right. \qquad (16) $$

where the magnetic field $\mathbf{H} = (H_x, H_y)^T$ and the electric field $E = E_z$. An incident plane wave takes the form $H_x^{\text{inc}}(x, y, t) = 0$, $H_y^{\text{inc}}(x, y, t) = -\cos(\omega t - kx)$ and $E_z^{\text{inc}}(x, y, t) = \cos(\omega t - kx)$ with a frequency $f = 300$ MHz [34]. The computational domain is a square $\Omega = [-2.6\lambda_0, 2.6\lambda_0]^2$ and the radius of the inner disk is 0.6 m (background wavelength $\lambda_0 = 1$ m). The range of the relative permittivity $\varepsilon_r$ is presented in Table I, i.e., the parameter domain we are interested in is $\mathcal{P} = \{\mu = \varepsilon_r \in [1, 5]\}$. To determine the number of sampling parameter points in $\mathcal{P}$, we create three training parameter sets uniformly with different sampling intervals $\Delta p = 1, 0.5$ and $0.05$, resulting in $N_p = 5, 9$ and $81$ points respectively. A DGTD solver [19] is used to generate a collection of high-fidelity solutions of the three training parameter sets with $N_t = 263$ samples in the last period $T_h^{tr} = \{49.0024, 49.006, \dots, 49.966\}$. We evaluate the proposed method on a test parameter set $\mathcal{P}_h^{te} = \{1.215, 2.215, 3.215, 4.215\}$ and a test time set $T_h^{te} = T_h^{tr}$. Note that the truncation parameter $R_u$ is selected carefully such that the projection coefficient tensor $G_u$ lose little information in the CP decomposition. Fig. 4 shows the average relative errors $\bar{e}_{u,\text{Pro}}$ and $\bar{e}_{u,\text{POD-CPD-CSI}}$ on the training set $\mathcal{P}_h^{tr} \times \mathcal{T}_h^{tr}$ with a fixed truncation parameter $\rho_t = 0.1$. The performance of the POD-CP ROM mainly depends on the size of training set as well as the truncation parameter $\rho_\mu$ used. So the training parameter set $\mathcal{P}_h^{tr} = \{1, 1.05, \dots, 5\}$ is adopted in the following discussion.

The convergence curves of the average relative errors $\bar{e}_{\mathbf{H},\text{Pro}}$ and $\bar{e}_{\mathbf{E},\text{Pro}}$ are plotted in Fig. 5, from which we can see that the tolerance $\rho_t$ has little impact on the performance of the two-step POD method. By decreasing the tolerance $\rho_\mu$ from $1 \times 10^{-1}$ to $1 \times 10^{-5}$, there is a significant improvement of the two-step POD performance. Therefore, we set $\rho_t = 1 \times 10^{-1}$ and $\rho_\mu = 1 \times 10^{-5}$. The number of POD basis functions $d_u$ generated by the two-step POD method are listed in Table II. The relative errors $e_{u,\text{CP}}$ for different truncation parameters $R_u$ adopted in the CP decomposition are plotted in Fig. 6.
<table><caption>TABLE VII
AVERAGE RELATIVE ERRORS OF THREE ROMS.</caption><thead><tr><th rowspan="2">Case</th><th colspan="2">Projection</th><th colspan="2">POD-ANN</th><th colspan="2">POD-CSI</th><th colspan="2">POD-CPD-CSI</th></tr><tr><th>e&#773;<sub>H</sub></th><th>e&#773;<sub>E</sub></th><th>e&#773;<sub>H</sub></th><th>e&#773;<sub>E</sub></th><th>e&#773;<sub>H</sub></th><th>e&#773;<sub>E</sub></th><th>e&#773;<sub>H</sub></th><th>e&#773;<sub>E</sub></th></tr></thead><tbody><tr><td>disk</td><td>1.11e-2</td><td>1.36e-2</td><td>3.24e-2</td><td>4.27e-2</td><td>1.16e-2</td><td>1.36e-2</td><td>1.19e-2</td><td>1.36e-2</td></tr><tr><td>multi-layer disk</td><td>7.37e-3</td><td>8.23e-3</td><td>2.60e-2</td><td>3.03e-2</td><td>8.01e-3</td><td>8.26e-3</td><td>8.03e-3</td><td>8.86e-3</td></tr><tr><td>sphere</td><td>6.64e-4</td><td>1.44e-3</td><td>3.97e-2</td><td>3.21e-2</td><td>6.73e-4</td><td>1.49e-3</td><td>6.89e-4</td><td>1.50e-3</td></tr><tr><td>multi-layer sphere</td><td>7.34e-4</td><td>4.02e-3</td><td>3.10e-2</td><td>4.03e-2</td><td>7.50e-4</td><td>4.12e-3</td><td>7.58e-4</td><td>4.15e-3</td></tr></tbody></table>

<table><caption>TABLE VIII
COMPARISON OF THE CPU TIME OF FOM AND THREE ROMS. THE UNIT OF TIME IS SECOND.</caption><thead><tr><th rowspan="2">Case</th><th colspan="2">FOM</th><th colspan="2">POD-ANN</th><th colspan="2">POD-CSI</th><th colspan="2">POD-CPD-CSI</th></tr><tr><th>DoFs</th><th>CPU time</th><th>training</th><th>test</th><th>training</th><th>test</th><th>training</th><th>test</th></tr></thead><tbody><tr><td>disk</td><td>30264</td><td>258</td><td>32450</td><td>0.02</td><td>9</td><td>0.25</td><td>19</td><td>0.04</td></tr><tr><td>multi-layer disk</td><td>37236</td><td>305</td><td>31784</td><td>0.03</td><td>7</td><td>0.19</td><td>16</td><td>0.05</td></tr><tr><td>sphere</td><td>91125</td><td>343</td><td>43287</td><td>0.03</td><td>23</td><td>0.34</td><td>35</td><td>0.07</td></tr><tr><td>multi-layer sphere</td><td>91125</td><td>3710</td><td>39749</td><td>0.03</td><td>27</td><td>0.39</td><td>38</td><td>0.06</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 4. The average relative errors of $\bar{e}_{H,\text{Pro}}$ and $\bar{e}_{H,\text{POD-CP}}$ (top), $\bar{e}_{E,\text{Pro}}$ and $\bar{e}_{E,\text{POD-CP}}$ (bottom) on the training set with different sampling intervals $\Delta p$ and truncation parameters $\rho_{\mu}$.</figcaption></figure>

We take $R_{H_x} = R_{H_y} = R_{E_z} = 50$ to minimize the information lost in CP decomposition while constructing as few interpolation models as possible.

<figure><img src="image_5.png" /><figcaption>Fig. 5. The average relative errors of $\bar{e}_{H,\text{Pro}}$ (top) and $\bar{e}_{E,\text{Pro}}$ (bottom) with different truncation parameters $\rho_t$ and $\rho_{\mu}$.</figcaption></figure>

For the parameter instance $\mu_3 = 3.215$, comparison of the time evolution of the field $H_y$ and $E_z$ at a fixed point is displayed in Fig. 7. Moreover, the total field of $|H_x|$ calculated by DGTD and the absolute errors between POD-ANN, POD-CSI, POD-CPD-CSI and DGTD respectively are shown in Fig. 8. The time evolution of the relative error $e_{u,\text{Pro}}$ and $e_{u,\text{POD-ANN}}$, $e_{u,\text{POD-CSI}}$ and $e_{u,\text{POD-CPD-CSI}}$ is plotted in

<figure><img src="image_6.png" /><figcaption>Fig. 6. The relative errors $\bar{e}_{u,\text{CP}}$ with different truncation parameters $R_u$.</figcaption></figure>

Fig. 9. The average projection error $\bar{e}_{u,\text{Pro}}$ and three ROMs are listed in Table VII. It can be seen that the accuracy of POD-CSI and POD-CP-CSI are similar, and both are better than POD-ANN.

In addition, Table VIII shows the average computing time of a time trajectory matrix (1) between DGTD and three ROMs. The online test time of the ROM achieves 5520 times faster than that of the DGTD solver, which demonstrates the high efficiency of the proposed method. The non-intrusive ROMs features an offline-online framework, i.e., offline training is completely decoupled from online testing. Given a new parameter, the RB solution is rapidly recovered without an access to the high-dimensional FOM, which greatly improves the efficiency of the ROM. In addition, POD-CSI and POD-CPD-CSI take less training time offline compared with POD-ANN [31]. Finally, as shown in Table IV, POD-CPD-CSI generates
less CSI models and takes less test time compared to POD-CSI.

<figure><img src="image_7.png" /><figcaption>Fig. 7. Comparison of the time evolution of the field $H_y$ (top) and $E_z$ (bottom) at a fixed point.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 8. The DGTD solution of $|H_x|$ (a) at $t = 49.966$ m and the absolute errors between POD-ANN (b), POD-CSI (c), POD-CPD-CSI (d) and DGTD respectively.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 9. Comparison of the time evolution of the relative errors $e_{u,Pro}$, $e_{u,POD-ANN}$, $e_{u,POD-CSI}$ and $e_{u,POD-CPD-CSI}$ for the field $H$ and $E$ with the testing-parameter instances $\varepsilon_3 = 3.215$.</figcaption></figure>

## B. Scattering of a plane wave by a multi-layer disk

Next, we consider the case of a multi-layer dielectric disk as shown in Fig. 11. A square $\Omega = [-3.2\lambda_0, 3.2\lambda_0]^2$ is the computational domain, and the radius and the range of material parameters of each layer are presented in Table I. The excitation is the same as the case in IV-A.

To determine the number of sampling parameter points in $\mathcal{P} = \mathcal{P}^1 \times \mathcal{P}^2 \times \mathcal{P}^3 \times \mathcal{P}^4$ with $\mu = (\varepsilon_{r,1}, \varepsilon_{r,2}, \varepsilon_{r,3}, \varepsilon_{r,4})$ and $\varepsilon_{r,i} \in \mathcal{P}^i$, we create two training parameter sets uniformly

with $\mathcal{N}_{\Delta p} = 2$ and 3 in each layer, resulting in $\mathcal{N}_p = 16$ and 81 points in total. For each parameter point, $\mathcal{N}_t = 253$ snapshots are calculated in time at the last period ($\Delta t = 0.0038$ m and $\mathcal{T}_h^{tr} = \{49.0002, 49.0041, \dots, 49.9669\}$). We evaluate the proposed method on a test set $\mathcal{P}_h^{te} \times \mathcal{T}_h^{te}$ where $\mathcal{T}_h^{te} = \mathcal{T}_h^{tr}$ and $\mathcal{P}_h^{te} = \{\mu_1, \mu_2, \mu_3, \mu_4\}$ with $\mu_1 = \{(5.1, 3.3, 2.2, 1.3)\}$, $\mu_2 = \{(5.3, 3.5, 2.2, 1.5)\}$, $\mu_3 = \{(5.5, 3.4, 2.5, 1.7)\}$ and $\mu_4 = \{(5.6, 3.7, 2.4, 1.6)\}$.

<figure><img src="image_10.png" /><figcaption>Fig. 10. The average relative errors of $\bar{e}_{H,Pro}$ and $\bar{e}_{H,POD-CPD-CSI}$ (top), $\bar{e}_{E,Pro}$ and $\bar{e}_{E,POD-CPD-CSI}$ (bottom) on the training set with different number of sampling points $\mathcal{N}_{\Delta p}$ and truncation parameters $\rho_{\mu}$.</figcaption></figure>

The average relative errors $\bar{e}_{u,Pro}$ and $\bar{e}_{u,POD-CPD-CSI}$ on the training set $\mathcal{P}_h^{tr} \times \mathcal{T}_h^{tr}$ with a fixed truncation parameter $\rho_t = 0.1$ are shown in Fig. 10, therefore we select $\mathcal{N}_{\Delta p} = 3$ parameter points uniformly in each layer. The parameter setting in the POD-CPD-CSI ROM is listed in TableIII.

The average computing time of a time trajectory matrix and average relative errors are presented in TableVIII and TableVII respectively.

<figure><img src="image_11.png" /><figcaption>Fig. 11. Configuration and scattering of a plane wave by a multi-layer disk (top) and multi-layer sphere (bottom).</figcaption></figure>

## C. Scattering of a plane wave by a dielectric sphere

We consider a plane wave that is polarized in the $x$ direction, propagates in the $+z$ direction, and is scattered by a dielectric sphere [9]. The plane wave is expressed as

$$ \mathbf{E}^{\text{inc}} = e^{-jks} e_x, \quad \mathbf{H}^{\text{inc}} = \frac{E_x^{\text{inc}}}{\eta} e^{-jks} e_y, \qquad (17) $$

where the magnitude of the electric part is 1 V/m with a frequency $f = 300$ MHz. The computational domain is a cube $\Omega = [-0.2\lambda_0, 2.2\lambda_0]^3$ with $\lambda_0 = 1$ m, and the range of material parameters is listed in Table I. During the offline phase, high-fidelity solutions are calculated by the Mie series method [9] on $\mathcal{N}_h = 91125$ nodes at $\mathcal{N}_p = 81$ uniformly distributed parameter points ($\mathcal{P}_h^{tr} = \{2, 2.05, \dots, 5.95, 6\}$) with
$N_t = 100$ points in a period. We evaluate the proposed method on a test parameter set $P_h^{te} = \{2.215, 3.215, 4.215, 5.215\}$ and test time set $T_h^{te} = T_h^{tr} = \{0, 0.0101, \dots, 1\}$. The number of basis functions and truncation tolerances of the two-step POD and CP decomposition are presented in TableIII.

The average computing time and average relative errors are presented in TableVIII and TableVII respectively.

### D. Scattering of a plane wave by a multi-layer sphere

Finally, we consider a plane wave scattered by a 3-layer dielectric sphere [9] as shown in Fig. 11. The computational domain is $\Omega = [-0.2\lambda_0, 2.2\lambda_0]^3$, and the radius and range of material parameters of each layer are listed in Table I. The excitation is the same as the case in IV-C. In the offline phase, full-order solutions are calculated by the Mie series method on $N_h = 91125$ nodes at $N_p = 125$ uniformly distributed parameter points. A test parameter set $P_h^{te} = \{(2.1, 8.4, 4.2), (2.2, 8.3, 4.4), (2.3, 8.2, 4.1), (2.4, 8.1, 4.3)\}$ and test time set $T_h^{te} = T_h^{tr} = \{0, 0.0101, \dots, 1\}$ are used to evaluate the proposed method. The number of basis functions and truncation tolerances of two-step POD and CP decomposition are presented in TableIII.

The computing time and average relative errors are presented in TableVIII and TableVII respectively.

## V. CONCLUSION

A data-driven RB method based on POD, CP decomposition and CSI is proposed to accelerate the solution of parameterized time-domain Maxwell equations. RB solutions can be obtained at a very low cost from the online phase. We verify the performance of the proposed POD-CP ROM with scattering by a medium disk and a dielectric sphere. As for scattering problems with complex geometric objects, more POD basis functions generated by the two-step POD might be required to capture the main dynamics of the full-order system.

We remark that although data-driven ROMs are effective and easy to implement when sufficient training data (experimental or computational data) is available, it may be still expensive to generate relatively large data in many application fields. In addition, data-driven ROMs fit observed data well, but generally have poor generalization performance for unseen data. A possible approach to address these issues is to take physical information into account. In fact, the governing partial differential equations (PDEs) behind the training data are clearly known, but often ignored by data-driven ROMs. Future research direction is to combine the data and physical information to improve the generalization capacity of ROMs when training data is limited.

## ACKNOWLEDGEMENT

The last author is supported by NSFC, 12101511 and Fundamental Research Funds for the Central Universities (Grant No. JBK22YJ37).

## REFERENCES

[1] Francesco Ballarin, Andrea Manzoni, Alfio Quarteroni, and Gianluigi Rozza. Supremizer stabilization of pod-galerkin approximation of parametrized steady incompressible navier-stokes equations. *International Journal for Numerical Methods in Engineering*, 102(5):1136-1161, 2015.

[2] Brett W. Bader, Tamara G. Kolda and others. Tensor Toolbox for MATLAB. www.tensortoolbox.org. September 21, 2022.

[3] T Bui-Thanh, Murali Damodaran, and Karen Willcox. Proper orthogonal decomposition extensions for parametric applications in compressible aerodynamics. In *21st AIAA Applied Aerodynamics Conference*, page 4213, 2003.

[4] J Douglas Carroll and Jih-Jie Chang. Analysis of individual differences in multidimensional scaling via an n-way generalization of “eckart-young” decomposition. *Psychometrika*, 35(3):283–319, 1970.

[5] Fabien Casenave, Alexandre Ern, and Tony Lelievre. A nonintrusive reduced basis method applied to aeroacoustic simulations. *Advances in Computational Mathematics*, 41(5):961–986, 2015.

[6] Wenqian Chen, Qian Wang, Jan S Hesthaven, and Chuhua Zhang. Physics-informed machine learning for reduced-order modeling of nonlinear problems. *Journal of Computational Physics*, 446:110666, 2021.

[7] M Couplet, C Basdevant, and P Sagaut. Calibrated reduced-order pod-galerkin system for fluid flow modelling. *Journal of Computational Physics*, 207(1):192–220, 2005.

[8] Stefania Fresca, Luca Dede', and Andrea Manzoni. A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized pdes. *Journal of Scientific Computing*, 87:1–36, 2021.

[9] G. Kevin Zhu. Sphere scattering. https://www.mathworks.com/matlabcentral/fileexchange/31119-sphere-scattering, 2022.

[10] Francisco J Gonzalez and Maciej Balajewicz. Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems. *arXiv preprint arXiv:1808.01346*, 2018.

[11] Mengwu Guo and Jan S Hesthaven. Data-driven reduced order modeling for time-dependent problems. *Computer methods in applied mechanics and engineering*, 345:75–99, 2019.

[12] Richard A Harshman et al. Foundations of the parafac procedure: Models and conditions for an "explanatory" multimodal factor analysis. 1970.

[13] Md Sakib Hasan, Syed K Islam, and Benjamin J Blalock. Modeling of soi four-gate transistor (g4fet) using multidimensional spline interpolation method. *Microelectronics journal*, 76:33–42, 2018.

[14] Jan S Hesthaven, Gianluigi Rozza, Benjamin Stamm, et al. Certified reduced basis methods for parametrized partial differential equations, volume 590. Springer, 2016.

[15] Jan S Hesthaven and Stefano Ubbiali. Non-intrusive reduced order modeling of nonlinear problems using neural networks. *Journal of Computational Physics*, 363:55–78, 2018.

[16] Jian-Ming Jin. Theory and computation of electromagnetic fields. John Wiley & Sons, 2011.

[17] Tamara G Kolda and Brett W Bader. Tensor decompositions and applications. *SIAM review*, 51(3):455–500, 2009.

[18] Ettore Lappano, Frank Naets, Wim Desmet, Domenico Mundo, and Eugene Nijman. A greedy sampling approach for the projection basis construction in parametric model order reduction for structural dynamics models. In *Proceedings of the 27th international conference on noise and vibration engineering.*, pages 3563–3571. KU Leuven; Leuven, Belgium, 2016.

[19] Kun Li, Ting-Zhu Huang, Liang Li, and Stéphane Lanteri. Pod-based model order reduction with an adaptive snapshot selection for a discontinuous galerkin approximation of the time-domain maxwell's equations. *Journal of Computational Physics*, 396:106–128, 2019.

[20] Kun Li, Ting-Zhu Huang, Liang Li, and Stéphane Lanteri. Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation. *J. Sci. Comput.*, 87(2):1–29, 2021.

[21] Kun Li, Ting-Zhu Huang, Liang Li, Stéphane Lanteri, Li Xu, and Bin Li. A reduced-order discontinuous galerkin method based on pod for electromagnetic simulation. *IEEE Transactions on Antennas and Propagation*, 66(1):242–254, 2017.

[22] YC Liang, HP Lee, SP Lim, WZ Lin, KH Lee, and CG1237 Wu. Proper orthogonal decomposition and its applications—part i: Theory. *Journal of Sound and vibration*, 252(3):527–544, 2002.

[23] Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis. Deepxde: A deep learning library for solving differential equations. *SIAM Review*, 63(1):208–228, 2021.
[24] Takaaki Murata, Kai Fukami, and Koji Fukagata. Nonlinear mode decomposition with convolutional neural networks for fluid dynamics. *Journal of Fluid Mechanics*, 882:A13, 2020.

[25] VVS Prakash, , and Raj Mittra. Characteristic basis function method: A new technique for efficient solution of method of moments matrix equations. *Microwave and optical technology letters*, 36(2):95-100, 2003.

[26] Alfio Quarteroni, Gianluigi Rozza, and Andrea Manzoni. Certified reduced basis approximation for parametrized partial differential equations and applications. *Journal of Mathematics in Industry*, 1:1–49, 2011.

[27] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. *Journal of Computational physics*, 378:686–707, 2019.

[28] Sadasiva Rao, Donald Wilton, and Allen Glisson. Electromagnetic scattering by surfaces of arbitrary shape. *IEEE Transactions on antennas and propagation*, 30(3):409–418, 1982.

[29] Jiming Song, Cai-Cheng Lu, and Weng Cho Chew. Multilevel fast multipole algorithm for electromagnetic scattering by large complex objects. *IEEE transactions on antennas and propagation*, 45(10):1488–1493, 1997.

[30] Eric Suter and Juan R Mosig. A subdomain multilevel approach for the efficient mom analysis of large planar antennas. *Microwave and Optical Technology Letters*, 26(4):270–277, 2000.

[31] Qian Wang, Jan S Hesthaven, and Deep Ray. Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem. *Journal of computational physics*, 384:289–307, 2019.

[32] Ren-Hong Wang. Multivariate spline functions and their applications, volume 529. Springer Science & Business Media, 2013.

[33] D Xiao, F Fang, CC Pain, and IM Navon. A parameterized non-intrusive reduced order model and error analysis for general time-dependent nonlinear partial differential equations and its applications. *Computer Methods in Applied Mechanics and Engineering*, 317:868–889, 2017.

[34] Ying Zhao, Liang Li, and Kun Li. Non-intrusive reduced order modeling of parametric electromagnetic scattering problems through gaussian process regression. *arXiv preprint arXiv:2103.12472*, 2021.

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

**Stéphane Lanteri** is a senior research scientist at Inria center of Université Côte d'Azur. Sophia He holds a PhD in Engineering Sciences from University of Nice-Sophia Antipolis (Year of defense: 1991) and a HDR (habilitation to advise doctoral theses) from the University of Nice/Sophia Antipolis (Year of defense: 2003 - Title: High-performance numerical methods on unstructured meshes with applications to compressible fluid dynamics). From May 1992 to October 1993, he was a postdoctoral fellow at the Center for Aerospace Structure, University of

Colorado at Boulder, under the supervision of Charbel Farhat. His current research interests are concerned with numerical modeling of physical problems related to computational electromagnetics and computational nanophotonics: design of high order discontinuous Galerkin type approximation methods on unstructured meshes and of parallel numerical algorithms for solving differential systems modeling electromagnetic wave interaction with complex media with a focus of nanoscale light-matter interaction. Since February 2020, he is the scientific head of the Atlantis project-team at the Inria Research Center of Université Côte d'Azur. He is also the scientific coordinator of the development of the DIOGENES software suite, which is dedicated to computational nanophotonics. Stéphane Lanteri has authored or co-authored more than 80 publications in international journals in applied mathematics, scientific computing, computational physics and more recently, optics and photonics. Web: https://slanteri.github.io/

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

**Xiao-Feng He** received the B.S. degree in information and computational science from Sichuan Agricultural University, Yaan, China, in 2020, and the M.S. degree in computational mathematics from University of Electronic Science and Technology of China, Chengdu, China, in 2023. His research interests include data-driven reduced-order modeling and numerical solvers in computational electromagnetics.

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

**Kun Li** received the Ph.D. degree in computational mathematics from the University of Electronic Science and Technology of China, Chengdu, China, in 2020. He is currently an assistant professor with the School of Mathematics, Southwestern University of Finance and Economics, Chengdu, China. His current research interests include the discontinuous Galerkin time-domain methods and model order reduction techniques for solving the electromagnetic, and local/nonlocal optical response problems.

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

**Liang Li** received the B.S., M.Sc., and Ph.D. degrees in computational mathematics from the University of Electronic Science and Technology of China (UESTC), Chengdu, Sichuan, China, in 2004, 2007, and 2009, respectively. During 2010, he was a Visiting Scholar with the Université Libre de Bruxelles, Bruxelles, Belgium. From 2010 to 2012, he was a Postdoctoral Fellow with Inria Sophia Antipolis-Méditerranée research center, France. He is currently a full professor with the School of Mathematical Sciences, UESTC. His current research

interests include iterative methods and preconditioning techniques for solving large linear systems, hybridizable discontinuous Galerkin methods, domain decomposition, model order reduction techniques, (algebraic) multigrid methods, and with applications to electromagnetics and nanophotonics.