37M05, 65M60, 65R20, 68U07, 78A50

# Fast _A Posteriori_ State Error Estimation for Reliable Frequency Sweeping
in Microwave Circuits via the Reduced-Basis Method

Valentín de la Rubia<EMAIL_ADDRESS>Sridhar Chellappa
<EMAIL_ADDRESS><EMAIL_ADDRESS>benner@mpi-
magdeburg.mpg.de Lihong Feng  Peter Benner

###### Abstract

We develop a compact, reliable model order reduction approach for fast
frequency sweeps in microwave circuits by means of the reduced-basis method.
Contrary to what has been previously done, special emphasis is placed on
certifying the accuracy of the reduced-order model with respect to the
original full-order model in an effective and efficient way. Previous works on
model order reduction accuracy certification rely on costly _a posteriori_
error estimators, which typically require expensive _inf-sup_ constant
evaluations of the underlying full-order model. This scenario is often too
time-consuming and unaffordable in electromagnetic applications. As a result,
less expensive and heuristic error estimators are commonly used instead. Very
often, one is interested in knowing about the full state vector, instead of
just some output quantities derived from the full state. Therefore, error
estimators for the full state vector become relevant. In this work, we detail
the frequency behavior of both the electric field and the state error when an
approximation to the electric field solution is carried out. Both field
quantities share the same frequency behavior. Based on this observation, we
focus on the efficient estimation of the electric field state error and
propose a fast evaluation of the reduced-order model state error in the
frequency band of analysis, minimizing the number of full-order model
evaluations. This methodology is of paramount importance to carry out a
reliable fast frequency sweep in microwave circuits. Finally, real-life
applications will illustrate the capabilities and efficiency of the proposed
approach.

###### keywords:

Computer aided engineering, design automation, error analysis, finite element
method, Galerkin method, microwave circuits, reduced basis method, reduced
order modeling.

We propose a fast method for _a posteriori_ state error estimation for
reduced-order modeling of microwave circuits. Our approach leverages a
similarity in the frequency domain behavior of the electric field and the
state error, which enables us to minimize the number of full-order model
evaluations during error estimation. The proposed fast state error estimator
is used in a greedy algorithm to adaptively build the reduced-order model. The
efficiency and reliability of our approach are demonstrated via its
application to four challenging models of microwave circuits.

##  1 Introduction

Microwave engineering relies on time-consuming electromagnetic simulations to
carry out robust electrical designs. The electromagnetic complexity in
microwave devices is such that _only_ detailed, full-wave simulations, i.e.,
solving Maxwell’s equations directly, can guide an engineer in pursuing the
target electrical design. As a result, most of a microwave engineer’s working
time is spent on waiting for an electromagnetic simulation to conclude, which
will assist him in taking an action to meet the specifications in an
electrical design. Current industrial needs have long since recognized that
one can no longer afford this design methodology. Different efforts in
computational electromagnetics (CEM) community have been carried out to speed
up this costly process, and most of them follow the model order reduction
(MOR) philosophy [25, 31, 21, 34, 6, 8, 43, 1, 40].

A reduced-order model (ROM) implies replacing a rather complex physical model
by a much simpler mathematical one that still maintains certain physical
aspects of the original model over a parameter domain. The computational
complexity of the ROM should be insignificant in comparison to the high
computational cost of the original full-order model (FOM). MOR has
demonstrated its robustness in reducing the complexity of parametric systems
[32, 20, 4]. However, the accuracy of the ROM is sometimes not guaranteed due
to lack of low-cost and computable error estimators. Although the ROM may be
valid for a certain parameter range, its validity over the entire parameter
domain is not guaranteed. As a result, it can not be used as a reliable
surrogate of the original FOM. This lack of providing accuracy guarantees
precludes the ROM from being used for industrial applications, where no _a
priori_ knowledge of the parameter range may be available. This is the worst
case scenario for MOR. It is quite usual that the proposition of corresponding
error estimation lags behind new MOR algorithms. To remedy this, a great
effort has been carried out in certifying the accuracy of the ROM, where
computationally expensive error estimation may be allowed. This is the case
for _inf-sup_ constant-based error estimators [18, 19, 37, 35]. The residual
norm divided by this costly _inf-sup_ constant [14] bounds the state error. As
already stated in [12], keeping the _inf-sup_ constant in the denominator of
the error estimator causes potential risk for many problems with small _inf-
sup_ constants. This is quite common in microwave circuits, where resonances
show up in the frequency band of analysis, dropping the _inf-sup_ constant
down to zero [10, 14]. Different strategies for _a posteriori_ error
estimation should be considered, reducing its computational cost to the same
order of the ROM, if possible. Residual norm-based error estimation can be
carried out without effort and this has been often used as a heuristic error
estimator [10, 11, 42, 9, 33, 38, 2, 13, 28]. Going back to the state error
estimation, recent works have focused on avoiding the _inf-sup_ constant
evaluation [36, 5]. There, additional dual or residual systems are solved to
obtain the error estimators, overcoming any time-consuming _inf-sup_ constant
calculation. However, despite the fact that they avoid computing the expensive
_inf-sup_ constant, both approaches need to solve additional dual or residual
systems, respectively.

In this work, we aim to further reduce the computational costs of the error
estimator proposed in [5], which was shown to be more efficient than that in
[36]. We study the frequency behavior of both the electric field and the state
error when an approximation to the electric field solution is carried out, and
detail a Fourier series representation in both cases. Both field quantities
share the same orthogonal Fourier series representation basis in frequency-
parameter systems. In comparison to what has been previously done for general
parametric systems, where a residual system needs to be solved independently,
we focus on the efficient determination of the electric field state error and
propose a fast evaluation of the ROM state error in the frequency band of
analysis, minimizing the number of FOM evaluations which plays a central role
in determining the efficiency in MOR. This methodology is of paramount
importance to carry out a reliable fast frequency sweep in microwave circuits.

This paper is organized as follows. In Section 2 we review the time-harmonic
Maxwell’s equations in variational form, solve for the electromagnetic field
in order to show its frequency behavior and detail the _inf-sup_ constant-
based standard error analysis. Section 3 deals with the proposed state error
estimation in frequency-parameter systems avoiding the _inf-sup_ constant.
Numerical simulations in Section 4 show the performance of the proposed
approach for reliable fast frequency sweeps in electromagnetics. Real-life
microwave circuits illustrate the capabilities and accuracy of the proposed
methodology. Finally, in Section 5, we provide conclusions.

##  2 Problem Statement

The electromagnetic phenomena in a given device are described by Maxwell’s
equations. Applying the Fourier transform to these, the fields in the
transform domain $i\omega$ can be found. They are

$\displaystyle\nabla\times\mathbf{E}$
$\displaystyle=-i\omega\mu\mathbf{H}\text{ in }\Omega,$ (1a)
$\displaystyle\nabla\times\mathbf{H}$
$\displaystyle=i\omega\varepsilon\mathbf{E}\text{ in }\Omega,$ (1b)
$\displaystyle\mathbf{n}\times\mathbf{E}$ $\displaystyle=\mathbf{0}\text{ on
}\Gamma_{\text{PEC}},$ (1c) $\displaystyle\mathbf{n}\times\mathbf{H}$
$\displaystyle=\mathbf{0}\text{ on }\Gamma_{\text{PMC}},$ (1d)
$\displaystyle\mathbf{n}\times\mathbf{H}$ $\displaystyle=\mathbf{J}\text{ on
}\Gamma\text{,}$ (1e)

where $\Omega\subset\mathbb{R}^{3}$ is a source-free, sufficiently smooth
bounded domain, $\mathbf{n}$ is the unit outward normal vector on the boundary
$\partial\Omega$ of $\Omega$. The boundary is divided into perfect electric
conductor (PEC), perfect magnetic conductor (PMC) and ports, i.e.,
$\partial\Omega=\Gamma_{\text{PEC}}\cup\Gamma_{\text{PMC}}\cup\Gamma$.
$\mathbf{E}$ and $\mathbf{H}$ are the electric and magnetic fields,
$\varepsilon$ and $\mu$ are, respectively, the permittivity and permeability
of the medium, which is assumed to be lossless, and the tangential field
$\mathbf{J}$ is the excitation current at the ports. Time-harmonic Maxwell’s
equations can be written in a classical weak formulation over an appropriate
admissible function space $\mathcal{H}$, viz.

find $\displaystyle\mathbf{E}\in\mathcal{H}\leavevmode\nobreak\ \text{such
that}$ (2) $\displaystyle
a(\mathbf{E},\mathbf{v})=f(\mathbf{v})\leavevmode\nobreak\
\forall\mathbf{v}\in\mathcal{H}\text{.}$

The bilinear form is

$a(\mathbf{E},\mathbf{v})=\int\limits_{\Omega}\left(\frac{1}{\mu}\nabla\times\mathbf{E}\cdot\nabla\times\mathbf{v}-\omega^{2}\varepsilon\mathbf{E}\cdot\mathbf{v}\right)dx\text{,}$
(3)

and the linear form

$f(\mathbf{v})=i\omega\int\limits_{\partial\Omega}\mathbf{J}\cdot\mathbf{v}\leavevmode\nobreak\
ds=i\omega\int\limits_{\Gamma}\mathbf{J}\cdot\mathbf{v}\leavevmode\nobreak\
ds\text{.}$ (4)

Here, the admissible space $\mathcal{H}$ is a subspace of the Hilbert space
$H(curl,\Omega)$ defined by:

$H(curl,\Omega)=\left\\{\mathbf{u}\in
L^{2}(\Omega,\mathbb{C}^{3})\leavevmode\nobreak\ |\leavevmode\nobreak\
\nabla\times\mathbf{u}\in L^{2}(\Omega,\mathbb{C}^{3})\right\\}\text{,}$ (5)

since $\mathcal{H}$ should take the boundary condition (1c) into account,
namely,

$\mathcal{H}=\left\\{\mathbf{u}\in H(curl,\Omega)\leavevmode\nobreak\
|\leavevmode\nobreak\
\mathbf{n}\times\mathbf{u}=\mathbf{0}\leavevmode\nobreak\ \text{on
}\Gamma_{\text{PEC}}\right\\}\text{.}$ (6)

Let us refer to the trace spaces, namely,

$\displaystyle H^{-\frac{1}{2}}(div,\partial\Omega)$
$\displaystyle=\\{\mathbf{n}\times\mathbf{u}\leavevmode\nobreak\
\text{on}\leavevmode\nobreak\ \partial\Omega\leavevmode\nobreak\
|\leavevmode\nobreak\ \mathbf{u}\in H(curl,\Omega)\\}$ (7) $\displaystyle
H^{-\frac{1}{2}}(curl,\partial\Omega)$
$\displaystyle=\\{\mathbf{n}\times\mathbf{u}\times\mathbf{n}\leavevmode\nobreak\
\text{on}\leavevmode\nobreak\ \partial\Omega\leavevmode\nobreak\
|\leavevmode\nobreak\ \mathbf{u}\in H(curl,\Omega)\\}\text{,}$

and point out that they are dual to each other with the following duality
pairing

$\langle\mathbf{u},\mathbf{v}\rangle=\int\limits_{\partial\Omega}\mathbf{u}\cdot\mathbf{v}\leavevmode\nobreak\
ds$ (8)

$\mathbf{u}\in H^{-\frac{1}{2}}(div,\partial\Omega)$ and $\mathbf{v}\in
H^{-\frac{1}{2}}(curl,\partial\Omega)$. It is now apparent that the excitation
current $\mathbf{J}$ belongs to $H^{-\frac{1}{2}}(div,\partial\Omega)$. We
refer to [16, 29] for a through explanation for all these spaces.

###  2.1 Field Frequency Dependency in Electromagnetics

Following [10, 24], where some frequency structure is shown in the solution to
the variational problem (2), we introduce the Helmholtz decomposition

$\mathcal{H}=\mathcal{H}(curl0,\Omega)\oplus\mathcal{V}\text{,}$ (9)

where

$\displaystyle\mathcal{H}(curl0,\Omega)=\\{\mathbf{u}\in\mathcal{H}\leavevmode\nobreak\
|\leavevmode\nobreak\ \nabla\times\mathbf{u}=\mathbf{0}\\},$ (10a)
$\displaystyle\mathcal{V}=\\{\mathbf{u}\in\mathcal{H}\leavevmode\nobreak\
|\leavevmode\nobreak\
(\varepsilon\mathbf{u},\mathbf{v})_{L^{2}(\Omega)}=0\leavevmode\nobreak\
\forall\mathbf{v}\in\mathcal{H}(curl0,\Omega)\\}\text{.}$ (10b)

$(\cdot,\cdot)_{L^{2}(\Omega)}$ is the inner product in
$L^{2}(\Omega,\mathbb{C}^{3})$. $\mathcal{H}(curl0,\Omega)$ denotes the
nullspace of the curl operator while $\mathcal{V}$ stands for its orthogonal
complement within the solution space $\mathcal{H}$ in the following inner
product

$(\mathbf{u},\mathbf{v})_{\mu,\varepsilon}=(\frac{1}{\mu}\nabla\times\mathbf{u},\nabla\times\mathbf{v})_{L^{2}(\Omega)}+(\varepsilon\mathbf{u},\mathbf{v})_{L^{2}(\Omega)}\text{.}$
(11)

It should be noted that both $\mathcal{H}(curl0,\Omega)$ and $\mathcal{V}$
spaces satisfy the PEC boundary condition on $\Gamma_{\text{PEC}}$.

The variational problem (2) can be solved by using the splitting
$\mathbf{E}=\mathbf{E}_{0}+\mathbf{e}$, $\mathbf{E}_{0}\in\leavevmode\nobreak\
\mathcal{H}(curl0,\Omega)$, $\mathbf{e}\in\mathcal{V}$. We refer to [10, 24]
for details. As a result, we can make the dependence of the solution to time-
harmonic Maxwell’s equations on frequency explicit, cf. [26, 7], i.e.,

$\displaystyle\text{if
}\omega^{2}\neq\omega_{n}^{2}\text{,}\leavevmode\nobreak\
\mathbf{E}=\mathbf{E}_{0}+\mathbf{e}=\frac{1}{i\omega}\mathbf{F}_{0}+i\omega\sum\limits_{n=1}^{\infty}\frac{A_{n}}{1-\frac{\omega^{2}}{\omega_{n}^{2}}}\mathbf{e}_{n},$
(12) $\displaystyle\text{if
}\omega^{2}=\omega_{n}^{2}\text{,}\leavevmode\nobreak\
\mathbf{E}=\mathbf{E}_{0}+\mathbf{e}=$
$\displaystyle=\frac{1}{i\omega}\mathbf{F}_{0}+i\omega\sum\limits_{\omega_{n}^{2}=\omega^{2}}a_{n}\mathbf{e}_{n}+i\omega\sum\limits_{\omega_{n}^{2}\neq\omega^{2}}\frac{A_{n}}{1-\frac{\omega^{2}}{\omega_{n}^{2}}}\mathbf{e}_{n}\text{.}$

$\mathbf{F}_{0}\in\mathcal{H}(curl0,\Omega)$ is related to the Riesz
representative for the electric field in statics. The set of eigenmodes
$\\{\mathbf{e}_{n}\leavevmode\nobreak\ |\leavevmode\nobreak\
n\in\leavevmode\nobreak\ \mathbb{N}\\}\subset\mathcal{V}$ stands for the
resonant modes in electrodynamics, along with their corresponding resonant
frequencies $\omega_{n}\in\mathbb{R}$, and forms a complete orthonormal system
in $\mathcal{V}$ with respect to the inner product (11) [24]. It should be
pointed out that $\mathcal{H}(curl0,\Omega)$ is orthogonal to $\mathcal{V}$
with respect to the same inner product (11). Getting to our point, (12)
details an orthogonal representation, i.e., a Fourier series for the electric
field where the frequency dependence is explicit. Further, $A_{n}$ are
coupling coefficients for the excitation current $\mathbf{J}$ to its
corresponding resonant mode $\mathbf{e}_{n}$ and are determined by

$A_{n}=\langle\mathbf{J},\mathbf{n}\times\overline{\mathbf{e}}_{n}\times\mathbf{n}\rangle\text{.}$
(13)

$\overline{\mathbf{e}}_{n}$ stands for the complex conjugate of
$\mathbf{e}_{n}$. Finally, $a_{n}$ are arbitrary coefficients since the
electric field is not unique at resonance.

###  2.2 Parametric Variational Problem and Standard _A Posteriori_ Error
Analysis

Taking frequency ($\omega$) as a parameter, the weak formulation for time-
harmonic Maxwell’s equations (2) turns into the following parametric
variational problem:

find $\displaystyle\mathbf{E}(\omega)\in\mathcal{H}\leavevmode\nobreak\
\text{such that}$ (14) $\displaystyle
a(\mathbf{E}(\omega),\mathbf{v};\omega)=f(\mathbf{v};\omega)\leavevmode\nobreak\
\forall\mathbf{v}\in\mathcal{H},\leavevmode\nobreak\
\forall\omega\in\mathcal{B}\text{.}$

where
$\mathcal{B}:=[\omega_{\text{min}},\omega_{\text{max}}]\subset\mathbb{R}$ is
the frequency band of interest, and the frequency-parameter bilinear and
linear forms $a(\cdot,\cdot;\omega)$ and $f(\cdot;\omega)$ are already defined
in (3) and (4), respectively. The well-posedness of the parametric problem
(14) relies on the behavior of the so-called _inf-sup_ constant
$\beta(\omega)$ as a function of frequency:

$\beta(\omega)=\adjustlimits{\inf}_{\mathbf{u}\in\mathcal{H}}{\sup}_{\mathbf{v}\in\mathcal{H}}\,\frac{|a(\mathbf{u},\mathbf{v};\omega)|}{\|\mathbf{u}\|_{\mathcal{H}}\|\mathbf{v}\|_{\mathcal{H}}}.$
(15)

For all $\omega\in\mathcal{B}$, $\beta(\omega)\geq\beta_{0}>0$ ensures the
well-posedness and uniqueness in the variational problem (14) [19].

This result gives rise to a standard _a posteriori_ error analysis. Provided
an approximate solution $\mathbf{\tilde{E}}(\omega)\in\mathcal{H}$ to the
variational problem (14) is found, the error in the approximation
$\|\mathbf{E}(\omega)-\leavevmode\nobreak\
\mathbf{\tilde{E}}(\omega)\|_{\mathcal{H}}$ can be bounded using the _inf-sup_
constant. Indeed, (15) can be rewritten as follows:

$\sup_{\mathbf{v}\in\mathcal{H}}\frac{|a(\mathbf{u},\mathbf{v};\omega)|}{\|\mathbf{u}\|_{\mathcal{H}}\|\mathbf{v}\|_{\mathcal{H}}}>\beta(\omega),\forall\mathbf{u}\in\mathcal{H}.$
(16)

In particular, this inequality still holds when replacing
$\mathbf{u}\in\mathcal{H}$ by the field
$\mathbf{E}(\omega)-\mathbf{\tilde{E}}(\omega)\in\mathcal{H}$, which gives
rise to an upper bound for the approximation error, namely,

$\|\mathbf{E}(\omega)-\mathbf{\tilde{E}}(\omega)\|_{\mathcal{H}}<\frac{1}{\beta(\omega)}\sup_{\mathbf{v}\in\mathcal{H}}\frac{|a(\mathbf{E}(\omega)-\mathbf{\tilde{E}}(\omega),\mathbf{v};\omega)|}{\|\mathbf{v}\|_{\mathcal{H}}}.$
(17)

However, this error bound not only involves the computation of the norm of the
residual functional

$\displaystyle r(\mathbf{\tilde{E}}(\omega),\mathbf{v};\omega):=$
$\displaystyle
f(\mathbf{v};\omega)-a(\mathbf{\tilde{E}}(\omega),\mathbf{v};\omega)$ (18)
$\displaystyle=$ $\displaystyle
a(\mathbf{E}(\omega)-\mathbf{\tilde{E}}(\omega),\mathbf{v};\omega),\forall\mathbf{v}\in\mathcal{H}\text{,}$

which can be determined in an efficient way as a function of frequency [11,
42, 18, 38, 33], but also the determination of the _inf-sup_ constant
throughout the frequency band of interest $\mathcal{B}$, which can be time-
consuming [19, 37, 14].

Furthermore, in microwave engineering, resonances appearing in $\mathcal{B}$
are responsible for the target electrical response. As a result, resonant
modes arise and the uniqueness of the solution is no longer valid in the band
of interest $\mathcal{B}$. The _inf-sup_ constant vanishes at the resonance
frequencies, giving rise to a near-infinity upper bound for the error in (17)
nearby resonances. In addition, the above error estimation leads to
unacceptable overestimation of the error even for well-conditioned problems
[35]. Given no better choices, the norm of the residual (18), which can be
straightforwardly computed, has been used as a _heuristic_ error estimator
[34, 9, 33, 2, 13, 28, 23]

##  3 State Error Estimation Avoiding the _Inf-Sup_ -Constant

The previous section has shown the main role the _inf-sup_ constant plays in
_a posteriori_ error estimation, as well as the incapability of _inf-sup_
constant-based error estimators to provide a tight error bound nearby resonant
frequencies. Unfortunately, the norm of the residual cannot provide a sharp
error estimation at or nearby resonance frequencies either. We shall elaborate
on this later in Section 5. As a result, we are in need of more efficient
state error estimators to certify the accuracy of the approximate field
solution to the frequency-parameter variational problem (14), even in the
presence of resonances. Recent efforts have moved towards this goal [36, 5].
There, instead of computing the _inf-sup_ constant, additional dual or
residual systems need to be solved to obtain the state error estimator. These
additional systems constitute an extra computational effort to certify the
accuracy of the approximate field solution. In this work, we focus on fast _a
posteriori_ state error estimator computation taking advantage of the
frequency dependency in the field solution (12) for the frequency-parameter
problem (14).

###  3.1 Field Error Frequency Dependency in Electromagnetics

Given an approximate solution $\mathbf{\tilde{E}}(\omega)\in\mathcal{H}$ to
(14), we can study the Fourier series representation of the state error

$\boldsymbol{\epsilon}(\omega):=\mathbf{E}(\omega)-\mathbf{\tilde{E}}(\omega)\in\mathcal{H}\text{,}$
(19)

making its frequency dependency explicit. The state error (19) satisfies the
frequency-parameter variational problem,

find
$\displaystyle\boldsymbol{\epsilon}(\omega)\in\mathcal{H}\leavevmode\nobreak\
\text{such that}$ (20) $\displaystyle
a(\boldsymbol{\epsilon}(\omega),\mathbf{v};\omega)=f^{\boldsymbol{\epsilon}}(\mathbf{v};\omega)\leavevmode\nobreak\
\forall\mathbf{v}\in\mathcal{H},\leavevmode\nobreak\
\forall\omega\in\mathcal{B}\text{.}$

The frequency-parameter bilinear form $a(\cdot,\cdot;\omega)$ is already
defined in (3). The frequency-parameter linear form
$f^{\boldsymbol{\epsilon}}(\cdot;\omega)$ is the residual functional
$r(\mathbf{\tilde{E}}(\omega),\cdot;\omega)$ detailed in (18),

$f^{\boldsymbol{\epsilon}}(\mathbf{v})=f(\mathbf{v};\omega)-a(\mathbf{\tilde{E}}(\omega),\mathbf{v};\omega)=i\omega\int\limits_{\partial\Omega}\mathbf{J}^{\boldsymbol{\epsilon}}\cdot\mathbf{v}\leavevmode\nobreak\
ds\text{,}$ (21)

which can be identified as a residual error current
$\mathbf{J}^{\boldsymbol{\epsilon}}\in\leavevmode\nobreak\
H^{-\frac{1}{2}}(div,\partial\Omega)$. By an analogous reasoning as the one in
Section 2, we get

$\displaystyle\text{if }\omega^{2}$
$\displaystyle\neq\omega_{n}^{2}\text{,}\leavevmode\nobreak\
\boldsymbol{\epsilon}(\omega)=\frac{1}{i\omega}\mathbf{F}^{\boldsymbol{\epsilon}}_{0}+i\omega\sum\limits_{n=1}^{\infty}\frac{A^{\boldsymbol{\epsilon}}_{n}}{1-\frac{\omega^{2}}{\omega_{n}^{2}}}\mathbf{e}_{n},$
(22) $\displaystyle\text{if }\omega^{2}$
$\displaystyle=\omega_{n}^{2}\text{,}$
$\displaystyle\boldsymbol{\epsilon}(\omega)$
$\displaystyle=\frac{1}{i\omega}\mathbf{F}^{\boldsymbol{\epsilon}}_{0}+i\omega\sum\limits_{\omega_{n}^{2}=\omega^{2}}a^{\boldsymbol{\epsilon}}_{n}\mathbf{e}_{n}+i\omega\sum\limits_{\omega_{n}^{2}\neq\omega^{2}}\frac{A^{\boldsymbol{\epsilon}}_{n}}{1-\frac{\omega^{2}}{\omega_{n}^{2}}}\mathbf{e}_{n}\text{.}$

$\mathbf{F}^{\boldsymbol{\epsilon}}_{0}\in\mathcal{H}(curl0,\Omega)$ is
related to the Riesz representative for the stationary error field. The same
set of eigenmodes $\\{\mathbf{e}_{n}\leavevmode\nobreak\ |\leavevmode\nobreak\
n\in\leavevmode\nobreak\ \mathbb{N}\\}\subset\mathcal{V}$ along with their
corresponding resonant frequencies $\omega_{n}$ as in (12) can be used.
$A^{\boldsymbol{\epsilon}}_{n}$ are coupling coefficients for the residual
error current $\mathbf{J}^{\boldsymbol{\epsilon}}$ to the corresponding
resonant mode $\mathbf{e}_{n}$, namely,

$A^{\boldsymbol{\epsilon}}_{n}=\langle\mathbf{J}^{\boldsymbol{\epsilon}},\mathbf{n}\times\overline{\mathbf{e}}_{n}\times\mathbf{n}\rangle\text{.}$
(23)

In addition, $a^{\boldsymbol{\epsilon}}_{n}$ are arbitrary coefficients since
there is no unique solution at resonance.

Having a closer look at equations (12) and (22), we can realize that the
solutions to both original and residual variational problems share the same
frequency-parameter behaviour and admit a similar Fourier series
representation with the same frequency pattern. We may be then tempted to use
the same representation basis to find an approximate solution to both original
and residual variational problems (14) and (20). However, this should be done
carefully to avoid stating that the state error
$\boldsymbol{\epsilon}(\omega)$ is identically zero even in the situation
where the approximate field $\mathbf{\tilde{E}}(\omega)$ may still be far away
from the true solution $\mathbf{E}(\omega)$. We will get back to this point
later in Section 3.2. In order to avoid this embarrassing scenario, [5]
proposes to approximate these variational problems by applying different
Galerkin projection spaces
$\mathcal{H}_{m},\mathcal{H}^{\boldsymbol{\epsilon}}_{m}:=\mathcal{H}_{m}+\mathcal{H}^{\mathbf{r}}_{m}\subset\mathcal{H}$,
to the original and residual problems, respectively, giving rise to
corresponding reduced systems which can be solved with ease, namely,

find
$\displaystyle\mathbf{\tilde{E}}(\omega)\in\mathcal{H}_{m}\leavevmode\nobreak\
\text{such that}$ (24a) $\displaystyle
a(\mathbf{\tilde{E}}(\omega),\mathbf{v};\omega)=f(\mathbf{v};\omega)\leavevmode\nobreak\
\forall\mathbf{v}\in\mathcal{H}_{m},\leavevmode\nobreak\
\forall\omega\in\mathcal{B}\text{,}$ and find
$\displaystyle\boldsymbol{\tilde{\epsilon}}(\omega)\in\mathcal{H}^{\boldsymbol{\epsilon}}_{m}\leavevmode\nobreak\
\text{such that}$ (24b) $\displaystyle
a(\boldsymbol{\tilde{\epsilon}}(\omega),\mathbf{v};\omega)=f^{\boldsymbol{\epsilon}}(\mathbf{v};\omega)\leavevmode\nobreak\
\forall\mathbf{v}\in\mathcal{H}^{\boldsymbol{\epsilon}}_{m},\leavevmode\nobreak\
\forall\omega\in\mathcal{B}\text{.}$

Algorithm 1 adaptively builds the reduced-basis spaces up, i.e., the Galerkin
projection spaces, in a greedy framework. As the number of iterations in this
procedure increases, $\boldsymbol{\tilde{\epsilon}}(\omega)$ approximates the
true state error $\boldsymbol{\epsilon}(\omega)$ better and better,
$\boldsymbol{\tilde{\epsilon}}(\omega)\approx\boldsymbol{\epsilon}(\omega)$,
and $\|\boldsymbol{\tilde{\epsilon}}(\omega)\|$ does perform as a sharp _a
posteriori_ state error estimator. We refer to [5] for the details. However,
it should be pointed out that, in Algorithm 1, the dimension of the Galerkin
projection space for the residual problem
$\mathcal{H}^{\boldsymbol{\epsilon}}_{m}$ may double the dimension in the
Galerkin projection space for the original problem $\mathcal{H}_{m}$ at each
iteration. Further, distinct sets of parameters are arranged to solve for both
the original and residual problems (14) and (20), which need to be done
independently, in spite of the same dynamics being observed in both
variational problems (see (12) and (22)). This may turn this procedure rather
time consuming. In this work, we focus on keeping a low computational effort,
taking advantage of the observations in Sections 2 and 3 to carry out a
reliable fast frequency sweep analysis.

Algorithm 1 Adaptive construction of the Galerkin projection spaces
$\mathcal{H}_{m},\mathcal{H}^{\boldsymbol{\epsilon}}_{m}$ for state error
estimation [5].

0: Frequency band of interest
$\mathcal{B}:=[\omega_{\text{min}},\omega_{\text{max}}]$, tolerance
$\texttt{tol}>0$ as the acceptable state error.

0: $\mathcal{H}_{m}$ to ensure tol state error in (24).

1: Initialize $\mathcal{H}_{m}=\\{\mathbf{0}\\}$,
$\mathcal{H}^{\mathbf{r}}_{m}=\\{\mathbf{0}\\}$, $\xi=\texttt{tol}+1$. Choose
different samples $\omega^{*}$ and $\omega^{*}_{\boldsymbol{\epsilon}}$
randomly taken from $\mathcal{B}$.

2: while $\xi>\texttt{tol}$ do

3: Solve for $\mathbf{E}(\omega^{*})$ in (14), orthonormalize and enrich
$\mathcal{H}_{m}$:
$\mathcal{H}_{m}=\mathcal{H}_{m}+\text{span}\\{\mathbf{E}(\omega^{*})\\}$.

4: Using $\mathcal{H}_{m}$, solve for the approximate field
$\mathbf{\tilde{E}}(\omega)$ in (24a).

5: Solve for $\boldsymbol{\epsilon}(\omega^{*}_{\boldsymbol{\epsilon}})$ in
(20), orthonormalize and enrich $\mathcal{H}^{\mathbf{r}}_{m}$:
$\mathcal{H}^{\mathbf{r}}_{m}=\mathcal{H}^{\mathbf{r}}_{m}+\text{span}\\{\boldsymbol{\epsilon}(\omega^{*}_{\boldsymbol{\epsilon}})\\}$.

6: Form $\mathcal{H}^{\boldsymbol{\epsilon}}_{m}$ and orthonormalize:
$\mathcal{H}^{\boldsymbol{\epsilon}}_{m}=\mathcal{H}_{m}+\mathcal{H}^{\mathbf{r}}_{m}$.

7: Using $\mathcal{H}^{\boldsymbol{\epsilon}}_{m}$, solve for the state error
estimation $\boldsymbol{\tilde{\epsilon}}(\omega)$ in (24b).

8: Choose the next sample $\omega^{*}$ from $\mathcal{B}$ as
$\omega^{*}=\arg\max\limits_{\omega\in\mathcal{B}}\|\boldsymbol{\tilde{\epsilon}}(\omega)\|_{\mathcal{H}}.$

9: Choose the next sample $\omega^{*}_{\boldsymbol{\epsilon}}$ from
$\mathcal{B}$ following
$\omega^{*}_{\boldsymbol{\epsilon}}=\arg\max\limits_{\omega\in\mathcal{B}}\|r(\boldsymbol{\tilde{\epsilon}}(\omega),\cdot;\omega)\|_{\mathcal{H}^{\prime}}.$
Here $r(\boldsymbol{\tilde{\epsilon}}(\omega),\cdot;\omega)$ is the residual
functional introduced by $\boldsymbol{\tilde{\epsilon}}(\omega)$ in (20),
namely,
$r(\boldsymbol{\tilde{\epsilon}}(\omega),\mathbf{v};\omega):=f^{\boldsymbol{\epsilon}}(\mathbf{v};\omega)-a(\boldsymbol{\tilde{\epsilon}}(\omega),\mathbf{v};\omega),\forall\mathbf{v}\in\mathcal{H}.$

10: $\xi=\|\boldsymbol{\tilde{\epsilon}}(\omega^{*})\|_{\mathcal{H}}.$

11: end while

12: Use $\mathcal{H}_{m}$ to solve (24a).

###  3.2 In-Band Eigenmodes in the Reduced-Basis Space

The works in [40, 10] suggest that a good approximation basis to the
frequency-parameter problem (14), i.e., a good reduced-basis space, should
include the resonant modes hit in the band of analysis $\mathcal{B}$. These
resonant modes constitute the dominant basis representing the electric field
$\mathbf{E}(\omega)$ in $\mathcal{B}$ (see (12)). In this work, we use this
basis not only for the original problem (14), but also for the residual
problem (20). For example,

$\displaystyle\mathcal{H}_{m}$
$\displaystyle=\text{span}\\{\mathbf{e}_{n}\leavevmode\nobreak\
|\leavevmode\nobreak\ \omega_{n}\in\leavevmode\nobreak\
\mathcal{B}\\}:=\mathcal{V}_{\mathcal{B}},$ (25a)
$\displaystyle\mathcal{H}^{\boldsymbol{\epsilon}}_{m}$
$\displaystyle=\mathcal{H}_{m}=\mathcal{V}_{\mathcal{B}}\text{.}$ (25b)

In this situation, we can find the solution for the reduced systems in (24) in
closed form and get some further insights, viz.,

$\displaystyle\mathbf{\tilde{E}}(\omega)=i\omega\sum\limits_{\omega_{n}^{2}\in\mathcal{B}_{2}}\frac{A_{n}}{1-\frac{\omega^{2}}{\omega_{n}^{2}}}\mathbf{e}_{n}\text{,
if }\omega^{2}\neq\omega_{n}^{2}\in\mathcal{B}_{2},$ (26a)
$\displaystyle\mathbf{\tilde{E}}(\omega)=i\omega\sum\limits_{\omega_{n}^{2}=\omega^{2}}a_{n}\mathbf{e}_{n}+i\omega\sum\limits_{\omega_{n}^{2}\in\mathcal{B}_{2}\setminus\\{\omega^{2}\\}}\frac{A_{n}}{1-\frac{\omega^{2}}{\omega_{n}^{2}}}\mathbf{e}_{n}\text{,}$
(26b) $\displaystyle\text{ if
}\omega^{2}=\omega_{n}^{2}\in\mathcal{B}_{2}\text{,}$ (26c)

and

$\displaystyle\boldsymbol{\tilde{\epsilon}}(\omega)=\mathbf{0}\text{, if
}\omega^{2}\neq\omega_{n}^{2}\in\mathcal{B}_{2},$ (27a)
$\displaystyle\boldsymbol{\tilde{\epsilon}}(\omega)=i\omega\sum\limits_{\omega_{n}^{2}=\omega^{2}}a^{\boldsymbol{\epsilon}}_{n}\mathbf{e}_{n}\text{,
if }\omega^{2}=\omega_{n}^{2}\in\mathcal{B}_{2}\text{.}$ (27b)

$\mathcal{B}_{2}$ stands for
$[\omega^{2}_{\text{min}},\omega^{2}_{\text{max}}]$. As (27) shows, this is
the worst-case scenario. Although the in-band eigenmode approximation basis in
the eigenspace $\mathcal{V}_{\mathcal{B}}$ is the best basis to capture the
fundamental dynamics of the electric field $\mathbf{E}(\omega)$ in the band of
analysis $\mathcal{B}$ and so might be the case for the residual system, it
turns out the opposite: the approximate state error
$\boldsymbol{\tilde{\epsilon}}(\omega)$ is, apart from the in-band resonances,
identically zero throughout the whole electromagnetic spectrum,
$\forall\omega\in\mathbb{R}$ (see (27)). Unfortunately, the actual error is
not zero but, on the contrary, we get zero error in (27). This is how Galerkin
approximation works, i.e., as far as the projection space is concerned, no
error can be identified within this space since all the error, which is not
identically zero, is orthogonal to the testing space and, therefore, remains
outside the projection space used in (25). As a result, we get zero error in
(27), apart from the in-band resonances. This is the rationale behind the use
of different projection spaces for the reduced problems (24) in [5], despite
the fact that two similar problems have to be solved independently at the same
time. Using two different projection spaces for the reduced systems increases
the chance to identify the true state error in the approximation. In this
work, we aim to achieve this goal while keeping the computational burden even
much lower.

###  3.3 Enhanced Reduced-Basis Space

In-band eigenmodes are not enough to ensure a good approximation to the
electric field $\mathbf{E}(\omega)$ in $\mathcal{B}$. In particular, only
resonant phenomena are strictly captured by the in-band eigenmode basis in the
eigenspace $\mathcal{V}_{\mathcal{B}}$, while other electromagnetic phenomena,
such as direct source to load couplings, are missing. This has already been
shown in Section 3.2. As a result, the reduced-basis space has to be enriched
by snapshots of the electric field $\mathbf{E}(w_{k})$ in the frequency band
$\mathcal{B}$, namely,

$\displaystyle\mathcal{H}_{m}$
$\displaystyle=\mathcal{V}_{\mathcal{B}}+\text{span}\\{\mathbf{E}(w_{k})\leavevmode\nobreak\
|\leavevmode\nobreak\ w_{k}\in\leavevmode\nobreak\ \mathcal{B}\\},$ (28a)
$\displaystyle\mathcal{V}_{\mathcal{B}}$
$\displaystyle=\text{span}\\{\mathbf{e}_{n}\leavevmode\nobreak\
|\leavevmode\nobreak\ \omega_{n}\in\leavevmode\nobreak\ \mathcal{B}\\},$ (28b)

giving rise to an enhanced reduced-basis space which ensures convergence to
the electric field $\mathbf{E}(\omega)$ in $\mathcal{B}$ within a fast setting
[10]. In our situation, one question still remains: what is the reduced-basis
space for the residual problem $\mathcal{H}^{\boldsymbol{\epsilon}}_{m}$ that
should be used to get an approximation $\boldsymbol{\tilde{\epsilon}}(\omega)$
to the state error while keeping a low computational effort? We have already
shown in Section 3.2 that using the same reduced-basis space considered for
the original problem, i.e.,
$\mathcal{H}^{\boldsymbol{\epsilon}}_{m}=\mathcal{H}_{m}$, while keeping the
computational burden low, is not a good choice (see (25)–(27)). Also, carrying
out a totally different reduced-basis space strategy for both reduced systems,
i.e., the original and residual systems, yields good approximation results but
the computational cost substantially increases [5]. We then follow a
compromise criterion: we allow the reduced-basis space for the residual
problem $\mathcal{H}^{\boldsymbol{\epsilon}}_{m}$ to be different from the
reduced-basis space for the original problem $\mathcal{H}_{m}$, but just for
only one additional basis vector. As a result, we keep the computational cost
low enough since only one additional solution needs to be carried out in the
whole process. This cheap approach may result in underestimation of
$\boldsymbol{\tilde{\epsilon}}(\omega)$ for the actual value of the state
error $\boldsymbol{\epsilon}(\omega)$. This is not critical, as will become
clear later. However, what is indeed essential is to choose an additional
basis vector that allows to monitor the state error by means of the reduced
residual system while not interfering the greedy algorithm with some unwanted
frequency modulation in the snapshot selection for the reduced original
problem. In other words, we should be careful to prevent the residual system
solution from biasing the greedy procedure to solve the reduced original
problem. Otherwise, the snapshots selected in the greedy process will end up
with a bad choice for approximation purposes.

We propose to add the stationary electric field $\mathbf{F}_{0}$ as the only
additional basis vector into $\mathcal{H}^{\boldsymbol{\epsilon}}_{m}$, as
$\mathbf{F}_{0}$ is missing in $\mathcal{H}_{m}$, so it contributes to the
error of the ROM. It is therefore reasonable to add $\mathbf{F}_{0}$ to the
reduced basis space $\mathcal{H}^{\boldsymbol{\epsilon}}_{m}$ for the
approximate error $\boldsymbol{\tilde{\epsilon}}(\omega)$. Furthermore, not
only is it orthogonal to the eigenspace $\mathcal{V}_{\mathcal{B}}$ but also
it has an almost-flat smooth influence on the fields $\mathbf{E}(\omega)$
throughout the frequency band of interest $\mathcal{B}$, due to its $1/\omega$
frequency behavior, cf. (12). As a consequence, its contribution to the error
is also small, not spoiling the desired property for
$\boldsymbol{\tilde{\epsilon}}(\omega)$, which should allow us to identify the
missing information in the reduced original system. Thus, we can ensure that
we are not modulating the adaptive sampling in the greedy algorithm. As a
result, the basis vector $\mathbf{F}_{0}$ added to the reduced-basis space for
the original problem $\mathcal{H}_{m}$ is a suitable candidate, since it has
the desired behavior as testing space
$\mathcal{H}^{\boldsymbol{\epsilon}}_{m}$ for the residual problem. Putting
everything together, we use the following spaces:

$\displaystyle\mathcal{H}_{m}$
$\displaystyle=\mathcal{V}_{\mathcal{B}}+\text{span}\\{\mathbf{E}(w_{k})\leavevmode\nobreak\
|\leavevmode\nobreak\ w_{k}\in\leavevmode\nobreak\ \mathcal{B}\\},$ (29a)
$\displaystyle\mathcal{V}_{\mathcal{B}}$
$\displaystyle=\text{span}\\{\mathbf{e}_{n}\leavevmode\nobreak\
|\leavevmode\nobreak\ \omega_{n}\in\leavevmode\nobreak\ \mathcal{B}\\},$ (29b)
$\displaystyle\mathcal{H}^{\boldsymbol{\epsilon}}_{m}$
$\displaystyle=\mathcal{H}_{m}+\text{span}\\{\mathbf{F}_{0}\\}\text{.}$ (29c)

Unfortunately, the actual value of the approximate state error
$\boldsymbol{\tilde{\epsilon}}(\omega)$ in this cheap statics-based reduced
residual problem solution, while accurately identifying the missing
information in the original problem, is _not_ a reliable indicator of the true
state error $\boldsymbol{\epsilon}(\omega)$ due to its implicit
underestimation, and can only be considered as a _rough_ estimator as compared
to the one in Algorithm 1. As a result, a different strategy should be carried
out to estimate the actual state error in the system.

We propose Algorithm 2 to carry out a fast _a posteriori_ state error
estimation for reliable frequency sweep analysis in microwave devices. The
reduced-basis space is adaptively built up by means of a greedy algorithm
based on this state error estimation. Further, to overcome the lack of
reliability of the rough error estimator
$\boldsymbol{\tilde{\epsilon}}(\omega)$ in case it is used as the stopping
criterion, a density stopping criterion [10] is preferred in Step 11. In
contrast to the error estimator $\boldsymbol{\tilde{\epsilon}}(\omega^{*})$ in
Step 10 of Algorithm 1, the true error at $\omega^{*}$ is computed in Step 11
of Algorithm 2 and is used as the stopping criterion. Note that the true error
$\mathbf{E}(\omega^{*})-\mathbf{\tilde{E}}(\omega^{*})$ at $\omega^{*}$ is
actually equivalent to the error $\mathbf{e}_{\perp}(\omega^{*})$ defined in
Algorithm 3 and is readily available from Step 10 in Algorithm 2, where
$\mathbf{E}(\omega^{*})$ is orthogonalized against the existing basis vectors
in $\mathcal{H}_{m}$ before being added to $\mathcal{H}_{m}$. In other words,
the new error indicator in Step 11 of Algorithm 2 is based on how much new
information the new greedy field sample $\mathbf{E}(\omega^{*})$ is adding to
the reduced-basis space. Eventually, the procedure stops whenever there is
nothing new to add, within a tolerance denoted by tol, and the reduced-basis
space can be considered as a dense enough approximation space which accurately
describes the field solution $\mathbf{E}(\omega)$ throughout the band of
interest $\mathcal{B}$.

In summary, we propose Algorithm 2 which makes use of two error indicators:
one is the rough error estimator $\boldsymbol{\tilde{\epsilon}}(\omega)$
employed in Step 8 based on a fixed $\mathcal{H}^{\mathbf{r}}_{m}$; the other
is the indicator used in Step 11, which is the true error at $\omega^{*}$
selected by the rough error estimator. It is then essential that these
adaptive greedy field snapshots are properly selected in the band of analysis.
Otherwise the procedure stops with no control on the actual error in the ROM.
Fortunately, the rough error estimator $\boldsymbol{\tilde{\epsilon}}(\omega)$
is capable of providing what is indeed the missing information in the system.
That means it catches the trend of the error change during the greedy
algorithm, although it may underestimate the actual error. Therefore, we
expect good performance in the proposed methodology. Finally, we present some
remarks about Algorithm 2.

Algorithm 2 Fast _a posteriori_ state error estimation for reliable frequency-
parameter MOR.

0: Frequency band of interest
$\mathcal{B}:=[\omega_{\text{min}},\omega_{\text{max}}]$, tolerance
$\texttt{tol}>0$ as the acceptable state error.

0: $\mathcal{H}_{m}$ to ensure tol state error in (24).

1: Solve for the in-band eigenmodes $\mathbf{e}_{n}$ in (14). Form the in-band
eigenspace $\mathcal{V}_{\mathcal{B}}$:
$\mathcal{V}_{\mathcal{B}}=\text{span}\\{\mathbf{e}_{n}\leavevmode\nobreak\
|\leavevmode\nobreak\ \omega_{n}\in\leavevmode\nobreak\ \mathcal{B}\\}$.

2: Solve for the stationary solution $\mathbf{F}_{0}$ in (14).

3: Initialize $\mathcal{H}_{m}=\mathcal{V}_{\mathcal{B}}$,
$\mathcal{H}^{\mathbf{r}}_{m}=\text{span}\\{\mathbf{F}_{0}\\}$,
$\xi=\texttt{tol}+1$. Choose a sample $\omega^{*}$ randomly taken from the end
points in $\mathcal{B}$.

4: Solve for $\mathbf{E}(\omega^{*})$ in (14) and update $\mathcal{H}_{m}$:
$\mathcal{H}_{m}=\mathcal{H}_{m}+\text{span}\\{\mathbf{E}(\omega^{*})\\}$.

5: while $\xi>\texttt{tol}$ do

6: Compute the approximate field $\mathbf{\tilde{E}}(\omega)$ in (24a).

7: Form $\mathcal{H}^{\boldsymbol{\epsilon}}_{m}$:
$\mathcal{H}^{\boldsymbol{\epsilon}}_{m}=\mathcal{H}_{m}+\mathcal{H}^{\mathbf{r}}_{m}$.

8: Compute the state error estimation $\boldsymbol{\tilde{\epsilon}}(\omega)$
in (24b).

9: Choose the next sample $\omega^{*}$ from $\mathcal{B}$ as
$\omega^{*}=\arg\max\limits_{\omega\in\mathcal{B}}\|\boldsymbol{\tilde{\epsilon}}(\omega)\|_{\mathcal{H}}.$

10: Compute $\mathbf{E}(\omega^{*})$ in (14) and update $\mathcal{H}_{m}$:
$\mathcal{H}_{m}=\mathcal{H}_{m}+\text{span}\\{\mathbf{E}(\omega^{*})\\}$.

11:
$\xi=\|\mathbf{e}_{\perp}(\omega^{*})\|_{\mathcal{H}}=\|\mathbf{E}(\omega^{*})\|_{\nu}.$

12: end while

13: Use $\mathcal{H}_{m}$ to solve (24a).

Algorithm 3 Evaluation of the new information added to the reduced-basis space
$\mathcal{H}_{m}$ by the field sample $\mathbf{E}$ [10].

0: Electric field $\mathbf{E}$ and reduced-basis space $\mathcal{H}_{m}$.

0: Norm $\|\mathbf{E}\|_{\nu}$ of the field $\mathbf{E}$ indicating the new
information added to the reduced-basis space $\mathcal{H}_{m}$.

1: Normalize $\mathbf{E}$:
$\mathbf{e}=\mathbf{E}/\|\mathbf{E}\|_{\mathcal{H}}$.

2: Project $\mathbf{e}$ onto $\mathcal{H}_{m}$. Decompose $\mathbf{e}$ into
its projection and its orthogonal complement onto $\mathcal{H}_{m}$:
$\mathbf{e}=\mathbf{e}_{\parallel}+\mathbf{e}_{\perp}$.

3: Set
$\|\mathbf{E}\|_{\nu}=\|\mathbf{e}-\mathbf{e}_{\parallel}\|_{\mathcal{H}}=\|\mathbf{e}_{\perp}\|_{\mathcal{H}}$.

######  Remark 1.

Once the in-band eigenmode basis is completed, a randomly chosen sample
between the end points in the band of analysis $\mathcal{B}$ is taken to
enrich the eigenbasis in Step 3. At this point, neither the residual norm nor
the state error norm can provide an answer to what is the best sample to
choose, since both estimators show a constant behavior throughout the whole
spectrum $\forall\omega\in\mathbb{R}$. We refer to [10] for the illustration
of the residual norm behavior when eigenbasis is used. As a result, based on
the new information arguments added to the eigenbasis, any of the end points
in the frequency band of interest should be sampled (see (12)). It should be
noted that, in (12), the further we are from the in-band eigenmodes, the more
linearly independent new information is found, until a new out of band
eigenmode eventually shows up. This is the main reason for sampling at the end
points of the frequency band of interest $\mathcal{B}$.

######  Remark 2.

Contrary to what is done in [10], the residual norm is not used to guide the
greedy algorithm at any step. Thus, better results are expected. Residual
information is problematic since its behavior suffers from the in-band
eigenresonances. The residual does not vanish at and nearby resonances, which
can definitely mislead the greedy algorithm. This will become apparent in
Section 4 throughout the numerical examples.

##  4 Numerical Results

In this section, we apply the proposed _a posteriori_ state error estimation
for reliable fast frequency sweeps of different challenging microwave
circuits, namely, a quad-mode dielectric resonator filter, an inline filter
with transmission zeros generated by frequency-dependent couplings, an inline
dielectric resonator filter and a combline diplexer. The capabilities and
reliability of the proposed procedure are demonstrated via these examples. The
in-house C++ code for finite element method (FEM) simulations uses a second-
order first family of Nédélec’s elements [30, 22], on meshes provided by Gmsh
[15]. All computations were carried out on a workstation with 3.00-GHz Intel
Xeon E5-2687W v4 processor and 256-GB RAM.

In our experiments, we define the true error ($\epsilon_{\text{true}}$) of the
ROM as the maximal error over the whole frequency band $\mathcal{B}$ using the
indicator $\mathbf{e}_{\perp}$ in Algorithm 3, namely,

$\epsilon_{\text{true}}=\max\limits_{\omega\in\mathcal{B}}\|\mathbf{e}_{\perp}(\omega)\|_{\mathcal{H}}\text{.}$
(30)

Note that the novelty of the proposed Algorithm 2 is twofold: 1) a rough error
estimator $\boldsymbol{\tilde{\epsilon}}(\omega)$ to save computational costs
as compared with the one in Algorithm 1; 2) the indicator
$\|\mathbf{e}_{\perp}(\omega)\|_{\mathcal{H}}$ to improve the reliability of
the rough estimator. To show that the proposed Algorithm 2 is more efficient,
we compare it with Algorithm 1, as well as with the greedy algorithm using the
residual norm-based error estimator. To this end, we define the indicators
$\|\mathbf{e}_{\perp}(\omega)\|_{\mathcal{H}}$ based on the above three
different estimators, respectively, i.e.,

$\epsilon_{\text{state}}=\|\mathbf{e}_{\perp}(\arg\max\limits_{\omega\in\mathcal{B}}\|\boldsymbol{\tilde{\epsilon}}(\omega)\|_{\mathcal{H}})\|_{\mathcal{H}}\text{,}$
(31)

where $\boldsymbol{\tilde{\epsilon}}(\omega)$ refers to either the one in
Algorithm 1 or the rough estimator in Algorithm 2. The indicator based on the
residual-norm is defined as

$\epsilon_{\text{res}}=\|\mathbf{e}_{\perp}(\arg\max\limits_{\omega\in\mathcal{B}}\|r(\mathbf{\tilde{E}}(\omega),\cdot;\omega)\|_{\mathcal{H}^{\prime}})\|_{\mathcal{H}}\text{.}$
(32)

For a fair comparison, we compare the results of Algorithm 2 with Algorithm 1,
where $\xi$ in Step 10 is replaced by the indicator $\epsilon_{\text{state}}$
in (31). We also compare Algorithm 2 with the residual-norm based greedy
algorithm, where the indicator in (32) is used as the stopping criterion. The
comparison between Algorithm 1 and Algorithm 2 is carried out for the first
two examples. The indicator $\mathbf{e}_{\perp}$ plays a role of paramount
importance in providing a fair comparison among the different strategies. It
should be pointed out that the true error $\epsilon_{\text{true}}$ in (30)
implies the computation of the field solution $\mathbf{E}(\omega)$ by means of
time-consuming FEM simulations throughout the whole frequency band of interest
$\mathcal{B}$. This can _only_ be carried out for academic purposes. Finally,
as a figure of merit, we use the metric of effectivity to gauge how close the
estimated error is to the true error:
$\texttt{eff}:=\frac{\epsilon}{\epsilon_{\text{true}}}$, $\epsilon$ refers to
either $\epsilon_{\text{state}}$ in (31) or $\epsilon_{\text{res}}$ in (32).
For each of the four microwave circuits considered, we evaluate the
performance of the error estimators using the indicators defined in (31) and
(32). A tolerance threshold of $\texttt{tol}=2\cdot 10^{-7}$ is used
throughout all the numerical examples, for all greedy algorithms corresponding
to their respective error estimators.

###  4.1 Quad-Mode Dielectric Resonator Filter

A quad-mode dielectric resonator filter in a single cylindrical cavity is
proposed in [27]. Fig. 1 shows the geometry of the filter as well as the mesh
used for its analysis. These structures are extremely attractive since
multiple resonant modes show up in a single cylindrical cavity due to the
dielectric resonator. At the same time, they are difficult to tune since all
dominant modes are coupled with each other, requiring multiple full-wave
electromagnetic analyses to carry out a good electrical design. Six tuning
screws are included in this filter. It is then of paramount importance to
accurately predict the electromagnetic behavior in a frequency band in an
efficient way. The filter detailed in Fig. 1 is a quad-mode filter, where a
four-pole passband filtering response, shown in Fig. 2, is obtained. However,
at the same time, there are additional resonant modes in the band of analysis,
$\mathcal{B}:=[3.4,4.2]$ GHz, giving rise to a more complicated response
including direct source to load coupling, affecting the position of the two
transmission zeros, rather than a typical four-pole frequency response. As a
result, a reliable ROM for fast frequency sweep analysis is essential.

An FEM system with 245,778 degrees of freedom is used to solve for the
electric field. Following Algorithm 2, a ROM is obtained by means of the
Reduced-Basis Method (RBM) giving rise to a reduced system of dimension 14 to
compute the frequency response detailed in Fig. 2. It is clear that there is
not even a competition ($14\ll 245,778$) between the computational effort
carried out by RBM and the one that would have been required by FEM to get the
same frequency response along the same frequency samples, nevermind
subsampling is considered. Solving a reduced system of dimension 14 many times
is totally effortless while solving many times an FEM system of dimension
245,778 is rather time-consuming. Good agreement is found between the FEM and
RBM results in Fig. 2.

Next, a comparison of the different MOR methodologies is carried out. Table 1
not only details the frequency samples adaptively chosen by each greedy
algorithm but also the estimated and true errors, at each iteration in the MOR
process. Fig. 3 depicts the convergence behavior for the different
methodologies and details the effectivity metric for the proposed _a
posteriori_ state error estimator. Good behavior is observed in the proposed
methodology in comparison to the true error. On the contrary, the residual
norm-based error estimation results in a greedy algorithm that underestimates
the error and prematurely stops the iterative procedure. The rationale behind
this is shown in Fig. 4. For a ROM of dimension 12 obtained by the proposed
approach, both the residual norm
$\|r(\mathbf{\tilde{E}}(\omega),\cdot;\omega)\|_{\mathcal{H}^{\prime}}$ and
the approximate state error
$\|\boldsymbol{\tilde{\epsilon}}(\omega)\|_{\mathcal{H}}$ are plotted versus
frequency. While the state error estimation has a smooth frequency behavior,
the residual norm is contaminated by the in-band resonant modes, creating an
eigenfrequency pollution. This makes the residual norm-based greedy algorithm
mislead its sampling once again around the eigenresonances, notoriously
deteriorating the new information added to the reduced-basis space. This is
the key advantage of using state error estimation, where this unwanted
behavior does not hold.

Finally, Table 2 compares the performance of the proposed approach with
Algorithm 1, where $\xi$ in Step 10 is replaced by the indicator in (31). It
should be pointed out that, even though the in-band eigenmodes are not imposed
in Algorithm 1, these in-band eigenmodes do show up in the greedy algorithm
after the first few iterations in the procedure. This indicates Algorithm 1 is
accurately working out since the in-band eigenmodes have been shown to be the
best choice from the theoretical point of view in Section 2. As far as
computational complexity is concerned, Algorithm 1 is more expensive than
Algorithm 2. As a matter of fact, the size of both ROMs when the same stopping
criterion is used is $14$ for Algorithm 2 and $16$ for Algorithm 1, whereas
the size of the residual ROMs defined in (24b) is $15$ and $30$ for Algorithms
2 and 1, respectively. This shows how many more efforts we need to carry out
in Algorithm 1 compared to Algorithm 2. As a result, the proposed approach
yields a fast _a posteriori_ state error estimation.

Table 1: Greedy Algorithm Frequency Samples for Different MOR Methodologies in the Quad-Mode Filter. GHz | $\epsilon_{\text{true}}$ | GHz | $\epsilon_{\text{state}}$ | GHz | $\epsilon_{\text{res}}$
---|---|---|---|---|---
3.6103 | 1. | 3.6103 | 1. | 3.6103 | 1.
3.6245 | 1. | 3.6245 | 1. | 3.6245 | 1.
3.6887 | 1. | 3.6887 | 1. | 3.6887 | 1.
3.7149 | 1. | 3.7149 | 1. | 3.7149 | 1.
3.9299 | 1. | 3.9299 | 1. | 3.9299 | 1.
4.0067 | 1. | 4.0067 | 1. | 4.0067 | 1.
4.1552 | 1. | 4.1552 | 1. | 4.1552 | 1.
4.2000 | $1.7\cdot 10^{-1}$ | 3.4000 | $1.6\cdot 10^{-1}$ | 3.4000 | $1.6\cdot 10^{-1}$
3.4000 | $7.1\cdot 10^{-2}$ | 4.2000 | $8.0\cdot 10^{-2}$ | 4.0067 | $2.8\cdot 10^{-7}$
3.7800 | $1.6\cdot 10^{-3}$ | 3.8820 | $8.3\cdot 10^{-4}$ | – | –
4.1000 | $9.4\cdot 10^{-5}$ | 4.0940 | $6.7\cdot 10^{-5}$ | – | –
3.4900 | $6.8\cdot 10^{-6}$ | 3.5930 | $2.5\cdot 10^{-6}$ | – | –
4.1800 | $1.2\cdot 10^{-6}$ | 4.1680 | $5.8\cdot 10^{-7}$ | – | –
3.8600 | $1.9\cdot 10^{-7}$ | 3.6930 | $1.5\cdot 10^{-7}$ | – | –
Table 2: State Error Algorithm Comparison in the Quad-Mode Filter. Algorithm 2 | Algorithm 1
---|---
GHz | $\epsilon_{\text{state}}$ | GHz | $\epsilon_{\text{state}}$
3.6103 | 1. | 3.4000 | $1.6\cdot 10^{-3}$
3.6245 | 1. | 3.6470 | $1.7\cdot 10^{-3}$
3.6887 | 1. | 4.0810 | $1.6$
3.7149 | 1. | 3.9300 | $2.6\cdot 10^{-3}$
3.9299 | 1. | 3.6250 | $1.2\cdot 10^{-3}$
4.0067 | 1. | 3.7150 | $3.2\cdot 10^{-4}$
4.1552 | 1. | 4.1550 | $9.8\cdot 10^{-7}$
3.4000 | $1.6\cdot 10^{-1}$ | 4.0070 | $5.5\cdot 10^{-8}$
4.2000 | $8.0\cdot 10^{-2}$ | – | –
3.8820 | $8.3\cdot 10^{-4}$ | – | –
4.0940 | $6.7\cdot 10^{-5}$ | – | –
3.5930 | $2.5\cdot 10^{-6}$ | – | –
4.1680 | $5.8\cdot 10^{-7}$ | – | –
3.6930 | $1.5\cdot 10^{-7}$ | – | –
ROM size | 14 | ROM size | 16
Residual ROM size | 15 | Residual ROM size | 30
Figure 1: Quad-mode dielectric resonator filter proposed in [27].
$3.4$$3.5$$3.6$$3.7$$3.8$$3.9$$4$$4.1$$4.2$$-80$$-60$$-40$$-20$$0$Frequency
(GHz)dB$|\text{S}_{\text{21}}|$$|\text{S}_{\text{11}}|$ Figure 2: Quad-mode
dielectric resonator filter scattering parameter response comparison with the
proposed approach. (–) RBM. ($\circ$) FEM.

$2$$4$$6$$8$$10$$12$$14$$10^{0}$$10^{-2}$$10^{-4}$$10^{-6}$IterationsError$\epsilon_{\text{res}}$$\epsilon_{\text{state}}$$\epsilon_{\text{true}}$
(a)

$2$$4$$6$$8$$10$$12$$14$$0.4$$0.6$$0.8$$1$$1.2$Iterationseff$\epsilon_{\text{state}}$
(b)

Figure 3: Quad-mode dielectric resonator filter error estimator results. (a)
Convergence of the greedy algorithm. (b) Effectivity (eff).
$3.4$$3.5$$3.6$$3.7$$3.8$$3.9$$4$$4.1$$4.2$$10^{-12}$$10^{-11}$$10^{-10}$$10^{-9}$$10^{-8}$$10^{-7}$Frequency
(GHz)Error Estimation Figure 4: Quad-mode dielectric resonator filter error
estimation frequency behavior for the ROM of order 12. (–) Residual error
$\|r(\mathbf{\tilde{E}}(\omega),\cdot;\omega)\|_{\mathcal{H}^{\prime}}$. (–)
State error estimator
$\|\boldsymbol{\tilde{\epsilon}}(\omega)\|_{\mathcal{H}}$.

###  4.2 Inline Filter with Frequency-Dependent Couplings

The next example is a fourth-order inline combline filter designed in [17],
where frequency dependent couplings are taken into account to provide finite
transmission zeros, even within an inline coupling route structure. The filter
geometry is depicted in Fig. 5 and the mesh for the electromagnetic analysis
is shown as well. Within four combline resonant cavities, two additional shunt
inductors and capacitors are included, which indeed give rise to additional
higher frequency resonances creating two finite transmission zeros near by the
filter passband. As a result, better rejection properties are allowed in this
inline structure. Frequency dependent coupling filter theory is continuously
developing [39, 41] and a large number of full-wave analyses are required to
design these kind of filtering responses. It is essential to carry out a
reliable fast frequency sweep analysis to meet the target electrical response
in the design optimization loop.

The frequency band $\mathcal{B}:=[1.4,2.9]$ GHz is taken into account for
analysis. We solve for the electric field in the band of interest
$\mathcal{B}$ by means of an FEM system with 105,690 degrees of freedom.
Following the procedure proposed in this work, we obtain a ROM of dimension 13
to sweep the frequency response of this filter via RBM. A comparison between
the filter response results obtained by FEM and RBM is shown in Fig. 6.
Reasonable agreement is achieved.

Next, we carry out the numerical tests to compare the different MOR
methodologies. Table 3 details the frequency samples adaptively chosen by each
greedy algorithm as well as the error estimator and true error at each
iteration in the MOR process. Contrary to what Algorithm 2 proposes, it should
be noted that the true error $\epsilon_{\text{true}}$-based greedy algorithm
does not choose an end point sample after the eigenbasis made of 6 in-band
resonant modes is built up. Fig. 7 depicts the convergence behavior for the
different methodologies and details the effectivity metric for the proposed _a
posteriori_ state error estimator. Reasonable performance is observed in the
proposed methodology in comparison to the true error. Once again, the residual
norm-based greedy algorithm prematurely aborts the iterative MOR procedure.

In addition, Table 4 compares the performance of the proposed approach to
Algorithm 1 where $\xi$ in Step 10 is replaced by the indicator in (31). Once
again, even though the in-band eigenmodes are not _a priori_ included in the
projection basis in Algorithm 1, these show up in the greedy algorithm after
the first few iterations. This shows Algorithm 1 is properly working out since
the in-band eigenmodes are expected to show up in a good approximation basis,
as has been discussed from a theoretical point of view. As far as
computational burden is concerned, Algorithm 1 is more time consuming than
Algorithm 2. As a matter of fact, the size of both ROMs when the same stopping
criterion is used is $12$ for Algorithm 2 and $14$ for Algorithm 1, whereas
the size of the residual ROMs defined in (24b) is $13$ and $27$ for Algorithms
2 and 1, respectively. This shows the additional effort that has to be carried
out in Algorithm 1 with respect to Algorithm 2. As a result, the proposed
approach gives rise to a fast _a posteriori_ state error estimation.

Table 3: Greedy Algorithm Frequency Samples for Different MOR Methodologies in the Inline Filter with Finite Transmission Zeros. GHz | $\epsilon_{\text{true}}$ | GHz | $\epsilon_{\text{state}}$ | GHz | $\epsilon_{\text{res}}$
---|---|---|---|---|---
1.6981 | 1. | 1.6981 | 1. | 1.6981 | 1.
1.7043 | 1. | 1.7043 | 1. | 1.7043 | 1.
1.7799 | 1. | 1.7799 | 1. | 1.7799 | 1.
1.8328 | 1. | 1.8328 | 1. | 1.8328 | 1.
2.7661 | 1. | 2.7661 | 1. | 2.7661 | 1.
2.7856 | 1. | 2.7856 | 1. | 2.7856 | 1.
2.3000 | $8.8\cdot 10^{-1}$ | 1.4000 | $6.7\cdot 10^{-1}$ | 1.4000 | $6.7\cdot 10^{-1}$
2.9000 | $7.7\cdot 10^{-2}$ | 2.9000 | $1.5\cdot 10^{-1}$ | 2.9000 | $1.5\cdot 10^{-1}$
1.4000 | $4.6\cdot 10^{-3}$ | 2.5156 | $2.8\cdot 10^{-3}$ | 2.7661 | $8.4\cdot 10^{-8}$
2.6200 | $2.9\cdot 10^{-4}$ | 2.0520 | $3.5\cdot 10^{-4}$ | – | –
1.9800 | $9.8\cdot 10^{-6}$ | 2.8187 | $7.6\cdot 10^{-6}$ | – | –
2.8500 | $4.9\cdot 10^{-7}$ | 1.7007 | $9.4\cdot 10^{-9}$ | – | –
1.7600 | $1.7\cdot 10^{-7}$ | – | – | – | –
Table 4: State Error Algorithm Comparison in the Inline Filter with Finite Transmission Zeros. Algorithm 2 | Algorithm 1
---|---
GHz | $\epsilon_{\text{state}}$ | GHz | $\epsilon_{\text{state}}$
1.6981 | 1. | 1.4000 | $3.6\cdot 10^{-3}$
1.7043 | 1. | 1.9430 | $1.6\cdot 10^{-2}$
1.7799 | 1. | 1.7260 | $6.1$
1.8328 | 1. | 2.7660 | $1.2\cdot 10^{-2}$
2.7661 | 1. | 1.6980 | $6.3\cdot 10^{-4}$
2.7856 | 1. | 2.7860 | $2.5\cdot 10^{-6}$
1.4000 | $6.7\cdot 10^{-1}$ | 1.8330 | $3.6\cdot 10^{-8}$
2.9000 | $1.5\cdot 10^{-1}$ | – | –
2.5156 | $2.8\cdot 10^{-3}$ | – | –
2.0520 | $3.5\cdot 10^{-4}$ | – | –
2.8187 | $7.6\cdot 10^{-6}$ | – | –
1.7007 | $9.4\cdot 10^{-9}$ | – | –
ROM size | 12 | ROM size | 14
Residual ROM size | 13 | Residual ROM size | 27
Figure 5: Inline filter with finite transmission zeros proposed in [17].
$1.4$$1.6$$1.8$$2$$2.2$$2.4$$2.6$$2.8$$-100$$-80$$-60$$-40$$-20$$0$Frequency
(GHz)dB$|\text{S}_{\text{21}}|$$|\text{S}_{\text{11}}|$ Figure 6: Inline
filter with finite transmission zeros scattering parameter response comparison
with the proposed approach. (–) RBM. ($\circ$) FEM.

$2$$4$$6$$8$$10$$12$$10^{0}$$10^{-2}$$10^{-4}$$10^{-6}$$10^{-8}$IterationsError$\epsilon_{\text{res}}$$\epsilon_{\text{state}}$$\epsilon_{\text{true}}$
(a)

$2$$4$$6$$8$$10$$12$$0$$1$$2$Iterationseff$\epsilon_{\text{state}}$ (b)

Figure 7: Inline filter with finite transmission zeros error estimator
results. (a) Convergence of the greedy algorithm. (b) Effectivity (eff).

###  4.3 Inline Dielectric Resonator Filter

A sixth order inline dielectric resonator filter with two transmission zeros
is depicted in Fig. 8. Cross-coupling between nonadjacent dielectric
resonators, appropriately arranging their orientations, is obtained by
exploiting multiple evanescent modes in the inline structure. This filter is
proposed in [3]. The $[2.14,2.2]$ GHz band is taken into account in the
analysis. A FEM discretization shown in Fig. 8 with 230,058 degrees of freedom
is used. Fig. 9 details the filter response in the band of analysis under both
time-consuming FEM simulation and fast RBM analysis. A ROM of dimension 10 is
used to get the fast frequency sweep results in Fig. 9. Good agreement is
obtained between both analyses. It should be pointed out that the FEM solution
of the FOM evaluated at a given frequency takes $5.610$ seconds. In the online
stage, the ROM resulting from RBM was evaluated for $1201$ different frequency
samples requiring $0.085$ seconds in total; this works out to $70$
microseconds to solve a single ROM and amounts to a speedup of nearly $80000$.

Next, we compare the different MOR techniques. Table 5 shows the frequency
samples adaptively chosen by each greedy algorithm as well as the error
estimator at each iteration in the MOR process. As expected in a sixth order
filter, 6 in-band eigenmodes are found in the frequency band of analysis. Fig.
10 details the convergence behavior for the different methodologies and shows
the effectivity metric for the proposed approach. A good behavior is observed
in the proposed _a posteriori_ state error estimator.

Table 5: Greedy Algorithm Frequency Samples for Different MOR Methodologies in the Inline Dielectric Resonator Filter. GHz | $\epsilon_{\text{true}}$ | GHz | $\epsilon_{\text{state}}$ | GHz | $\epsilon_{\text{res}}$
---|---|---|---|---|---
2.1635 | 1. | 2.1635 | 1. | 2.1635 | 1.
2.1640 | 1. | 2.1640 | 1. | 2.1640 | 1.
2.1663 | 1. | 2.1663 | 1. | 2.1663 | 1.
2.1709 | 1. | 2.1709 | 1. | 2.1709 | 1.
2.1768 | 1. | 2.1768 | 1. | 2.1768 | 1.
2.1788 | 1. | 2.1788 | 1. | 2.1788 | 1.
2.1400 | $3.4\cdot 10^{-1}$ | 2.1400 | $3.4\cdot 10^{-1}$ | 2.1400 | $3.4\cdot 10^{-1}$
2.2000 | $7.8\cdot 10^{-3}$ | 2.2000 | $7.8\cdot 10^{-3}$ | 2.1709 | $5.4\cdot 10^{-8}$
2.1530 | $2.8\cdot 10^{-5}$ | 2.1697 | $5.7\cdot 10^{-6}$ | – | –
2.1915 | $3.6\cdot 10^{-7}$ | 2.1874 | $1.7\cdot 10^{-7}$ | – | –
2.1785 | $1.6\cdot 10^{-7}$ | – | – | – | –
Figure 8: Inline dielectric resonator filter designed in [3].
$2.14$$2.15$$2.15$$2.16$$2.16$$2.17$$2.17$$2.18$$2.18$$2.19$$2.19$$2.2$$2.2$$-100$$-80$$-60$$-40$$-20$$0$Frequency
(GHz)dB$|\text{S}_{\text{21}}|$$|\text{S}_{\text{11}}|$ Figure 9: Inline
dielectric resonator filter scattering parameter response comparison with the
proposed approach. (–) RBM. ($\circ$) FEM.

$2$$4$$6$$8$$10$$10^{0}$$10^{-2}$$10^{-4}$$10^{-6}$IterationsError$\epsilon_{\text{res}}$$\epsilon_{\text{state}}$$\epsilon_{\text{true}}$
(a)

$2$$4$$6$$8$$10$$0.2$$0.4$$0.6$$0.8$$1$Iterationseff$\epsilon_{\text{state}}$
(b)

Figure 10: Inline dielectric resonator filter error estimator results. (a)
Convergence of the greedy algorithm. (b) Effectivity (eff).

###  4.4 Combline Diplexer

The last real-life application is an $11^{\text{th}}$ order combline diplexer
with star-junction designed in [44]. The geometry of this diplexer is shown in
Fig. 11. The frequency band of analysis is $\mathcal{B}:=[2.2,3.0]$ GHz. An
FEM system with 270,446 degrees of freedom arises whereas the application of
the proposed methodology in Algorithm 2 gives rise to a ROM of dimension 20 by
means of RBM. The scattering parameter response for this diplexer is detailed
in Fig. 12. Good agreement is found between FEM and RBM results. It should be
pointed out that further tuning is needed to obtain the target equiripple
response. The FOM solution for this example requires $4.630$ seconds. In the
online stage, the ROM was evaluated at $1601$ different frequencies taking
$0.446$ seconds. Thus, the time to solve a single ROM is $278$ microseconds,
at a speedup of around $16650$.

As expected in an $11^{\text{th}}$ order diplexer, 11 in-band eigenmodes are
found in the frequency band of interest. Table 6 details the different
frequency samples for each methodology. A comparison for the different MOR
strategies is shown in Fig. 13, where the convergence of the estimated and
true errors at each iteration as well as the effectivity metric is detailed.
Once again, the residual norm-based greedy algorithm prematurely stops, misled
by oversampling nearby the eigenresonances. A reasonable performance is
observed in the proposed _a posteriori_ state error estimator.

Table 6: Greedy Algorithm Frequency Samples for Different MOR Methodologies in the Combline Diplexer. GHz | $\epsilon_{\text{true}}$ | GHz | $\epsilon_{\text{state}}$ | GHz | $\epsilon_{\text{res}}$
---|---|---|---|---|---
2.4623 | 1. | 2.4623 | 1. | 2.4623 | 1.
2.4914 | 1. | 2.4914 | 1. | 2.4914 | 1.
2.5233 | 1. | 2.5233 | 1. | 2.5233 | 1.
2.5592 | 1. | 2.5592 | 1. | 2.5592 | 1.
2.5803 | 1. | 2.5803 | 1. | 2.5803 | 1.
2.5981 | 1. | 2.5981 | 1. | 2.5981 | 1.
2.6163 | 1. | 2.6163 | 1. | 2.6163 | 1.
2.6419 | 1. | 2.6419 | 1. | 2.6419 | 1.
2.6766 | 1. | 2.6766 | 1. | 2.6766 | 1.
2.7102 | 1. | 2.7102 | 1. | 2.7102 | 1.
2.7274 | 1. | 2.7274 | 1. | 2.7274 | 1.
3.0000 | $8.6\cdot 10^{-1}$ | 3.0000 | $8.6\cdot 10^{-1}$ | 3.0000 | $8.6\cdot 10^{-1}$
2.2000 | $4.5\cdot 10^{-1}$ | 2.2000 | $4.5\cdot 10^{-1}$ | 2.2000 | $4.5\cdot 10^{-1}$
2.8600 | $6.0\cdot 10^{-2}$ | 2.5094 | $1.1\cdot 10^{-2}$ | 2.6766 | $3.1\cdot 10^{-9}$
2.3500 | $7.2\cdot 10^{-3}$ | 2.9067 | $1.1\cdot 10^{-2}$ | – | –
2.9500 | $1.1\cdot 10^{-3}$ | 2.6060 | $9.2\cdot 10^{-5}$ | – | –
2.2600 | $1.6\cdot 10^{-4}$ | 2.2835 | $2.4\cdot 10^{-4}$ | – | –
2.7800 | $1.5\cdot 10^{-5}$ | 2.9712 | $3.7\cdot 10^{-5}$ | – | –
2.4200 | $1.1\cdot 10^{-6}$ | 2.2137 | $7.9\cdot 10^{-7}$ | – | –
2.4100 | $1.9\cdot 10^{-7}$ | 2.6766 | $1.1\cdot 10^{-8}$ | – | –
Figure 11: Combline diplexer designed in [44].
$2.2$$2.3$$2.4$$2.5$$2.6$$2.7$$2.8$$2.9$$3$$-100$$-80$$-60$$-40$$-20$$0$Frecuency
(GHz)dB$|\text{S}_{\text{11}}|$$|\text{S}_{\text{21}}|$$|\text{S}_{\text{31}}|$
Figure 12: Combline diplexer scattering parameter response comparison with the
proposed approach. (–) RBM. ($\circ$) FEM.

$2$$5$$8$$11$$14$$17$$20$$10^{0}$$10^{-2}$$10^{-4}$$10^{-6}$$10^{-8}$IterationsError$\epsilon_{\text{res}}$$\epsilon_{\text{state}}$$\epsilon_{\text{true}}$
(a)

$2$$5$$8$$11$$14$$17$$20$$0$$1$$2$Iterationseff$\epsilon_{\text{state}}$ (b)

Figure 13: Combline diplexer error estimator results. (a) Convergence of the
greedy algorithm. (b) Effectivity (eff).

##  5 Conclusions

A compact and reliable MOR method for fast frequency sweeps in microwave
circuits by means of the reduced-basis method has been detailed. A compact
basis including the in-band resonant modes hit in the frequency band of
interest for both the reduced basis approximation of the electric field and
its corresponding state error has been proposed. This allows to efficiently
solve for both the reduced original and residual problems, thus minimizing the
additional computational effort. The benefits of using proper state error
estimators avoiding time-consuming _inf-sup_ constant evaluations has also
been highlighted. As a result, a fast _a posteriori_ state error estimator for
the ROM has been obtained. Real-life microwave devices, including a quad-mode
dielectric resonator filter and a combline diplexer, have shown the
capabilities and reliability of the proposed methodology.

## Acknowledgments

Sridhar Chellappa is supported by the International Max Planck Research School
for Advanced Methods in Process and Systems Engineering (IMPRS-ProEng).

## References

  * [1] L. Balewski, G. Fotyga, M. Mrozowski, M. Mul, P. Sypek, D. Szypulski, and A. Lamecki.  Step on it! Bringing fullwave finite-element microwave filter design up to speed.  IEEE Microwave Magazine, 21(3):34–49, 2020.  doi:10.1109/MMM.2019.2958165.
  * [2] R. Baltes, A. Schultschik, O. Farle, and R. Dyczij-Edlinger.  A finite-element-based fast frequency sweep framework including excitation by frequency-dependent waveguide mode patterns.  IEEE Trans. Microw. Theory Techn., 65(7):2249–2260, 2017.  doi:10.1109/TMTT.2017.2679181.
  * [3] Simone Bastioli and Richard V Snyder.  Inline pseudoelliptic $\text{TE}_{01\delta}$-mode dielectric resonator filters using multiple evanescent modes to selectively bypass orthogonal resonators.  IEEE Trans. Microw. Theory Techn., 60(12):3988–4001, 2012.  doi:10.1109/TMTT.2012.2222659.
  * [4] Peter Benner, Serkan Gugercin, and Karen Willcox.  A survey of projection-based model reduction methods for parametric dynamical systems.  SIAM Rev., 57(4):483–531, 2015.  doi:10.1137/130932715.
  * [5] S. Chellappa, L. Feng, V. de la Rubia, and P. Benner.  Inf-sup-constant-free state error estimator for model order reduction of parametric systems in electromagnetics.  e-prints 2104.12802, arXiv, 2021.  math.NA.  URL: https://arxiv.org/abs/2104.12802.
  * [6] Lorenzo Codecasa, Gian Guido Gentili, and Marco Politi.  Exploiting port responses for wideband analysis of multimode lossless devices.  IEEE Trans. Microw. Theory Techn., 68(2):555–563, 2020.  doi:10.1109/TMTT.2019.2952853.
  * [7] G. Conciauro, M. Guglielmi, and R. Sorrentino.  Advanced Modal Analysis.  Wiley, New York, NY, USA, 2000.  doi:10.5555/554736.
  * [8] Qi I. Dai, Y. H. Lo, W. C. Chew, Y. G. Liu, and L. J. Jiang.  Generalized modal expansion and reduced modal representation of 3-D electromagnetic fields.  IEEE Trans. Antennas Propag., 62(2):783–793, 2014.  doi:10.1109/TAP.2013.2292083.
  * [9] Valentín de la Rubia.  Reliable reduced-order model for fast frequency sweep in microwave circuits.  Electromagnetics, 34(3-4):161–170, 2014.  doi:10.1080/02726343.2014.877735.
  * [10] Valentín de la Rubia and Michal Mrozowski.  A compact basis for reliable fast frequency sweep via the reduced-basis method.  IEEE Trans. Microw. Theory Techn., 66(10):4367–4382, 2018.  doi:10.1109/TMTT.2018.2865957.
  * [11] Valentín de la Rubia, Ulrich Razafison, and Yvon Maday.  Reliable fast frequency sweep for microwave devices via the reduced-basis method.  IEEE Trans. Microw. Theory Techn., 57(12):2923–2937, 2009.  doi:10.1109/TMTT.2009.2034208.
  * [12] Lihong Feng and Peter Benner.  A new error estimator for reduced-order modeling of linear parametric systems.  IEEE Trans. Microw. Theory Techn., 67(12):4848–4859, 2019.  doi:10.1109/TMTT.2019.2948858.
  * [13] Grzegorz Fotyga, Martyna Czarniewska, Adam Lamecki, and Michał Mrozowski.  Reliable greedy multipoint model-order reduction techniques for finite-element analysis.  IEEE Antennas Wireless Propag. Lett., 17(5):821–824, 2018.  doi:10.1109/LAWP.2018.2817391.
  * [14] Sara García, Valentín de la Rubia, and Michal Mrozowski.  Reduced basis approximations in microwave filters and diplexers: Inf-sup constant behavior.  In 2017 IEEE MTT-S International Conference on Numerical Electromagnetic and Multiphysics Modeling and Optimization for RF, Microwave, and Terahertz Applications (NEMO), pages 275–277, 2017.  doi:10.1109/NEMO.2017.7964258.
  * [15] Christophe Geuzaine and Jean-François Remacle.  Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities.  Internat. J. Numer. Methods Engrg., 79(11):1309–1331, 2009.  doi:10.1002/nme.2579.
  * [16] V. Girault and P.-A. Raviart.  Finite Element Methods for Navier-Stokes Equations.  Springer-Verlag, Berlin, Germany, 1986.  doi:10.5555/2408702.
  * [17] Yuxing He, Giuseppe Macchiarella, Gang Wang, Wentao Wu, Liguo Sun, Lu Wang, and Rong Zhang.  A direct matrix synthesis for in-line filters with transmission zeros generated by frequency-variant couplings.  IEEE Trans. Microw. Theory Techn., 66(4):1780–1789, 2018.  doi:10.1109/TMTT.2018.2791940.
  * [18] Martin W Hess and Peter Benner.  Fast evaluation of time–harmonic Maxwell’s equations using the reduced basis method.  IEEE Trans. Microw. Theory Techn., 61(6):2265–2274, 2013.  doi:10.1109/TMTT.2013.2258167.
  * [19] Martin W Hess, Sara Grundel, and Peter Benner.  Estimating the inf-sup constant in reduced basis methods for time-harmonic Maxwell’s equations.  IEEE Trans. Microw. Theory Techn., 63(11):3549–3557, 2015.  doi:10.1109/TMTT.2015.2473157.
  * [20] J. S. Hesthaven, G. Rozza, and B. Stamm.  Certified Reduced Basis Methods for Parametrized Partial Differential Equations.  SpringerBriefs in Mathematics. Springer, Cham, 2016.  doi:10.1007/978-3-319-22470-1.
  * [21] Amit Hochman, Jorge Fernández Villena, Athanasios G Polimeridis, Luís Miguel Silveira, Jacob K White, and Luca Daniel.  Reduced-order models for electromagnetic scattering problems.  IEEE Trans. Antennas Propag., 62(6):3150–3162, 2014.  doi:10.1109/TAP.2014.2314734.
  * [22] P Ingelstrom.  A new set of H(curl)-conforming hierarchical basis functions for tetrahedral meshes.  IEEE Trans. Microw. Theory Techn., 54(1):106–114, 2006.  doi:10.1109/TMTT.2005.860295.
  * [23] Mohamed Jemai and Ammar B Kouki.  New adaptive multi-expansion frequencies approach for SP-MORe techniques with application to the well-conditioned asymptotic waveform evaluation.  IEEE Trans. Microw. Theory Techn., 65(10):3709–3719, 2017.  doi:10.1109/TMTT.2017.2694819.
  * [24] A. Kirsch and F. Hettlich.  The Mathematical Theory of Time-Harmonic Maxwell’s Equations.  Springer, Cham, Switzerland, 2015.  doi:10.1007/978-3-319-11086-8.
  * [25] Yves Konkel, Ortwin Farle, Alexander Sommer, Stefan Burgard, and Romanus Dyczij-Edlinger.  A posteriori error bounds for Krylov-based fast frequency sweeps of finite-element systems.  IEEE Trans. Magn., 50(2):441–444, 2014.  doi:10.1109/TMAG.2013.2285442.
  * [26] K. Kurokawa.  An Introduction to the Theory of Microwave Circuits.  Academic, New York, NY, USA, 1969.
  * [27] Mohammad Memarian and Raafat R Mansour.  Quad-mode and dual-mode dielectric resonator filters.  IEEE Trans. Microw. Theory Techn., 57(12):3418–3426, 2009.  doi:10.1109/TMTT.2009.2034310.
  * [28] A. Monje-Real and V. de la Rubia.  Electric field integral equation fast frequency sweep for scattering of nonpenetrable objects via the reduced-basis method.  IEEE Trans. Antennas Propag., 68(8):6232–6244, 2020.  doi:10.1109/TAP.2020.2992882.
  * [29] P. Monk.  Finite Element Methods for Maxwell’s Equations.  Oxford University Press, New York, NY, USA, 2003.  doi:10.1093/acprof:oso/9780198508885.001.0001.
  * [30] Jean-Claude Nédélec.  Mixed finite elements in R3.  Numer. Math., 35(3):315–341, 1980.  doi:10.1007/BF01396415.
  * [31] Julio L Nicolini, Dong-Yeop Na, and Fernando L Teixeira.  Model order reduction of electromagnetic particle-in-cell kinetic plasma simulations via proper orthogonal decomposition.  IEEE Plasma Sci., 47(12):5239–5250, 2019.  doi:10.1109/TPS.2019.2950377.
  * [32] A. Quarteroni, A. Manzoni, and F. Negri.  Reduced Basis Methods for Partial Differential Equations, volume 92 of La Matematica per il 3+2.  Springer International Publishing, 2016.  ISBN: 978-3-319-15430-5.  doi:10.1007/978-3-319-15431-2.
  * [33] Michał Rewieński, Adam Lamecki, and Michał Mrozowski.  A goal-oriented error estimator for reduced basis method modeling of microwave devices.  IEEE Microw. Wireless Compon. Lett., 25(4):208–210, 2015.  doi:10.1109/LMWC.2015.2400937.
  * [34] Michal Rewienski, Adam Lamecki, and Michal Mrozowski.  Greedy multipoint model-order reduction technique for fast computation of scattering parameters of electromagnetic systems.  IEEE Trans. Microw. Theory Techn., 64(6):1681–1693, 2016.  doi:10.1109/TMTT.2016.2560167.
  * [35] Andreas Schmidt, Dominik Wittwar, and Bernard Haasdonk.  Rigorous and effective a-posteriori error bounds for nonlinear problems - application to RB methods.  Adv. Comput. Math., 46(2):1–30, 2020.  doi:10.1007/s10444-020-09741-x.
  * [36] Kathrin Smetana, Olivier Zahm, and Anthony T Patera.  Randomized residual-based error estimators for parametrized equations.  SIAM J. Sci. Comput., 41(2):A900–A926, 2019.  doi:10.1137/18M120364X.
  * [37] A. Sommer, O. Farle, and R. Dyczij-Edlinger.  Certified dual-corrected radiation patterns of phased antenna arrays by offline–online order reduction of finite-element models.  J. Comput. Phys., 299:22–44, 2015.  doi:10.1016/j.jcp.2015.06.024.
  * [38] A. Sommer, O. Farle, and R. Dyczij-Edlinger.  A new method for accurate and efficient residual computation in adaptive model-order reduction.  IEEE Trans. Magn., 51(3):1–4, 2015.  doi:10.1109/TMAG.2014.2352812.
  * [39] Lukasz Szydlowski, Adam Lamecki, and Michal Mrozowski.  Coupled-resonator filters with frequency-dependent couplings: Coupling matrix synthesis.  IEEE Microw. Wireless Compon. Lett., 22(6):312–314, 2012.  doi:10.1109/LMWC.2012.2197386.
  * [40] D. Szypulski, G. Fotyga, V. de la Rubia, and M. Mrozowski.  A subspace-splitting moment-matching model-order reduction technique for fast wideband fem simulations of microwave structures.  IEEE Trans. Microw. Theory Techn., 68(8):3229–3241, 2020.  doi:10.1109/TMTT.2020.3001627.
  * [41] Stefano Tamiazzo and Giuseppe Macchiarella.  Synthesis of cross-coupled filters with frequency-dependent couplings.  IEEE Trans. Microw. Theory Techn., 65(3):775–782, 2017.  doi:10.1109/TMTT.2016.2633258.
  * [42] W. Wang, G. N. Paraschos, and M. N. Vouvakis.  Fast frequency sweep of FEM models via the balanced truncation proper orthogonal decomposition.  IEEE Trans. Antennas Propag., 59(11):4142–4154, 2011.  doi:10.1109/TAP.2011.2164184.
  * [43] Li Xue and Dan Jiao.  Rapid modeling and simulation of integrated circuit layout in both frequency and time domains from the perspective of inverse.  IEEE Trans. Microw. Theory Techn., 68(4):1270–1283, 2020.  doi:10.1109/TMTT.2020.2966699.
  * [44] Ping Zhao and Ke-Li Wu.  An iterative and analytical approach to optimal synthesis of a multiplexer with a star-junction.  IEEE Trans. Microw. Theory Techn., 62(12):3362–3369, 2014.  doi:10.1109/TMTT.2014.2364222.


sec:experiments

Setup Implementation details. We implement the presented framework to train
and test on the proposed dataset. We employ Transformer encoder
vaswani2017attention for both the local and global state encoders, as well as
Transformer decoders for the policy networks. Each Transformer consists of $3$
layers with $8$ attention heads. We share parameters for policy networks
$\phi_{(1)}\dots\phi_{(K)}$. We set the strategic reasoning depth $K=3$ unless
otherwise stated. For all the methods, we provide $1$s motion history and
predict future $1$s motion. For additional experimental details, please refer
to the appendix.

Evaluation metrics. We compute the Mean Per Joint Position Error (MPJPE)
between the predicted future motion and ground-truth. In order to disentangle
the results of trajectory prediction and pose prediction, we further calculate
the mean root position error and mean local pose error (MPJPE after root
alignment). We report all the errors in millimeters.

Evaluation and comparison

In this section, we evaluate the proposed dataset and method, comparing with
the prior works. We conduct experiments with the following baseline methods
including a naive baseline and multiple state-of-the-art approaches: 1) Frozen
is a naive baseline that simply replicates the last frame of the input motion.
2) Social-STGCNN mohamed2020social is a multi-agent trajectory prediction
method using spatio-temporal graph convolutional networks. For trajectory
prediction methods, we train and test using the root (global) trajectories and
replicate the last frame for the root-relative (local) pose. 3) SocialVAE
socialvae2022 is a multi-agent trajectory prediction method using timewise
variational autoencoder. 4) History Repeats Itself (HRI) hri is a single-
person motion prediction method that allows absolute coordinate inputs. 5)
Multi-Range Transformers (MRT) mrt predicts multi-person motion utilizing a
local-range encoder and a global-range encoder. 6) Social Motion Transformer
(SoMoFormer) vendrow-somoformer employs the joints of all people as queries
and predicts their future motion in parallel.

[width=]figs/baseline.pdf

Figure : Comparison of baseline methods on CMU-Mocap [mocap] and our dataset.
We compute the $1$s mean local pose prediction errors. The baseline methods
are trained and tested on two datasets separately.

###  Quantitative results

We first benchmark three baseline methods to assess the motion prediction
difficulty of our dataset in comparison to the commonly-used existing dataset
[mocap]. Figure  reveals that for the CMU-Mocap dataset [mocap], merely
repeating the last-frame pose (_Frozen_) leads to significantly lower errors
than the model predictions. This finding indicates that the motions in this
dataset are relatively “static”, making it less challenging and informative
for the social motion prediction task. In contrast, the dynamic nature of our
dataset gives rise to considerably higher prediction errors of the naive
_Frozen_ baseline. The results demonstrate that the proposed dataset serves as
a more convincing benchmark for social motion prediction and introduces more
challenges to the field. Furthermore, we conduct a quantitative evaluation of
the proposed framework’s performance and compare it with state-of-the-art
approaches. As illustrated in Table , our method achieves competitive results
with SOTA approaches. It is worth noting that the single-person method HRI
[hri] demonstrates respectable short-term motion prediction accuracy; however,
its performance over longer horizons is constrained by the absence of global
context, particularly in root trajectory prediction. In contrast, our approach
significantly surpasses SOTA methods in long-term social motion prediction,
highlighting its effectiveness in modeling long-range human interactions. In
the subsequent section, we will demonstrate the critical role that social
context plays in predicting long-term human motion via examples.

Table : Performance comparison with the baseline methods. We compute the mean
prediction errors of global pose, local pose, and root, respectively.

! milliseconds Global Local Root 400 600 800 1000 400 600 800 1000 400 600 800
1000 Frozen 119.2 165.7 207.3 244.5 62.8 82.8 98.9 112.1 99.8 139.3 175.0
207.6 Social-STGCNN [mohamed2020social] 188.2 246.4 298.8 346.5 62.8 82.8 98.9
112.1 173.1 225.6 273.7 317.9 SocialVAE [socialvae2022] 84.0 127.2 171.7 215.1
62.8 82.8 98.9 112.1 52.6 89.7 130.3 170.8 HRI [hri] 50.0 106.8 145.5 189.2
39.5 62.9 78.4 94.1 40.1 88.0 119.6 157.7 MRT [mrt] 66.9 103.2 140.2 176.4
49.4 68.8 85.4 99.5 50.6 79.0 109.3 140.1 SoMoFormer [vendrow-somoformer] 53.2
88.8 124.9 160.0 42.3 62.6 79.7 93.9 42.5 70.7 100.8 131.2 Ours 54.6 86.2
119.3 152.5 43.7 60.8 74.6 86.6 41.7 66.9 94.8 124.0

[width=]figs/327_latest.pdf

Figure : Qualitative comparison with existing methods. Left two columns are
past motion and right three columns are future motion. We highlight the
prominent player with bright colors. Real motion are shown in
[HTML]6ebc85green, and model predictions are shown in [HTML]8e39c6purple.

[width=]figs/levels_combined_v7.pdf

Figure : Actions produced by policy networks at different levels. Left two
columns are past motions and right three columns are future motions.

###  Qualitative results

In order to better understand the performance improvement of our method, we
visualize its motion prediction results and compare with the baseline methods
qualitatively. Figure  illustrates an representative result in which the input
motion shows that a defensive player (highlighted) is switching the target of
defense. Our approach correctly predicts its future direction and orientations
as the offensive player at the other side will likely be covered by the
teammate. In contrast, other approaches fail to reason this, predicting that
the player would continue to move along the input trajectory or simply stop at
the middle. In Figure , three players to be noticed are highlighted in
[HTML]b471a6rose, [HTML]8e39c6purple and [HTML]278acbblue (we also call them
this way below). Additionally, we visualize the actions produced by the
intermediate policy networks to interpret the decision-making process. Shown
on the top row, player [HTML]8e39c6purple at level $1$ policy network simply
follows the past trend to receive the ball. As displayed in the second row,
defender [HTML]b471a6rose at level $2$ gains knowledge from the first level
that player [HTML]8e39c6purple is receiving the ball, and implements a steal
attempt. Nevertheless, such aggressive actions could create defensive gaps.
Defender [HTML]278acbblue at level $3$ senses such gaps caused by defender
[HTML]b471a6rose from level $2$, therefore executes cooperative defense
upwards to fill the gap.

##  Ablation study

Finally, we conduct ablation studies to investigate the influence of each
component. In Table  (a)-(e), we change the policy network depth $K$. $K=0$
corresponds to producing the final action directly without any recursive
reasoning. We find that introducing the recursive reasoning process ($K>0$)
clearly lowers the prediction errors, proving the effectiveness of our
framework design based on cognitive hierarchy. From (f) to (h), we keep $K=3$
and study the influence of training regularizations. (f) removes the GAIL
regularization for intermediate levels $1\dots K-1$, and (g) completely
removes GAIL. Training without GAIL increases the long-term trajectory
prediction errors, while slightly lowers the local pose prediction error. We
also observe that models without GAIL tend to produce motions that appear less
natural (e.g. limb twisting), leading to inferior visual effects. In (h), we
use different policy network instances for levels $1\dots K$, and it turns out
that sharing the policy networks leads to better generalization. In addition
to the qualitative visualization of actions generated by policy networks at
different reasoning levels, we also provide a quantitative comparison of their
distance from ground-truth. Table  demonstrates that policy networks at
different levels learn different actions, and policy networks at higher
reasoning levels produce actions with less prediction errors. We hypothesize
that, through the proposed recursive reasoning process, policy networks at
different levels display a range of cognitive types, with behavior spanning
from relatively naive to substantially rational. Higher-level policies more
closely approximate an equilibrium strategy, which is usually in alignment
with expert demonstrations. We provide several example motion sequences in the
video demo to further illustrate this observation.

Table : Ablation study of cognition levels and architecture designs.

! milliseconds Global Local Root 400 600 800 1000 400 600 800 1000 400 600 800
1000 (a) $K=0$ 60.1 94.3 129.2 163.7 47.4 66.2 82.0 95.3 45.7 72.7 101.8 131.8
(b) $K=1$ 55.2 88.0 121.9 155.5 44.3 62.1 77.0 89.7 42.3 68.0 96.3 125.4 (c)
$K=2$ 54.8 87.2 120.7 154.2 43.5 60.6 75.0 87.3 41.7 67.6 95.9 125.3 (d) $K=4$
54.5 86.3 119.7 153.5 43.6 60.5 74.7 86.6 42.0 67.3 95.5 125.2 (e) $K=3$, full
54.6 86.2 119.3 152.5 43.7 60.8 74.6 86.6 41.7 66.9 94.8 124.0 (f) _w/o_ mid
GAIL 54.9 87.5 121.4 155.4 43.7 60.4 74.5 86.4 41.8 67.8 96.6 126.9 (g) _w/o_
GAIL 54.6 87.1 121.1 155.1 43.4 60.4 74.7 86.8 41.9 68.1 97.0 127.0 (h) _w/o_
weight sharing 57.1 91.0 125.8 160.2 44.4 62.2 77.1 89.7 43.6 70.3 99.5 129.5

Table : Quantitative comparison of actions produced by policy networks at
different levels.

! milliseconds Global Local Root 400 600 800 1000 400 600 800 1000 400 600 800
1000 Level-$1$ 76.8 122.4 170.3 217.7 50.1 69.0 85.2 99.1 61.4 100.4 143.2
186.9 Level-$2$ 61.0 99.9 144.6 192.5 47.1 66.0 82.5 96.6 46.5 78.0 116.2
158.8 Level-$3$ 54.6 86.2 119.3 152.5 43.7 60.8 74.6 86.6 41.7 66.9 94.8 124.0