Departamento de Física Teórica and IPARCOS, Universidad Complutense de Madrid,
Plaza de las Ciencias 1, 28040 Madrid, Spain

# Thermal hadron resonances and Ward identities: results for the QCD phase
diagram

A. Gómez Nicolains:x  J. Ruiz de Elvirains:x  A. Vioque-Rodríguezins:x
ins:xins:x

###### Abstract

We review recent work regarding the role of light scalar resonances at finite
temperature for chiral symmetry and $U(1)_{A}$ restoration. The results
obtained are based on unitarized Chiral Perturbation Theory and Ward
Identities and are directly connected with presently open problems within the
QCD phase diagram.

##  1 Motivation

Hadronic matter under extreme conditions such as finite temperature and/or
chemical potentials in connection with the QCD phase diagram has been
thoroughly studied lately, through theoretical models based on effective
theories including Chiral Perturbation Theory (ChPT) [1, 2, 3, 4, 5, 6, 7, 8,
9, 10, 11, 12], lattice simulations [13, 14, 15, 16, 17, 18, 19, 20, 21, 22,
23] and experimental results from Heavy-Ion Collisions [24, 25]. Although most
of the relevant Physics within the QCD phase diagram has been understood,
there are still some open problems [12, 20, 21]. Among them, we will be
interested here in the relation between the Chiral Symmetry Restoration (CSR)
and $U(1)_{A}$ asymptotic restoration, which affects the very nature of the
transition and has many theoretical and phenomenological consequences
depending on how effective $U(1)_{A}$ restoration is at the QCD chiral
transition temperature $T_{c}\simeq$ 155 MeV [1, 9, 10, 12, 16, 17, 18, 19,
26, 27, 28, 29, 30, 31, 32, 33]. In addition, as discussed below, an important
role in this context is played by thermal resonances, i.e., resonances
generated within the heavy-ion or thermal environment, whose spectral
properties can be considerably modified due to the interactions among
particles in the thermal bath.

The main observables regarding CSR are the light quark condensate
$\left\langle{\bar{q}q}\right\rangle_{l}=\langle\bar{u}u+\bar{d}d\rangle\equiv\langle\bar{q}_{l}q_{l}\rangle$,
the order parameter of the chiral transition, and the light scalar
susceptibility $\chi_{S}$, defined as the two-point light $\bar{q}q$
correlator at vanishing four-momentum. On the one hand, the inflection point
of $\left\langle{\bar{q}q}\right\rangle_{l}$ at the crossover-like transition
temperature $T_{c}$ signals CSR for physical masses, while the condensate
vanishes at a lower $T_{c}$ in the light chiral limit $m_{u}=m_{d}=0$, where
the transition is of second order, and $\chi_{S}$ peaks for physical masses
and diverges in the chiral limit at $T_{c}$ [13, 14, 34]. On the other hand,
screening masses and susceptibilities in different channels provide also
relevant information, as they would become degenerate at chiral and/or
$U(1)_{A}$ restoration. Such chiral or $U(1)_{A}$ partners are particularly
useful regarding the effective theory description, as commented below.
Consider for instance the lightest meson states in the scalar/pseudoscalar
nonets for isospin $I=0,1/2,1$, whose corresponding quark bilinears degenerate
as

$\displaystyle\pi^{a}=i\bar{q}_{l}\gamma_{5}\tau^{a}q_{l}$
$\displaystyle\stackrel{{\scriptstyle SU(2)_{A}}}{{\longleftrightarrow}}$
$\displaystyle\sigma_{l}=\bar{q}_{l}q_{l},$
$\displaystyle\qquad\qquad{\big{\updownarrow}\mbox{\tiny$U(1)_{A}$}}$
$\displaystyle\qquad{\big{\updownarrow}\mbox{\tiny$U(1)_{A}$}}$
$\displaystyle\delta^{a}=\bar{q}_{l}\tau^{a}q_{l}$
$\displaystyle\stackrel{{\scriptstyle SU(2)_{A}}}{{\longleftrightarrow}}$
$\displaystyle\eta_{l}=i\bar{q}_{l}\gamma_{5}q_{l},$ $\displaystyle
K^{a}=i\bar{q}\gamma_{5}\lambda^{a}q$ $\displaystyle\stackrel{{\scriptstyle
SU(2)_{A},U(1)_{A}}}{{\longleftrightarrow}}$
$\displaystyle\kappa^{a}=\bar{q}\lambda^{a}q,$

where $\delta$ stands for the state corresponding to the $a_{0}(980)$ and
$\sigma_{l},\eta_{l}$ would correspond to the light-quark part of the
$f_{0}(500)/f_{0}(980)$ and $\eta/\eta^{\prime}$, respectively. The states
connected through a $SU(2)_{A}$ transformation are chiral partners, i.e.,
those whose associated screening masses and susceptibilities should become
degenerate for $T\geq T_{c}$, while those connected by an $U(1)_{A}$
transformation would degenerate asymptotically for large temperatures.

Lattice simulations yield the following degeneration pattern: for $N_{f}=2+1$
light flavors. i.e., including strangeness, chiral partners degenerate
reasonably around $T_{c}$ while $U(1)_{A}$ partner degeneration takes place at
a significantly higher temperature, typically around $1.3\,T_{c}$ from which
the system would reach a $O(4)\times U(1)_{A}$ restoration pattern [17, 22].
The case of the $K-\kappa$ sector deserves some additional comments: while
that pair degenerates both by $SU(2)_{A}$ and $U(1)_{A}$ transformations,
lattice results indicate that both symmetries need to be restored for
$K-\kappa$ degeneration. In the particular $N_{f}=2$ case, results are
qualitatively different, being compatible with a strong $U(1)_{A}$ restoration
at $T_{c}$ in the chiral limit [16, 18, 19, 22]. An important deal of
theoretical work in the last few years has been devoted to explaining the
differences in the behavior of these two cases, as summarized below.

##  2 Ward Identities and partner degeneration

Ward Identities derived formally from the QCD generating functional have
proven to be a very valuable tool to explore the degeneration of chiral and
$U(1)_{A}$ partners [7, 8, 9, 10, 12, 33, 35, 36]. Let us summarize here some
relevant results in those works. The following WI is particularly revealing in
that respect:

$\chi_{P}^{\eta_{l}\eta_{s}}(T)=-\frac{m_{l}}{2m_{s}}\left[\chi_{P}^{\pi\pi}(T)-\chi_{P}^{\eta_{l}\eta_{l}}(T)\right]=-\frac{2}{m_{l}m_{s}}\chi_{top}(T),$
(1)

where $\chi_{P}^{\phi_{a}\phi_{b}}(T)$ are pseudoscalar susceptibilities
corresponding to $\langle\phi_{a}\phi_{b}\rangle$ correlators,
$\eta_{s}=i\bar{s}\gamma_{5}s$ and $\chi_{top}$ is the topological
susceptibility corresponding to the $U(1)_{A}$ gluon anomaly operator. Now,
one can perform a $SU(2)_{A}$ transformation rotating the bilinears
$\eta_{l}\rightarrow\delta$ keeping $\eta_{s}$ invariant and since
$\langle\delta\eta_{s}\rangle$ is a scalar-pseudoscalar correlator, it
vanishes by parity. Therefore, $\chi_{P}^{\eta_{l}\eta_{s}}$ vanishes if
chiral symmetry is completely restored and, since
$\chi_{P}^{\pi\pi}-\chi_{P}^{\eta_{l}\eta_{l}}$ in (1) vanishes when both
chiral and $U(1)_{A}$ are restored, the above WI is consistent with $U(1)_{A}$
restoration at exact CSR. Such a result is compatible with the $N_{f}=2$
lattice results in the chiral limit since in that degenerate limit CSR is
exact at at $T_{c}$. Nevertheless, the effect of both strangeness and a
nonzero physical quark mass will generate a gap between chiral and $U(1)_{A}$
restoration, as observed in $N_{f}=2+1$ simulations. Such behavior has been
checked in a $U(3)$ ChPT calculation, showing in particular that the
transition temperatures of chiral and $U(1)_{A}$ degeneration of different
partners tend to coincide in the light chiral limit [10]. In addition, (1)
establishes that $\chi_{top}(T)$ behaves as an order parameter of $O(4)\times
U(1)_{A}$ restoration, as confirmed by theoretical analyses [33, 37, 38] and
lattice simulations [39]. Two additional important WIs are those pertaining to
the $I=1/2$ sector since they clarify the role of strangeness:

$\chi_{P}^{K}(T)=-\frac{\left\langle{\bar{q}q}\right\rangle_{l}(T)+2\langle\bar{s}s\rangle(T)}{m_{l}+m_{s}},\quad\chi_{S}^{\kappa}(T)=\frac{\left\langle{\bar{q}q}\right\rangle_{l}(T)-2\langle\bar{s}s\rangle(T)}{m_{s}-m_{l}}.$
(2)

Thus, taking into account that, due to CSR,
$\left\langle{\bar{q}q}\right\rangle_{l}(T)$ decreases with $T$ much more
abruptly than $\langle\bar{s}s\rangle(T)$, the identities (2) imply, on the
one hand, that $\chi_{P}^{K}(T)$ should decrease monotonically, while
$\chi_{S}^{\kappa}(T)$ peaks above $T_{c}$ [35]. The behavior of
$\chi_{S}^{\kappa}(T)$ below such peak is controlled by CSR, the slope
increasing with stronger CSR effects, like reducing the pion mass. Above the
peak, $U(1)_{A}$ restoration operates, driving $\chi_{S}^{\kappa}$ towards
degeneration with $\chi_{P}^{K}$. In Fig. 1, we plot the lattice
$\chi^{K,\kappa}$ susceptibilities “reconstructed” from the WI (2) from
condensate data, since there are no direct lattice determinations of
susceptibilities available in those channels. One can clearly observe the
behavior we have just commented on.

Figure 1: Susceptibilities in the $K-\kappa$ sector extracted from the WIs in
(2) and the quark condensate lattice values in [14, 15].

##  3 Thermal Resonances

The role of thermal effects in light hadron resonances has been recently
proved to be crucial to explain many properties of the QCD phase diagram, in
particular those related to CSR and the $U(1)_{A}$ symmetry that we are
interested in here. Thus, in [6, 7, 11] it has been shown that the thermal
$f_{0}(500)$ generated from the Inverse Amplitude Method (IAM) unitarization
of $\pi\pi$ scattering at finite temperature [40] dominates the scalar
susceptibility through the saturation relation

$\frac{\chi_{S}(T)}{\chi_{S}(0)}=\frac{M_{S}^{2}(0)}{M_{S}^{2}(T)},$ (3)

where $M_{S}^{2}(T)$ plays the role of the $f_{0}(500)$ thermal mass,
degenerating with the pion mass as $T$ approaches CSR [12], and is given by
$M_{S}^{2}=M_{p}^{2}-\Gamma_{p}^{2}/4$, being
$s_{p}=(M_{p}-i\Gamma_{p}/2)^{2}$ the pole of the resonance in the $s$-complex
plane in the second Riemann sheet across the thermal unitarity cut. It is
worth pointing out that at finite temperature, scattering of the incoming and
outgoing particles with thermal bath components modifies the $T=0$ unitarity
relation in two main aspects [7, 36, 40], namely, phase space factors are
corrected with Bose-Einstein thermal distribution functions and new Landau-
like pure thermal cuts appear for the case of unequal masses.

Figure 2: Unitarized scalar susceptibilties $\chi_{S}$ (left) and
$\chi_{S}^{\kappa}$ (right) within the IAM thermal resonance approach. In the
left panel, the lattice points come from [13], and $\chi_{S}(0)$ is chosen to
match the ChPT value in [41]. The bands correspond to the LECs uncertainties.

In Figure 2 we see that the unitarized saturated scalar susceptibility
reproduces the expected transition peak and is compatible with lattice results
below and around the peak within the uncertainty range of the Low Energy
Constants (LEC) of the ChPT NLO Lagrangian, unlike the purely perturbative
ChPT result [41], which is monotonically increasing and lies well above the
lattice points for temperatures approaching the transition. The same approach
has been followed for the thermal $K_{0}^{*}(700)$ or $\kappa$ resonance,
generated from $\pi K$ scattering at finite temperature [35, 36]. In this
case, the unitarized saturation approach reproduces the expected peak for
$\chi_{S}^{\kappa}$ explained in section 2. In Fig. 2, we show the result for
$\chi_{S}^{\kappa}$ within this approach, confirming the expected behavior
towards the light chiral limit, seen also in Fig. 1. In the $SU(3)$ limit
$m_{s}\rightarrow m_{l}$, one finds $\chi_{S}^{\kappa}\rightarrow\chi_{S}$ as
expected. Thermal resonances in the vector channel $\rho(770)$, $K^{*}(890)$
have also been studied, with results compatible with heavy-ion phenomenology
[6, 36].

##  4 Conclusions

Thermal resonances and Ward Identities play a crucial role in understanding
key properties of the QCD phase diagram, like those related to chiral and
$U(1)_{A}$ restoration reviewed here. The lightest scalar thermal resonances
saturate the corresponding susceptibilities in those channels, giving rise to
results compatible with lattice and WI when those resonances are generated
from unitarized meson-meson scattering.

###### Acknowledgements.

Work partially supported by research contracts PID2019-106080GB-C21 and
PID2022-136510NB-C31 (Spanish “Ministerio de Ciencia e Innovación”) and by the
European Union Horizon 2020 research and innovation program under grant
agreement No 824093. JRE acknowledges support from the Swiss National Science
Foundation, proj.No. PZ00P2_174228 and the Ramón y Cajal program
(RYC2019-027605-I) of the Spanish MINECO. A. V-R acknowledges support from a
fellowship of the UCM predoctoral program.

## References

  * [1] R. D. Pisarski and F. Wilczek, Phys. Rev. D, 29 (1984) 338.
  * [2] A. Bochkarev and J. I. Kapusta, Phys. Rev. D, 54 (1996) 4066.
  * [3] P. Gerber and H. Leutwyler, Nucl. Phys. B, 321 (1989) 387.
  * [4] A. Schenk, Phys. Rev. D, 47 (1993) 5138.
  * [5] R. Rapp and J. Wambach, Adv. Nucl. Phys., 25 (2000) 1.
  * [6] A. Dobado et al.,Phys. Rev. C, 66 (2002) 055201.
  * [7] A. Gómez Nicola, J. Ruiz de Elvira and R. Torres Andres, Phys. Rev. D, 88 (2013) 076007.
  * [8] A. Gómez Nicola and J. Ruiz de Elvira, JHEP, 03 (2016) 186.
  * [9] A. Gómez Nicola and J. Ruiz de Elvira, Phys. Rev. D, 97 (2018) 074016.
  * [10] A. Gómez Nicola and J. Ruiz de Elvira, Phys. Rev. D, 98 (2018) 014020.
  * [11] S. Ferreres-Solé et al., Phys. Rev. D, 99 (2019) 036018.
  * [12] A. Gómez Nicola, Eur. Phys. J. ST, 230 (2021) 1645.
  * [13] Y. Aoki et al., JHEP, 0906 (2009) 088.
  * [14] A. Bazavov et al. [HotQCD Collaboration], Phys. Rev. D, 85 (2012) 054503.
  * [15] A. Bazavov et al. [HotQCD Collaboration], Phys. Rev. D, 90 (2014) 094503.
  * [16] G. Cossu et al, Phys. Rev. D, 87 (2013) 114514.
  * [17] M. I. Buchoff et al., Phys. Rev. D, 89 (2014) 054514.
  * [18] B. B. Brandt et al, JHEP, 1612 (2016) 158.
  * [19] A. Tomiya et al, Phys. Rev. D, 96 (2017) 034509.
  * [20] C. Ratti, Rept. Prog. Phys., 81 (2018) 084301.
  * [21] A. Bazavov et al. [USQCD Collaboration], Eur. Phys. J. A, 55 (2019) 194.
  * [22] A. Bazavov, et al.. Phys. Rev. D, 100 (2019) 094510.
  * [23] S. Aoki et al. [JLQCD], Phys. Rev. D, 103 (2021) 074506.
  * [24] L. Adamczyk et al. [STAR Collaboration], Phys. Rev. C, 96 (2017) 044904.
  * [25] A. Andronic et al., Nature, 561 (2018) 321.
  * [26] M. Ishii, H. Kouno and M. Yahiro, Phys. Rev. D, 95 (2017) 114022.
  * [27] E. V. Shuryak, Comments Nucl. Part. Phys., 21 (1994) 235.
  * [28] J. I. Kapusta, D. Kharzeev and L. D. McLerran, Phys. Rev. D, 53 (1996) 5028.
  * [29] T. D. Cohen, Phys. Rev. D, 54 (1996) R1867.
  * [30] E. Meggiolaro and A. Morda, Phys. Rev. D, 88 (2013) 096010.
  * [31] A. Pelissetto and E. Vicari, Phys. Rev. D, 88 (2013) 105018.
  * [32] V. Dick, et al., Phys. Rev. D, 91 (2015) 094504.
  * [33] A. Gómez Nicola, J. Ruiz De Elvira and A. Vioque-Rodríguez, JHEP, 11 (2019) 086.
  * [34] H. T. Ding et al., Phys. Rev. Lett., 123 (2019) 062002.
  * [35] A. Gómez Nicola et al.,Eur. Phys. J. C, 81 (2021) 637.
  * [36] A. Gómez Nicola, J. Ruiz de Elvira and A. Vioque-Rodríguez, JHEP, 08 (2023) 148.
  * [37] V. Azcoiti, Phys. Rev. D, 94 (2016) 094505.
  * [38] G. Grilli di Cortona, E. Hardy, J. Pardo Vega and G. Villadoro, JHEP, 1601 (2016) 034.
  * [39] C. Bonati et al., JHEP, 1603 (2016) 155.
  * [40] A. Gómez Nicola, F. J. Llanes-Estrada and J. Pelaez, Phys. Lett. B, 550 (2002) 55.
  * [41] A. Gómez Nicola, J. R. Pelaez and J. Ruiz de Elvira, Phys. Rev. D, 87 (2013) 016001;


# Data-driven reduced-order modeling for parametric PDE eigenvalue problem
using Gaussian process regression

Abdul Halim  Daniele Boffi  F. Bertrand

###### Abstract

In this article, we propose a data-driven reduced basis (RB) method for
parametric eigenvalue problems. The method is based on offline and online
paradigm. In offline stage we generate snapshots and construct the RB basis
using POD approach. Gaussian process regressions (GPR) are used for
approximating the eigenvalues and projection coefficients of the eigenvectors
in the RB space. All GPR are trained in the offline stage using the data
generated in offline stage. The output corresponding to new values can be
obtained in online stage using the trained GPR. The model is used to solve
affine and non-affine parameter dependent eigenvalue problems. Numerical
results of the proposed non-intrusive model proves the robustness of the
model.

###### keywords:

Reduced basis method, Gaussian process regression, Eigenvalue problem, Proper
orthogonal decomposition, non-intrusive method.

##  1 Introduction

Eigenvalue problems occurs naturally in many areas of science and engineering.
It has lot of applications in fluid dynamics [BoffiCodina16, Akhan18],
structural mechanics [Bertrand19, Zheng16], electromagnetism [Boffi19],
nuclear reactor physics [Buchanetal13, reactorref1, GermanRagusa19], in
photonics [Photonic1, photonic2, Gilbert19], optics [optics01, optics04] for
optimal control and optimal design, population dynamics [Pdynamics], in
acoustic problem [acoustic13], solid state physics [Pau07a] for calculating
band structure. In nuclear physics, the eigenvalue problem represents the
energetic diffusion of the neutron where the first eigenvalue represents the
criticality of the reactor and the corresponding eigenvalue represent neutron
flux. For photonic crystal fibers, one has to solve two eigenvalue problem in
order to calculate the photonic band gap [Gilbert19, photonic2].

A lot of work has been done for solving eigenvalue problem using finite
element method (FEM) and adaptive FEM, after the pioneer works by Fix and
Strang [StrangFix]. There are few works for dealing eigenvalue problem in
mixed form [Boffi10, Boffietal13, Boffietal97, Boffietal00]. Many eigenvalue
problems are parameter dependent where the parameters may represent the
material properties, physical domain and initial/bounday conditions. In case
of parametric many-query context where we need to solve the problems for
different parameters will be too much expensive. In past few decades, reduced
order modeling (ROM) is used to reduced the dimensionality of the problem. The
idea is to reduce the high-fidelity (HF) problem into a problem with much
smaller dimension by projection or any other means. Among ROMs, RB methods are
efficient for solving parametric PDEs. RB basis are constructed in two ways
namely POD and greedy approach. In POD approach the HF is solved for pre-
selected values of parameters and then SVD is applied to the snapshot matrix
and in greedy approach the snapshots are the basis, where the parameters are
chosen based on some optimal criterion. In general, Galerkin projection is
used to calculate the RB solution in the reduced space.

The parametric eigenvalue problems are subject of the articles
[Machielsetal00, GermanRagusa19, Buchanetal13, ourrom22, Fumagallietal16,
Horgeretal17, Pau07a, Pau07b]. All of these works are proposed for solving the
affine parameter dependent problem using projection based Galerkin method. The
RB space is constructed either by POD or greedy approach. The works in the
initial days used POD based approach [Machielsetal00, GermanRagusa19,
Buchanetal13]. In [Machielsetal00], authors have proposed a RB method for
calculating the first eigenvalues and also established a posteriori estimate
for the first eigenvalues but they have used the POD method for construction
of basis. In [Buchanetal13], authors have introduce the fictitious time
parameter and used POD approach for finding first eigenpair of the parameter
dependent eigenvalue problem. Further, we have extended them for finding later
eigenpairs in [ourrom22] and investigated the behaviour of the eigenvalues
with different mesh. In [GermanRagusa19], a system of parametric eigenvalue
problem has been considered to find the first eigenpairs only.

In [Fumagallietal16], authors have derived a posteriori error estimate for the
first eigenvalues and eigenvectors used them for greedy approach. In
[Horgeretal17] authors have established a posteriori estimate for the repeated
eigenvalues and the estimate is used for greedy algorithm and they have
presented the numerical result for elasticity eigenvalue problem. Finding
multiple eigenvalues and eigenvectors simultaneously is the focus of the work
[Pau07a, Pau07b]. The authors of the [Pau07a, Pau07b] used projection based
Galerkin method and adopted the greedy approach for construction of the RB-
basis.

The affine dependence of the problem is the main key behind the success of the
RB methods which helps to decouple the problem into online-offline stage. But
this decoupling is not possible/difficult for non-affine and non-linear
problems. Although empirical interpolation method (EIM) and its discretized
version are developed to tackle these issues but these methods are not enough
as they are problem dependent and difficult to generalize. So, we want a model
which is problem independent. We will propose a data-driven model which will
depend only on input-output pairs.

In this current work, we solve the eigenvalue problem using data-driven
technique used in [mguo18, mguo19] for parametric PDEs, through Gaussian
process regression. In the offline phase we solve the high-fidelity (HF)
eigenvalue problem for some selected parameters. Then using POD technique we
construct the basis of the reduced space and with the help of the reduced
basis we will get reduced solution corresponding to the selected parameters.
We construct a Gaussian process regression model using the selected parameters
as input and the eigenvalues as outputs. If the dimension of the reduced space
is $N$ then we need to construct another $N$ GPR corresponding to each
coefficients of the reduced eigenvector. Our model is capable to calculate
multiple eigenvectors simultaneously also. In the online stage, for a new
parameter we get eigenvalue by evaluating the GPR model trained using
eigenvalues and the projection coefficients are obtained by evaluating other
GPRs. Multiplying the RB basis with the corresponding coefficients, we get an
approximate eigenvector of the HF problem. The paper is organised as follows:
in section 2, we have discussed Gaussian process regression and its training.
In section 3, we have discussed the POD approach, the generation of training
data proposed data-driven reduced order model(DD-ROM). Numerical results for
different eigenvalues are discussed in section 4.

##  2 Gaussian process regression

A Gaussian process (GP) is a collection of random variables, any finite number
of which follow joint Gaussian distribution. A Gaussian process is specified
by its mean function $m(x)$ and covariance function $\kappa(x,x^{\prime})$. We
write:

$f\sim GP(m,\kappa)$

and say that the function $f$ is distributed as a GP with mean function $m$
and covariance function $\kappa$.

Regression is a process of approximate an unknown output function from a set
of input-output data. In Gaussian process regression, we approximate the
unknown function by the mean of a Gaussian process. In GPR, we assume a priori
that the input-output pairs are generated from a Gaussian process (GP)
involving some parameters in the mean and variance function. Then using the
training data, we find the parameters of the mean and the variance of the GP
and get the posterior GP. The posterior GP is used to make predictions for new
inputs. The advantage of using the GPR is that we are able to (quantify) bound
the approximate function in terms of the variance. In our case, the
approximate function will be a function from the parameter space
$(\mathcal{P})$ to set of real numbers. The output of the functions will be
either the eigenvales or the projection co-efficients of the eigenvectors onto
the RB space. First, we will discuss the GPR of single valued parameter and
assume that there are no errors in the output data.

Let the output data be $y_{n}=U({\mu_{n}})$ for $n=1,2,\dots,n_{s}$. We want
to use these data to approximate the function $U$ by a Gaussian process
regression. Let $X=[U({\mu_{1}}),\dots,U({\mu_{n_{s}}})]$ and
$\boldsymbol{y}=\\{y_{1},y_{2},\dots,y_{n_{s}}\\}$. We want to predict the
value $y_{\ast}$ of $U$ at point ${\mu}_{\ast}$. Let $U$ be a Gaussian random
field with mean $m:\mathscr{P}\to\mathbb{R}$ and covariance function
$C:\mathscr{P}\times\mathscr{P}\to\mathbb{R}$, the random vector
$(\boldsymbol{y},y_{\ast})$ follow Gaussian distribution with mean
$(m(X),m({\mu}_{\ast}))$ and covariance matrix

$\begin{bmatrix}C(X,X)&C(X,{\mu}_{\ast})^{T}\\\
C(X,{\mu})&C({\mu}_{\ast},{\mu}_{\ast})\end{bmatrix}.$

Then [Sullivan15], the conditional distribution of $U({\mu}_{\ast})$ with
given observation $U(X)=\boldsymbol{y}$ follows Gaussian distribution with
mean

$m^{\star}({\mu}_{\ast})=m({\mu}_{\ast})+C({\mu}_{\ast},X)C(X,X)^{-1}(\boldsymbol{y}-m(X))$

and variance

$K^{\star}({\mu}_{\ast})=C({\mu}_{\ast},{\mu}_{\ast})-C({\mu}_{\ast},X)C(X,X)^{-1}C(X,{\mu}_{\ast}).$

For example, let $\mathcal{P}=[0,1]$ and $U$ follow a Gaussian process with
mean $m(x)=ax+b$ and variance

$C(x,x^{\prime}):=exp(-\frac{(x-x^{\prime})^{2}}{2l^{2}})$

where $l>0$ is the correlation length. Then we need to find the suitable
hyperparameter $\theta=(a,b,l)$ using the training data and by maximizing the
log-likelihood function.

###  2.1 Gaussian process regression with vector valued input and noisy output

In this subsection, we will consider the vector valued input parameters and
noisy output data. Let
$\mathcal{D}=\\{(\boldsymbol{x}_{i},y_{i}):i=1,2,\dots,n_{s}\\}$ be the
training set where $\boldsymbol{x}_{i}\in\mathscr{P}\subset\mathbb{R}^{d}$ and
$y_{i}$ are the corresponding outputs. Let the input-output pairs follow some
regression function $U:\mathscr{P}\to\mathbb{R}$, whose prior is defined as a
Gaussian process. Let the regression function be a GP corrupted by independent
Gaussian noise [mguo19], that is

$\displaystyle
y=U(\boldsymbol{x})+\epsilon,\quad\epsilon\sim\mathcal{N}(0,\sigma^{2}),\quad
U(\boldsymbol{x})\sim
GP(m(\boldsymbol{x}),\kappa(\boldsymbol{x},\boldsymbol{x}^{\prime}))$

where
$m(\boldsymbol{x})=\theta_{1}H_{1}(\boldsymbol{x})+\dots+\theta_{k}H_{k}(\boldsymbol{x})$
for some basis $\\{H_{1}(\boldsymbol{x}),\dots,H_{k}(\boldsymbol{x})\\}$
defined on $\mathscr{P}$. For our numerical result we used linear basis that
is $m(\boldsymbol{x})=\theta_{1}+\theta_{2}\boldsymbol{x}$. One choice of the
kernels is automatic relevance determination squared exponential

$\kappa(\boldsymbol{x},\boldsymbol{x}^{\prime})=\sigma_{1}^{2}exp\Big{(}-\frac{1}{2}\sum\limits_{i=1}^{d}\frac{(x_{i}-x_{i}^{\prime})^{2}}{l_{i}^{2}}\Big{)}$

where $l_{i}>0$ are the correlation length in the coordinate diections. Given
a finite number of points in input domain, the outputs follow a joint
Gaussian/Normal distribution that is

$\boldsymbol{y}|\boldsymbol{X}\sim\mathcal{N}(m(\boldsymbol{X}),\boldsymbol{K}_{y}),\quad\boldsymbol{K}_{y}=cov(\boldsymbol{y}|\boldsymbol{X})=\kappa(\boldsymbol{X},\boldsymbol{X})+\sigma^{2}\mathbb{I}_{n_{s}}$
(2.1)

where
$\boldsymbol{y}=\\{y_{1},y_{2},\dots,y_{n_{s}}\\}^{T},\boldsymbol{X}=[\boldsymbol{x}_{1}|\boldsymbol{x}_{1}|\dots|\boldsymbol{x}_{N_{s}}]$
and $\mathbb{I}_{n_{s}}$ is the unit matrix of dimension $n_{s}$. Let
$\boldsymbol{y}^{\ast}$ be a set of functional values corresponding to the
test set inputs $\boldsymbol{X}^{\ast}$ then the goal is to predict the output
$\boldsymbol{y}^{\ast}(\boldsymbol{X}^{\ast})$ for a new test input set
$\boldsymbol{X}^{\ast}$. From the properties of Gaussian distribution the
random vector $(\boldsymbol{y},\boldsymbol{y}^{\ast})^{T}$ follow the
following joint normal distribution

$\begin{bmatrix}\boldsymbol{y}\\\
\boldsymbol{y}^{\ast}\end{bmatrix}\sim\mathcal{N}\Bigg{(}\begin{bmatrix}m(\boldsymbol{X})\\\
m(\boldsymbol{X}^{\ast})\end{bmatrix},\begin{bmatrix}\boldsymbol{K}_{y}&\kappa(\boldsymbol{X},\boldsymbol{X}^{\ast})\\\
\kappa(\boldsymbol{X}^{\ast},\boldsymbol{X})&\kappa(\boldsymbol{X}^{\ast},\boldsymbol{X}^{\ast})\\\
\end{bmatrix}\Bigg{)}.$ (2.2)

Since we know the values of $\boldsymbol{y}$ at the training set, we are
interested in the conditional distribution of $\boldsymbol{y}^{\ast}$ given
$\boldsymbol{y}$ which is expressed as:

$\boldsymbol{y}^{\ast}|\boldsymbol{y}\sim\mathcal{N}\bigg{(}m(\boldsymbol{X}^{\ast})+\kappa(\boldsymbol{X}^{\ast},\boldsymbol{X})\boldsymbol{K}_{y}^{-1}(\boldsymbol{y}-m(\boldsymbol{X})),\kappa(\boldsymbol{X}^{\ast},\boldsymbol{X}^{\ast})-\kappa(\boldsymbol{X}^{\ast},\boldsymbol{X})K_{y}^{-1}\kappa(\boldsymbol{X},\boldsymbol{X}^{\ast})\bigg{)}$

This is the posterior distribution for a specific set of test cases. Thus the
posterior function follow a new GP [Gprpaper]:

$U^{\ast}(\boldsymbol{\mu})\sim
GP(m^{\ast}(\boldsymbol{\mu})),\kappa^{\ast}(\boldsymbol{\mu},\boldsymbol{\mu}^{\prime})),$

where

$m^{\ast}(\boldsymbol{\mu})=m(\boldsymbol{\mu})+\kappa(\boldsymbol{\mu},\boldsymbol{X})\boldsymbol{K}_{y}^{-1}(\boldsymbol{y}-m(\boldsymbol{X})),\quad\kappa^{\ast}(\boldsymbol{\mu},\boldsymbol{\mu}^{\prime})=\kappa(\boldsymbol{\mu},\boldsymbol{\mu}^{\prime})-\kappa(\boldsymbol{\mu},\boldsymbol{X})K_{y}^{-1}\kappa(\boldsymbol{X},\boldsymbol{\mu}^{\prime}).$

Note that the posterior variance is smaller than the prior variance which is
achieved using the data.

The hyperparameters
$\boldsymbol{\theta}=(\theta_{1},\dots,\theta_{k},l_{1},\dots,l_{d},\sigma,\sigma_{1})$
plays an important role for the prediction. As the random vector
$\boldsymbol{y}|\boldsymbol{X}$ follow the mutivariate Gaussian distribution,
so form (2.1) the likelihood function is given by:

$\displaystyle
p(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\theta})=\frac{1}{\sqrt{(2\pi)^{N_{s}}|\boldsymbol{K}_{y}(\boldsymbol{\theta})|}}exp\Big{(}-\frac{1}{2}(\boldsymbol{y}-m(\boldsymbol{X}))^{T}\boldsymbol{K}_{y}^{-1}(\boldsymbol{\theta})(\boldsymbol{y}-m(\boldsymbol{X}))\Big{)}$

where $|\boldsymbol{K}_{y}(\boldsymbol{\theta})|$ denotes the determinant of
the matrix $\boldsymbol{K}_{y}$. The optimal hyperparameters is obtained my
maximizing the log likelihood function that is

$\displaystyle\boldsymbol{\theta}_{opt}$
$\displaystyle=arg\max\limits_{\boldsymbol{\theta}}logp(\boldsymbol{y}|\boldsymbol{X},\boldsymbol{\theta})$
$\displaystyle=arg\max\limits_{\boldsymbol{\theta}}\Big{\\{}-\frac{1}{2}(\boldsymbol{y}-m(\boldsymbol{X}))^{T}K_{y}^{-1}(\boldsymbol{y}-m(\boldsymbol{X}))-\frac{1}{2}log|K_{y}(\boldsymbol{\theta})|-\frac{N_{s}}{2}log(2\pi)\Big{\\}}$

For more details on GPR we refer the readers to [Gprbk, Gprpaper, Sullivan15,
mguo19].

##  3 Problem settings and data-driven reduced order model

Let us consider the parametric PDE eigenvalue problem in variational form:
given $\mu\in\mathscr{P}$ find $u(\mu)\in V$ and $\lambda(\mu)\in\mathbb{R}$
such that

$a(u(\mu),v;\mu)=\lambda(\mu)b(u(\mu),v;\mu)\quad\forall v\in V$ (3.1)

where $a,b:V\times V\to\mathbb{R}$ are bilinear forms and $V$ is a suitable
Hilbert space. Let $V_{h}\subset V$ be a finite dimensional finite element
space and $dim(V_{h})=N_{h}$. The discrete problem reads: for any given
$\mu\in\mathscr{P}$ find $u_{h}(\mu)\in V_{h}$ and
$\lambda_{h}(\mu)\in\mathbb{R}$ such that

$a(u_{h}(\mu),v_{h};\mu)=\lambda_{h}(\mu)b(u_{h}(\mu),v_{h};\mu)\quad\forall
v_{h}\in V_{h}$ (3.2)

Let $\\{\phi^{j}\\}_{j=1}^{N_{h}}$ denote the basis of $V_{h}$ and
$u_{h}(\mu)=\sum\limits_{j=1}^{N_{h}}u^{j}_{h}(\mu)\phi^{j}$ then the problem
is equivalent to the solution of a linear system

$\mathbb{A}_{h}(\mu)\boldsymbol{u}_{h}(\mu)=\lambda_{h}\mathbb{B}_{h}(\mu)\boldsymbol{u}_{h}(\mu)$
(3.3)

where
$\boldsymbol{u}_{h}(\mu)=(u^{1}_{h}(\mu),\dots,u^{N_{h}}_{h}(\mu))^{T}\in\mathbb{R}^{N_{h}}$
is the vector containing the nodal values and

$(\mathbb{A}_{h}(\mu))_{i,j}=a(\phi^{j},\phi^{i};\mu),\quad(\mathbb{B}_{h}(\mu))_{i,j}=b(\phi^{j},\phi^{i};\mu),\quad
1\leq i,j\leq N_{h}$.

###  3.1 Projection based reduced order modeling

Let $V_{N}\subset V_{h}$ be the reduced space of dimension $N$ and
$V_{N}=span\\{\zeta_{1},\dots,\zeta_{N}\\}$. There are two methods to obtain
reduced basis namely POD and greedy approach. In this paper, we use POD
approach which is discussed in details later in subsection 3.2. The reduced
order problem reads [Quarteronietal16]: for any given $\mu\in\mathscr{P}$ find
$u_{N}(\mu)\in V_{N}$ and $\lambda_{N}(\mu)\in\mathbb{R}$ such that

$a(u_{N}(\mu),v_{N};\mu)=\lambda_{N}(\mu)b(u_{N}(\mu),v_{N};\mu)\quad\forall
v_{N}\in V_{N}.$ (3.4)

Let $u_{N}(\mu)=\sum\limits_{j=1}^{N}u^{j}_{N}(\mu)\zeta_{j}$. Putting
$v_{N}=\zeta_{n}$ we get linear system is for the reduced model as

$\mathbb{A}_{N}({\mu})\boldsymbol{u}_{N}({\mu})=\lambda_{N}(\mu)\mathbb{B}_{N}({\mu})\boldsymbol{u}_{N}({\mu})$
(3.5)

where
$\boldsymbol{u}_{N}(\mu)=(u^{1}_{N}(\mu),\dots,u^{N}_{N}(\mu))^{T}\in\mathbb{R}^{N}$
is the vector containing projection co-efficients of the eigenvector
$u_{h}(\mu)$ and

$(\mathbb{A}_{N}({\mu}))_{n,m}=a(\zeta_{m},\zeta_{n};{\mu}),\quad(\boldsymbol{f}_{N}({\mu}))_{n}=f(\zeta_{n};{\mu}),\quad
1\leq m,n\leq N.$

If we assume the affine parametric dependence of the bilinear forms then we
can write $\mathbb{A}_{N}({\mu})$ in terms of $\mathbb{A}_{h}({\mu})$,
otherwise we can not write it. So for non-affine and non-linear problem data-
driven model can be helpful. Let the bilinear forms are affinely dependent
that is it can be written as follows:

$\displaystyle a(w,v;{\mu})$
$\displaystyle=\sum\limits_{q=1}^{Q_{a}}\theta_{a}^{q}({\mu})a_{q}(w,v),\quad\forall
v,w\in V,{\mu}\in\mathscr{P},$ (3.6) $\displaystyle b(w,v;{\mu})$
$\displaystyle=\sum\limits_{q=1}^{Q_{b}}\theta_{b}^{q}({\mu})b_{q}(w,v),\quad\forall
v,w\in V,{\mu}\in\mathscr{P}.$ (3.7)

where the bilinear forms $a_{q}(\cdot,\cdot),\quad q=1,2,\dots,Q_{a}$ and
$b_{q^{\prime}}(\cdot),\quad q^{\prime}=1,2,\dots,Q_{b}$ are parameter
independent. That is we are able to write the bilinear forms as a sum of
products of few parameter dependent functions $\theta_{q}$ and parameter
independent bilinear forms $a_{q}$. Then the matrix $\mathbb{A}_{h}$ and the
vector $\mathbb{B}_{h}$ of the high fidelity problem can be written as

$\mathbb{A}_{h}({\mu})=\sum\limits_{q=1}^{Q_{a}}\theta_{a}^{q}({\mu})\mathbb{A}_{h}^{q},\quad\mathbb{B}_{h}({\mu})=\sum\limits_{q=1}^{Q_{b}}\theta_{b}^{q}({\mu})\mathbb{B}_{h}^{q},$

and the matrix $\mathbb{A}_{N}$ and $\mathbb{B}_{N}$ of the RB problem can be
written as

$\mathbb{A}_{N}({\mu})=\sum\limits_{q=1}^{Q_{a}}\theta_{a}^{q}({\mu})\mathbb{A}_{N}^{q},\quad\mathbb{B}_{N}({\mu})=\sum\limits_{q=1}^{Q_{b}}\theta_{b}^{q}({\mu})\mathbb{B}_{N}^{q},$

Let $\mathbb{V}\in\mathbb{R}^{N_{h}\times N}$ be the transformation matrix of
the basis of $V_{h}$ and $V_{N}$. More precisely, let us write the basis of
$V_{N}$ in terms of the basis of $V_{h}$ as

$\zeta_{m}=\sum\limits_{i=1}^{N_{h}}\zeta_{m}^{(i)}\phi^{i},\quad 1\leq m\leq
N,$ (3.8)

Define the transformation matrix $\mathbb{V}\in\mathbb{R}^{N_{h}\times N}$ as
$\mathbb{V}=[\boldsymbol{\zeta}_{1}|\dots|\boldsymbol{\zeta}_{N}]$ where
$\boldsymbol{\zeta}_{m}=(\zeta_{m}^{(1)},\zeta_{m}^{(2)},\dots,\zeta_{m}^{(N_{h})})^{T}$.
Thus the matrix $\mathbb{A}_{N}$ and $\mathbb{B}_{N}$ of the RB problem can be
written as [Quarteronietal16]

$\mathbb{A}_{N}({\mu})=\sum\limits_{q=1}^{Q_{a}}\theta_{a}^{q}({\mu})\mathbb{V}^{T}\mathbb{A}_{h}^{q}\mathbb{V},\quad\mathbb{B}_{N}({\mu})=\sum\limits_{q=1}^{Q_{b}}\theta_{b}^{q}({\mu})\mathbb{V}^{T}\mathbb{B}_{h}^{q}\mathbb{V}$
(3.9)

In the offline stage we calculate and store the matrices
$\mathbb{V}^{T}\mathbb{A}_{h}^{q}\mathbb{V}$ and
$\mathbb{V}^{T}\mathbb{B}_{h}^{q}\mathbb{V}$. In the online stage, for a new
parameter ${\mu}$ we only evaluate the functions $\theta_{a}^{q}$ and
$\theta_{b}^{q}$ and then with the help of stored matrices we will get the
linear system (3.9) for the ROM. But if the bilinear/linear form are not
affine that is if we can not separate the parameter and space dependent
function then we can not use the matrices HF problem the to get the system of
RB method, as they will be parameter dependent.  That’s why we propose a data-
driven model which is discussed in section 3.3.

###  3.2 Proper orthogonal decomposition (POD)

Let $\mathbb{S}$ be the matrix of the $n_{s}$ snapshots, containing the
eigenvectors at the sample points $\mu_{1},\dots,\mu_{n_{s}}$. The sample
points can be selected using uniform tensorial sampling or latin hypercube
samplinr or Smolyak sparse sampling. The SVD gives

$\mathbb{S}=\mathbb{U}\Sigma\mathbb{Z}^{\top}$

with

$\mathbb{U}=[\boldsymbol{\zeta}_{1},\dots,\boldsymbol{\zeta}_{N_{h}}]\in\mathbb{R}^{N_{h}\times
N_{h}}\quad\text{and}\quad\mathbb{Z}=[\boldsymbol{\psi}_{1},\cdots,\boldsymbol{\psi}_{n_{s}}]\in\mathbb{R}^{n_{s}\times
n_{s}}$

orthogonal matrices and
$\Sigma=\mathrm{diag}(\sigma_{1},\dots,\sigma_{r},0,\dots,0)\in\mathbb{R}^{N_{h}\times
n_{s}}$ with $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{r}>0$, where $r$
is the rank of the matrix $\mathbb{S}$. The first $N$ left eigenvecors will be
the basis i.e. $\\{\boldsymbol{\zeta}_{1},\dots,\boldsymbol{\zeta}_{N}\\}$
will be the basis of the low-dimensional subspace $V_{N}$. We will not use the
SVD directly to find the left singular vectors unless $n_{s}\geq N_{h}$. Next,
we describe the process of calculating the POD basis. We can write

$\displaystyle\mathbb{S}\boldsymbol{\psi}_{i}=\sigma_{i}\boldsymbol{\zeta}_{i}\quad\text{and}\quad\mathbb{S}^{\top}\boldsymbol{\zeta}_{i}=\sigma_{i}\boldsymbol{\psi}_{i},\quad
i=1,\dots,r$ (3.10)

or, equivalently,

$\displaystyle\mathbb{S}^{\top}\mathbb{S}\boldsymbol{\psi}_{i}=\sigma^{2}_{i}\boldsymbol{\psi}_{i}\quad\text{and}\quad\mathbb{S}\mathbb{S}^{\top}\boldsymbol{\zeta}_{i}=\sigma^{2}_{i}\boldsymbol{\zeta}_{i},\quad
i=1,\dots,r.$ (3.11)

The matrix $\mathbb{C}=\mathbb{S}^{T}\mathbb{S}$ is called correlation matrix.
The POD basis $\mathbb{V}$ of dimension $N\leq n_{s}$ is defined as the set of
the first $N$ left singular vectors
$\\{\boldsymbol{\zeta}_{1},\dots,\boldsymbol{\zeta}_{N}\\}$ of $\mathbb{S}$
or, equivalently, the set of vectors

$\displaystyle\boldsymbol{\zeta}_{j}=\frac{1}{\sigma_{j}}\mathbb{S}\boldsymbol{\psi}_{j},\qquad
1\leq j\leq N$ (3.12)

obtained from the first $N$ eigenvectors
$\\{\boldsymbol{\psi}_{1},\dots,\boldsymbol{\psi}_{N}\\}$ of the correlation
matrix $\mathbb{C}$. So, instead of using the SVD of the matrix $\mathbb{S}$
we solve the eigenvalue problem
$\mathbb{C}\boldsymbol{\psi}=\sigma^{2}\boldsymbol{\psi}$ which is of size
$n_{s}$; notice that $n_{s}$ is usually much smaller than $N_{h}$.

The POD dimension $N$ is the smallest positive integer such that

$I(N)=\frac{\sum\limits_{i=1}^{N}\sigma_{i}^{2}}{\sum\limits_{i=1}^{r}\sigma_{i}^{2}}\geq
1-\varepsilon.$ (3.13)

where $\varepsilon$ is the tolerance. For the numerical computation of this
paper we have used $\varepsilon=10^{-8}$.

###  3.3 Data-driven model

In this subsection, we describe the data-driven model for calculating first
eigenpair which can be extended to other eigenpair also, by considering the
eigenvector of interest in the snapshot matrix. Next, subsection is dedicated
to deal with finding multiple eigenvectors simultaneously.

First, select the sample parameters $\mu_{1},\mu_{2},\dots,\mu_{n_{s}}$ using
sampling techniques like tensorial sampling, Monte Carlo sampling, Latin
hypercube sampling, Smolyak sparse sampling with Clenshaw-Curits points from
the parameter space. We calculate the first eigenpairs
$\big{(}\boldsymbol{u}_{h}(\mu_{1}),\lambda_{1}(\mu_{1})\big{)},\big{(}\boldsymbol{u}_{h}(\mu_{2}),\lambda_{1}(\mu_{2})\big{)},\cdots,\big{(}\boldsymbol{u}_{h}(\mu_{n_{s}}),\lambda_{1}(\mu_{n_{s}})\big{)}$
of high fidelity eigenvalue problem (3.2) in the finite dimensional space
$V_{h}$ of dimension $N_{h}$. These eigenvectors are the snapshots
corresponding to selected parameters $\mu_{1},\mu_{2},\dots,\mu_{n_{s}}$.
Consider the snapshot matrix
$\mathbb{S}=[\boldsymbol{u}_{h}(\mu_{1})|\cdots|\boldsymbol{u}_{h}(\mu_{n_{s}})]$,
containing the high-fidelity (HF) solutions as columns. Applying the SVD or
using POD approach to the snapshot matrix $\mathbb{S}$, we construct the
transform matrix $\mathbb{V}$ by taking first few dominating left singular
vectors of $\mathbb{S}$. The columns of $\mathbb{V}$ are the POD basis. If
$\boldsymbol{u}_{h}$ is the vector containing the nodal values of the solution
$u_{h}\in V_{h}$ and $\boldsymbol{u}_{N}$ is the vector containing the
projection co-efficients of the corresponding reduced solution $u_{N}\in
V_{N}$. Then the HF solutions and the RB solutions are related as:

$\textbf{u}_{N}=\mathbb{V}^{T}\textbf{u}_{h},\quad\boldsymbol{u}_{h}=\mathbb{V}\boldsymbol{u}_{N}.$

Thus the columns of
$\mathbb{S}_{N}:=\mathbb{V}^{T}\mathbb{S}=[\mathbb{V}^{T}\boldsymbol{u}_{h}(\mu_{1})|\cdots|\mathbb{V}^{T}\boldsymbol{u}_{h}(\mu_{n_{s}})]\approx[\boldsymbol{u}_{N}(\mu_{1})|\cdots|\boldsymbol{u}_{N}(\mu_{n_{s}})]\in\mathbb{R}^{N\times
n_{s}}$ are solutions of the reduced problem at
$\mu_{1},\mu_{2},\dots,\mu_{n_{s}}$. Hence the $i^{th}$ row of the matrix
$\mathbb{S}_{N}$ contains the $i^{th}$ coefficients of the reduced solutions
at $\mu_{1},\mu_{2},\dots,\mu_{n_{s}}$. So, $i^{th}$ row of $\mathbb{S}_{N}$
serves as the outputs of the $i^{th}$ coefficient at the input points
$\mu_{1},\mu_{2},\dots,\mu_{n_{s}}$.

Let $P_{\lambda}:\mathcal{P}\to\mathbb{R}$ is the mapping defined by
$\mu\mapsto\lambda(\mu)$, we want to approximate the map $P_{\lambda}$ by a
non-linear GP-based regression $\hat{P}_{\lambda}$ using the training data
$\\{\big{(}\mu_{i},\lambda_{1}(\mu_{i})\big{)}:i=1,2,\dots,n_{s}\\}$.

Similarly, let $P_{N}:\mathcal{P}\to\mathbb{R}^{N}$ is the mapping defined by
$\mu\mapsto\boldsymbol{u}_{N}(\mu)$, we want to approximate the map $P_{N}$ by
a non-linear GP-based regression $\hat{P}_{N}$ using the training data
$\\{\big{(}\mu_{i},\mathbb{V}^{T}\boldsymbol{u}_{h}(\mu_{i})\big{)}:i=1,2,\dots,n_{s}\\}$.
The GP-based regression models are constructed during the offline stage. In
the online stage the eigenvalues at a new parameter $\mu_{\ast}$ is obtained
by evaluating the regression $\hat{P}_{\lambda}$ at that point and the
corresponding eigenvector is obtained as
$u_{N}(\mu_{\ast})=\sum\limits_{j=1}^{N}\hat{P}^{j}_{N}(\mu_{\ast})\zeta_{j}$
where we assume that $\hat{P}_{N}=(\hat{P}^{1}_{N},\dots,\hat{P}^{N}_{N})$
with $\hat{P}^{j}_{N}$ is the regression corresponding to the mapping
$P^{j}_{N}:\mathcal{P}\to\mathbb{R}$ defined by $\mu\mapsto u^{j}_{N}(\mu)$
and the corresponding training data :

$\displaystyle\\{\big{(}\mu_{i},{u}_{N}^{j}(\mu_{i})\big{)}:i=1,2,\dots,n_{s},j=1,2,\dots,N$

where we use the notations of (3.5).

Let
$\boldsymbol{u}_{N}(\mu_{\ast})=(\hat{P}^{1}_{N}(\mu_{\ast}),\dots,\hat{P}^{N}_{N}(\mu_{\ast}))^{T}$
be the vector containing the projection co-efficients obtained using the GPR,
then pre-multiplying it by $\mathbb{V}$, we get the corresponding approximate
eigenvector $\boldsymbol{u}_{h}(\mu_{\ast})$ of the EVP (3.2). In a similar
way considering other eigenpair we can construct the corresponding GPR in
offline phase and evaluating those GPR in online phase we can get the other
eigenpair for a new parameter. For numerical implementation of GPR we have
used MATLAB command $fitrgp$.

###  3.4 Data-driven model for finding multiple eigenvectors simultaneously

Let, we are interested to find $n_{e}$ number of eigenvalues and eigenvectors
for a given EVP. Let $\mu_{1},\dots,\mu_{n_{s}}$ are sample points, then we
need to find all the $n_{e}$ eigenpairs at these sample points. For each
eigenvalues we need to train one GPR using the training data
$\\{\big{(}\mu_{i},\lambda_{j}(\mu_{i})\big{)}:i=1,2,\dots,n_{s}\\}$, for
$j=1,2,\dots,n_{e}$. But we can find the eigenvectors simultaneously by
considering the following snapshot matrix with

$\mathbb{S}_{s}=[\boldsymbol{s}_{1}|\dots|\boldsymbol{s}_{n_{s}}]\in\mathbb{R}^{n_{e}N_{h}\times
n_{s}}$

where

$\displaystyle\boldsymbol{s}_{j}$
$\displaystyle=\begin{bmatrix}\boldsymbol{u}_{1,h}(\mu_{j})\\\
\boldsymbol{u}_{2,h}(\mu_{j})\\\ \vdots\\\
\boldsymbol{u}_{n_{e},h}(\mu_{j})\end{bmatrix}$ (3.14)

Then we proceed as earlier with the modified snapshot matrix, we will get a
transformation matrix $\mathbb{V}_{s}\in\mathbb{R}^{n_{e}N_{h}\times N}$,
where $N$ is the number of POD basis corresponding to the new snapshot matrix,
obtained using the criterion (3.13). Then we get the matrix
$\mathbb{S}_{s,N}=\mathbb{V}_{s}^{T}\mathbb{S}_{s}\in\mathbb{R}^{N\times
n_{s}}$, each row will serve as the training data. Let
$P^{s}_{N}:\mathcal{P}\to\mathbb{R}^{N}$ is the mapping defined that maps the
$\mu$ to reduced co-efficeints. We want to approximate the map $P^{s}_{N}$ by
a non-linear GP-based regression $\hat{P}^{s}_{N}$ using the training data
$\\{\big{(}\mu_{i},\mathbb{V}_{s}^{T}\boldsymbol{s}(\mu_{i})\big{)}:i=1,2,\dots,n_{s}\\}$.
The GP-based regression models are constructed during the offline stage. In
the online stage all the GPR $\hat{P}^{s}_{N}$ are evaluated at parameter
$\mu_{\ast}$, collection all the coefficients form a vector
${\boldsymbol{\tilde{u}}_{N}}=(\hat{P}^{s}_{1}(\mu_{\ast}),\dots,\hat{P}^{s}_{N}(\mu_{\ast}))^{T}.$
Then the vector $\mathbb{V}_{s}{\boldsymbol{\tilde{u}}_{N}}$ contains first
$n_{e}$ eigenvectors at $\mu_{\ast}$. We can extract all the eigenvectors from
this vector.

The advantage of this method is that we need to train only $N$ GPR. So for
finding $n_{e}$ number of eigenpair, we need to training $n_{e}+N$ GPR but if
we do it separately we need $n_{e}(N+1)$ GPR.

##  4 Numerical result

First, we test our data-driven model on some affine parametric dependent
eigenvalue problems and then apply it to nonlinear and non-affine parameter
dependent eigenvalue problems. Finally, we test our model on the eigenvalue
problem with crossing eigenvalues.

###  4.1 Results of affine parameter dependent eigenvalue problem

In this subsection, we considered simple harmonic oscillator eigenvalue
problem both in 1-D and 2-D. The analytical solutions of these problems are
known. Note that the bilinear forms for these problems can be written affinely
and we can use RB methods. Here we will use data driven model to verify that
the data-driven model is working well for this case.

####  4.1.1 Simple harmonic oscillator eigenvalue problem in 1D

Let us consider the eigenvalue problem defined on
$\Omega=(-L,L)\subset\mathbb{R}$

$\displaystyle-\frac{1}{2}\frac{d^{2}u}{dx^{2}}+\frac{1}{2}\mu^{2}x^{2}u=\lambda
u$ (4.1)

with initial conditions $u(-L)=u(L)=0$.

The eigenvalues eigenvectors of this problem are

$\displaystyle\lambda_{n}(\mu)$ $\displaystyle=(n+\frac{1}{2})\mu\quad
n=0,1,2,3,\dots$ (4.2) $\displaystyle u_{n}(x;\mu)$
$\displaystyle=\frac{\pi^{-\frac{1}{4}}}{\sqrt{2^{n}n!}}e^{-\mu
x^{2}/2}H_{n}(\sqrt{\mu}x)$ (4.3)

where $H_{n}$ are the Hermite polynomials. First few Hermite polynomials are
$H_{0}(x)=1,H_{1}(x)=2x,H_{2}(x)=2x^{2}-2,H_{3}(x)=8x^{3}-12x,x_{4}(x)=16x^{4}-48x^{2}+12,H_{5}(x)=32x^{5}-160x^{3}+120x.$

The variational form of the above problem (4.1) is

$a(u(\mu),v;\mu)=\lambda(\mu)b(u(\mu),v),\quad\forall v\in H_{0}^{1}(\Omega)$
(4.4)

with

$a(u(\mu),v;\mu)=a_{1}(u,v)+\mu^{2}a_{2}(u,v)$

where

$a_{1}(u,v)=\frac{1}{2}\int_{\Omega}\frac{du}{dx}\frac{dv}{dx},\quad
a_{2}(u,v)=\frac{1}{2}\int_{\Omega}x^{2}uv,\quad b(u,v)=\int_{\Omega}uv.$

Numerical result

Figure 1: First eigenvalue plot with 95 percentage confidence interval with
mesh size $h=0.05$ for the problem (4.1) corresponding to 41 and 21 sample
points respectively.

Figure 2: First eigenvectors of the EVP (4.1), at $\mu=1.5,3.5,5.5$
respectively for the case of 21 sample points and mesh size $h=0.05$.

For the numerical test, we choose $\Omega=(-10,10)$ and the parameter space
$\mathcal{P}=[1,9]$. We have calculated the results of full order model (FOM)
with mesh size $h=0.05,0.02,0.01$ respectively. 41 uniformly distributed
sample parameters are selected from the parameter domain $[1,9]$. The first
eigenvectors are calculated at the selected sample points and put them into a
matrix column wise to form the snapshot matrix. We choose
$\mu=2.5,3.5,\dots,8.5$ as test points and calculated the first eigenvalues
and eigenvectors at these test points using data-driven (DD) reduced order
model. The first eigenvalues of FEM and the DD model at the test points are
reported in Table 1. Note that the sample points used for snapshot matrix and
the test points used for prediction are different. From the table one can see
that the results are matching with the FEM results.

Next, we choose 21 uniformly distributed points from the sample space $[1,9]$
and calculated the first eigenpairs at these sample points. Using the
eigenvalues at these points we trained a GPR model. Using the GPR model, we
predict the eigenvalues at the test points $\mu=2.5,3.5,4.5,5.5,6.5,7.5,8.5$.
The predicted eigenvalues are reported in Table 2 and compared them with the
first eigenvalues obtained by FEM. In this case also the eigenvalues are
matching with FEM eigenvalues. If we compare these results with the results of
Table 1 then one can see that for coarse mesh the DD results in Table 1 are
better than the results in Table 2. Thus if we take more points in the
parameter domain then the results of the DD model is better even if the mesh
is not fine enough.

$h$ | $\mu\to$ | 2.5 | 3.5 | 4.5 | 5.5 | 6.5 | 7.5 | 8.5
---|---|---|---|---|---|---|---|---
0.05 | DD | 1.2508 | 1.7528 | 2.2548 | 2.7568 | 3.2588 | 3.7608 | 4.2628
| FEM | 1.2511 | 1.7522 | 2.2537 | 2.7555 | 3.2577 | 3.7603 | 4.2632
0.02 | DD | 1.2501 | 1.7504 | 2.2508 | 2.7511 | 3.2514 | 3.7517 | 4.2521
| FEM | 1.2502 | 1.7504 | 2.2506 | 2.7509 | 3.2512 | 3.7516 | 4.2521
0.01 | DD | 1.2500 | 1.7501 | 2.2502 | 2.7503 | 3.2504 | 3.7504 | 4.2505
| FEM | 1.2500 | 1.7501 | 2.2501 | 2.7502 | 3.2503 | 3.7504 | 4.2505
Table 1:  First eigenvalues of (4.1) corresponding to $\mu=2.5,3.5,4.5,5.5,6.5,7.5,8.5$ respectively using Data-Driven (DD) and FEM with different mesh size of the domain when we consider 41 uniformly distributed sample points. $h$ | $\mu\to$ | 2.5 | 3.5 | 4.5 | 5.5 | 6.5 | 7.5 | 8.5
---|---|---|---|---|---|---|---|---
0.05 | DD | 1.2511 | 1.7529 | 2.2547 | 2.7565 | 3.2584 | 3.7602 | 4.2620
| FEM | 1.2511 | 1.7522 | 2.2537 | 2.7555 | 3.2577 | 3.7603 | 4.2632
0.02 | DD | 1.2502 | 1.7505 | 2.2508 | 2.7510 | 3.2513 | 3.7516 | 4.2519
| FEM | 1.2502 | 1.7504 | 2.2506 | 2.7509 | 3.2512 | 3.7516 | 4.2521
0.01 | DD | 1.2500 | 1.7501 | 2.2502 | 2.7503 | 3.2503 | 3.7504 | 4.2505
| FEM | 1.2500 | 1.7501 | 2.2501 | 2.7502 | 3.2503 | 3.7504 | 4.2505
Table 2:  First eigenvalues of (4.1) corresponding to
$\mu=2.5,3.5,4.5,5.5,6.5,7.5,8.5$ respectively using Data-Driven (DD) and FEM
with different mesh size of the domain when we consider 21 uniformly
distributed sample points.

Figure 3: First six RB coefficients with confidence interval for the case of
21 sample points and mesh size $h=0.05$.

In Figure 1, we have plotted the first eigenvalue for the DD-model with 95
percentage confidence interval for both the cases with 41, 21 points
respectively. we can see that the confidence interval is very narrow that
means, the predictions by the DD model is very accurate. In Figure 2, we have
shown the first eigenvector at test points $\mu=1.5,3.5,5.5$ respectively,
obtained by our DD model with 21 sample points and compared them with
corresponding FEM eigenvectors. From the Figure 2, it is evident that the DD
eigenvectors are coinciding with the FEM eigenvectors.

In Figure 3, we have shown the GPR plot of the first six projection
coefficients of the RB solutions. We have shown the test data (blue star) and
also plot the GPR mean function in magenta colour curve. For each GPR, we also
plot the 95 percentage confidence interval. Note that in the test set we take
two points namely $\mu=0.5,9.5$ which lies outside of the parameter interval
so at the end points the confidence interval is large. So our DD-model will
work best for all parameter within the parameter space.

####  4.1.2 Simple harmonic oscillator eigenvalue problem in 2D

Let us consider the eigenvalue problem defined on
$\Omega\subset\mathbb{R}^{2}$

$\left\\{\begin{array}[]{ll}-\frac{1}{2}\Delta
u(x,y)+\frac{1}{2}\mu^{2}(x^{2}+y^{2})u(x,y)&=\lambda
u(x,y)\quad\text{in}\quad\Omega=(-\frac{\pi}{2},\frac{\pi}{2})^{2}\\\
u&=0\quad\text{on}\quad\partial\Omega\end{array}\right.$ (4.5)

The analytic eigenvalues and eigenvectors of this problem are given by

$\displaystyle\lambda_{m,n}(\mu)$ $\displaystyle=(m+n+1)\mu\quad
n=0,1,2,3,\dots$ (4.6) $\displaystyle u_{m,n}(x,y;\mu)$
$\displaystyle=\frac{1}{\sqrt{2^{m+n}m!n!\pi}}e^{-\mu(x^{2}+y^{2})/2}H_{m}(\sqrt{\mu}x)H_{n}(\sqrt{\mu}y)$
(4.7)

where $H_{n}$ are the Hermite polynomials.

The variational form of the above problem (4.5) is

$a(u(\mu),v;\mu)=\lambda(\mu)b(u(\mu),v)\quad\forall v\in H_{0}^{1}(\Omega)$
(4.8)

with

$a(u,v;\mu)=a_{1}(u,v)+\mu^{2}a_{2}(u,v)$

where

$a_{1}(u,v)=\frac{1}{2}\int_{\Omega}\nabla u\cdot\nabla v,\quad
a_{2}(u,v)=\frac{1}{2}\int_{\Omega}(x^{2}+y^{2})uv,\quad
b(u,v)=\int_{\Omega}uv.$

Numerical Result:

h |  | $\delta\mu$
---
| Model
---
| $\mu=1.5$
---
| $2.5$
---
| $3.5$
---
| $4.5$
---
| $5.5$
---
| $6.5$
---
| $7.5$
---
0.1 | 0.4 | DD | 1.6683 | 2.5374 | 3.5040 | 4.5205 | 5.5275 | 6.5297 | 7.5457
|  | FEM | 1.6490 | 2.5317 | 3.5136 | 4.5167 | 5.5241 | 6.5334 | 7.5441
0.1 | 0.2 | DD | 1.6633 | 2.5318 | 3.5075 | 4.5214 | 5.5248 | 6.5302 | 7.5461
|  | FEM | 1.6490 | 2.5317 | 3.5136 | 4.5167 | 5.5241 | 6.5334 | 7.5441
0.05 | 0.4 | DD | 1.6668 | 2.5334 | 3.4970 | 4.5084 | 5.5088 | 6.5044 | 7.5125
|  | FEM | 1.6475 | 2.5279 | 3.5062 | 4.5043 | 5.5058 | 6.5081 | 7.5107
0.05 | 0.2 | DD | 1.6618 | 2.5279 | 3.5004 | 4.5090 | 5.5062 | 6.5052 | 7.5128
|  | FEM | 1.6475 | 2.5279 | 3.5062 | 4.5043 | 5.5058 | 6.5081 | 7.5107
Table 3:  First eigenvalues corresponding to different values of $\mu$ of 2D-
Oscillation Problem. Note that the exact first eigenvalue is equal to $\mu$.

For Numerical purpose we choose $\Omega=(-\frac{\pi}{2},\frac{\pi}{2})^{2}$
and the parameter space as $[1,9]$. We select two sets of uniformly
distributed parameters from the parameter space with parameter spacing
$\delta\mu=0.2,0.4$ respectively. For $\delta\mu=0.2$ we have 41 uniformly
distributed points and get 21 uniformly distributed points when
$\delta\mu=0.4$.

In Table 3, we have reported the first eigenvalue using the data-driven model
corresponding to the mesh size $h=0.1$ and $0.05$ at test points
$\mu=1.5,2.5,3.5,4.5,5.5,6.5,7.5$ and compared them with the results obtained
by P1-FEM. Note that the analytic first eigenvalue is equal to $\mu$. From the
Table 3 one can see that the results matching with the result of FEM. In
Figure 4, we have shown the eigenvalues used for the training and the
predicted eigenvalues for all the four cases corresponding to $h=0.1,0.05$ and
$\delta\mu=0.2,0.4$. Also, we have indicated the confidence interval for the
predicted eigenvalues.The confidence interval is very narrow.

In Figure 5, we have presented the training set and the predicted results for
the GPR corresponding to the six projected coefficients corresponding to first
eigenvectors for the 2D oscillator problem for the case
$h=0.05,\delta\mu=0.4$. The training sets are plotted as star and the
predicted values in magenta line and the 95 percentage of confidence interval
is shown as shaded region. In this case also the test points are
$0.5:0.1:9.5$. The end points of the test set lies outside of the parameter
interval, so the confidence interval is little wider at the end points.

The first eigenvectors obtained by our GPR-based DD model for the case
$h=0.05$ and $\delta\mu=0.4$, at the test points $\mu=1.5,3.5,5.5,7.5$ are
presented in Figure 6, and compared the with the eigenvectors obtained using
P1-FEM. Also, we have presented the error between the FEM eigenvectors and the
DD-based eigenvectors. We can see that the error is of order $10^{-4}.$

(a) $h=0.1,\delta\mu=0.4$

(b) $h=0.1,\delta\mu=0.2$

(c) $h=0.05,\delta\mu=0.4$

(d) $h=0.05,\delta\mu=0.2$

Figure 4: The GPR plot for the first eigenvalues with 95 percentage confidence
interval for the problem (4.5).

Figure 5: The GPR corresponding to the six RB coefficients for the problem
(4.5) with 21 sample points and $h=0.05$.

(a) FEM $(u_{1,h}$)

(b) DD ROM $(\tilde{u}_{1,h})$

(c) Error$(u_{1,h}-\tilde{u}_{1,h}$)

Figure 6: First eigenvectors of (4.5) at the test points
$\mu=1.5,3.5,5.5.5,7.5$ respectively with 21 sample points and $h=0.05$.

###  4.2 Nonlinear eigenvalue problem

Consider the nonlinear eigenvalue problem of the following form: given any
$\mu\in\mathscr{P}$ find $\lambda(\mu)\in\mathbb{R}$ and $u(\mu)\in V$
satisfying

$\left\\{\begin{aligned}
&a(u(\mu),v)+\mu^{2}\int_{\Omega}g(u(\mu))v=\lambda(\mu)m(u(\mu),v),\quad\forall
v\in V\\\ &(u(\mu),u(\mu))=1,\end{aligned}\right.$ (4.9)

where $g$ is a nonlinear function $u$. To apply Newton’s method for this
nonlinear problem let us write

$u^{k+1}=u^{k}+\delta u,\quad\lambda^{k+1}=\lambda^{k}+\delta\lambda$

then (4.9) reduces to

$\left\\{\begin{aligned} &a(u^{k}+\delta
u,v)+\mu^{2}\int_{\Omega}g(u^{k}+\delta
u)v=(\lambda^{k}+\delta\lambda)(u^{k}+\delta u,v),\quad\forall v\in V\\\
&(u^{k}+\delta u,u^{k}+\delta u)=1,\end{aligned}\right.$ (4.10)

Using Taylor’s theorem to the function $g$ and neglecting the higher order
terms we get

$\left\\{\begin{aligned} &a(u^{k},v)+a(\delta
u,v)+\mu^{2}\int_{\Omega}g(u^{k})v+\mu^{2}\int_{\Omega}\delta
ug^{\prime}(u^{k})v=\lambda^{k}(u^{k},v)+\delta\lambda(u^{k},v)+\lambda^{k}(\delta
u,v),\quad\forall v\in V\\\ &(u^{k}+,u^{k})+2(\delta
u,u^{k})=1,\end{aligned}\right.$ (4.11)

On simplifying we get,

$\left\\{\begin{aligned} &a(\delta u,v)+\mu^{2}\int_{\Omega}\delta
ug^{\prime}(u^{k})v-\lambda^{k}(\delta
u,v)-\delta\lambda(u^{k},v)=\lambda^{k}(u^{k},v)-a(u^{k},v)-\mu^{2}\int_{\Omega}g(u^{k})v,\quad\forall
v\in V\\\ &2(\delta u,u^{k})=1-(u^{k},u^{k}),\end{aligned}\right.$ (4.12)

Let $V_{h}\subset V$ is the finite dimensional space. Then the Galerkin FE
method corresponding to the non-liner problem reads: given
$(u_{h}^{k},\lambda^{k})$ find $(u_{h}^{k+1},\lambda^{k+1})=(u_{h}^{k}+\delta
u_{h},\lambda^{k}+\delta\lambda_{h})$ such that

$\left\\{\begin{aligned} &a(\delta u_{h},v_{h})+\mu^{2}\int_{\Omega}\delta
u_{h}g^{\prime}(u_{h}^{k})v_{h}-\lambda^{k}(\delta
u_{h},v_{h})-\delta\lambda_{h}(u_{h}^{k},v_{h})=\lambda^{k}(u_{h}^{k},v_{h})-a(u_{h}^{k},v_{h})-\mu^{2}\int_{\Omega}g(u_{h}^{k})v_{h},\quad\forall
v_{h}\in V_{h}\\\ &2(\delta
u_{h},u_{h}^{k})=1-(u_{h}^{k},u_{h}^{k}),\end{aligned}\right.$ (4.13)

Let $\\{\varphi_{j}:j=1,\dots,N_{h}\\}$ be the basis of $V_{h}$. Now
expressing $u_{h}^{k}=\sum\limits_{j=1}^{N_{h}}u_{j,h}^{k}\varphi_{j}$ and
$\delta u_{h}=\sum\limits_{j=1}^{N_{h}}d_{j}\varphi_{j}$ we get the matrix
form of the above discretized equation as:

$\displaystyle\begin{bmatrix}A&-\boldsymbol{c}\\\ 2\boldsymbol{c}^{T}&0\\\
\end{bmatrix}\begin{bmatrix}\boldsymbol{d}\\\
\delta\lambda_{h}\end{bmatrix}=\begin{bmatrix}\boldsymbol{b}\\\
1-s\end{bmatrix}$

where the matrix $A=A_{1}+\mu^{2}M_{g^{\prime}}-\lambda^{k}M$ with

$A_{1}(i,j)=a(\varphi_{j},\varphi_{i}),\quad
M(i,j)=(\varphi_{j},\varphi_{i}),\quad
M_{g^{\prime}}(i,j)=\int_{\Omega}g^{\prime}(u^{k})\varphi_{j}\varphi_{i},$

and $\boldsymbol{b},\boldsymbol{c}$ are the vectors satisfying
$\boldsymbol{c}=M\boldsymbol{u}^{k}_{h}(\mu)$,
$\boldsymbol{b}=\lambda^{k}M\boldsymbol{u}_{h}^{k}-A_{1}\boldsymbol{u}_{h}^{k}-\mu^{2}\boldsymbol{g}$
with $\boldsymbol{g}=(g_{1},\dots,g_{N_{h}})^{T}$ and
$g_{j}=\int_{\Omega}g(u^{k}_{h})\varphi_{j}$ and the scalar
$s=\boldsymbol{u}_{h}^{k}(\mu)^{T}M\boldsymbol{u}_{h}^{k}(\mu)$ with
$\boldsymbol{u}_{h}^{k}=(u_{1,h}^{k},\dots,u_{N_{h},h}^{k})^{T}$. Solving the
above equation we will get $\delta u_{h}$ and $\delta\lambda$. Thus the
updated eigenvalue is $\lambda^{k+1}=\lambda^{k}+\delta\lambda$ and updated
eigenvectors $u^{k+1}(\mu)=u^{k}(\mu)+\delta u_{h}(\mu)$. For non-linear
eigenvalue problem we only find the first eigenvalue and first eigenvector.

Numerical Result:

For numerical purpose, we consider the problem on
$\Omega=(0,1)\subset\mathbb{R}$ and the parameter space is
$\mathscr{P}=[1,9]$. For given $\mu\in\mathscr{P}$, we solve the problem (4.9)
for the following particular choice:

$a(w,v)=\int_{\Omega}\frac{dw}{dx}\frac{dv}{dx},\quad
m(w,v)=\int_{\Omega}wvdx,\quad g(w)=|w|^{7/3}w.$

We select two sets of sample parameters for calculating snapshots. First we
take 41 uniformly distributed points and then 21 uniform points in the
parameter space $\mathscr{P}$, that is the parameter step size will be
$\delta\mu=0.2,0.4$ respectively. First eigenvalues at the test parameters
$\mu=1.5:1:6.5$ obtained by our GPR-based data-driven model is reported in
Table 4 and compared them with the first eigenvalues obtained using P1-FEM. We
have reported the results with mesh size $h=0.05$ and $h=0.01$. The
eigenvalues obtained by FEM and the data-driven model are almost same. Note
that because of non-linearity, the FEM methods is iterative in nature and it
takes few iterations to get the solution to converge but in the data-driven
model we are getting directly the solution by evaluating the regression
functions.

h |  | $\delta\mu$
---
| Method
---
| 1.5
---
| 2.5
---
| 3.5
---
| 4.5
---
| 5.5
---
| 6.5
---
0.05 | 0.4 | FEM | 13.4620 | 19.5299 | 28.1425 | 39.0260 | 52.0268 | 67.0749
|  | DD | 13.4926 | 19.6184 | 28.1278 | 38.9666 | 51.9936 | 67.0948
0.05 | 0.2 | FEM | 13.4620 | 19.5299 | 28.1425 | 39.0260 | 52.0268 | 67.0749
|  | DD | 13.4807 | 19.5773 | 28.1241 | 38.9963 | 52.0197 | 67.0857
0.01 | 0.4 | FEM | 13.4553 | 19.5367 | 28.1596 | 39.0516 | 52.0624 | 67.1239
|  | DD | 13.4863 | 19.6256 | 28.1447 | 38.9918 | 52.0289 | 67.1440
0.01 | 0.2 | FEM | 13.4553 | 19.5367 | 28.1596 | 39.0516 | 52.0624 | 67.1239
|  | DD | 13.4743 | 19.5842 | 28.1409 | 39.0218 | 52.0552 | 67.1348
Table 4: 1st eigenvalues of the nonlinear problem (4.9) obtained using DD
model with mesh size $h=0.05,0.01$.

In Figure 7, we have shown the first eigenvalue plot using FEM and our data-
driven model when the mesh parameter $h=0.05$ and $\delta\mu=0.2,0.4$
respectively. Also, we have indicated the confidence interval for the GPR
corresponding to the first eigenvalue. The confidence interval is very narrow,
that is the eigenvalues obtained by FEM and data-driven model are very close
to each other. That is the eigenvalues we get using the GPR-based DD model is
very accurate. Also, in Figure 8, we have plotted the mean function of the GPR
corresponding to the four projected coefficients of the RB solutions with 95
percent confidence of interval for the choice $h=0.05$ and the sample points
are $\mu_{tr}=1:0.4:9$ that is the snapshot matrix contains all the first
eigenvectors at these sample points.

(a) $\mu_{tr}=1:0.2:9$

(b) $\mu_{tr}=1:0.4:9$

Figure 7: The GPR for the first eigenvalues of the nonlinear problem (4.9)
with different number of snapshots with $h=0.05$.

Figure 8: The four GPR corresponding to the four co-efficents of the first
reduced eigenvector and their confidence interval of the EVP (4.9), when the
snapshots are calculated at $\mu=1:0.4:9$ and mesh size $h=0.05$.

(a) $\mu=1.5$

(b) $\mu=2.5$

(c) $\mu=3.5$

(d) $\mu=4.5$

(e) $\mu=5.5$

(f) $\mu=6.5$

Figure 9: First eigenvectors of the nonlinear problem (4.1) using DD model and
FEM with $h=0.05$ and $\delta\mu=0.4$.

In Figure 9, we have shown the first eigenvector obtained by our DD-model at
the test points $\mu=1.5:1:6.5$ and also plotted the corresponding eigenvector
using FEM for the case $h=0.05$ and $\delta\mu=0.4$. From the Figure 9, one
can see that the eigenvectors using FEM and DD-model are coinciding to each
other at all the test points.

###  4.3 Non-affine parameter dependent eigenvalue problem

In this section, we consider the eigenvalue problem (EVP) that are not affine
parameter dependent. For these problem, we can not apply the standard
projection based reduced method. But as our model is based on the data, we can
easily get the reduced solution using the trained GPR.

####  4.3.1 Non-affine eigenvalue problem with single parameter

Let us consider the following eigenvalue problem

$\left\\{\begin{array}[]{ll}-\Delta
u(\mu)=\lambda(\mu)e^{-\mu(x^{2}+y^{2})}u(\mu)&\textrm{ in
}\Omega=(0,1)^{2}\\\ u(\mu)=0&\textrm{ on }\partial\Omega\end{array}\right.$
(4.14)

where the parameter $\mu\in[1,8]$. Note that the RHS of the equation is non-
affine. For the training our GPR model, we have calculated the eigensolutions
of the EVP at the sample points $\mu=1:0.4:8$ and $\mu=1:0.2:8$ respectively.
Also, we have experimented with mesh size $h=0.05,0.01$. Using the technique
describe in section 3.3, we constructed one GPR for each projection
coefficients of the RB solution and one GPR for the eigenvalue. We have
reported the eigenvalues obtained using our model at the test points
$\mu=1.5,3.5,5.5,7.5$ in Table 5 for all the four cases corresponding to
$h=0.05,0.01$ and $\delta\mu=0.4,0.2$. To compared the eigenvalues we also
reported the P1-FEM based eigenvalues. Note that the test points are different
from the sample points. Eigenvalues are very close to the FEM eigenvalues.

In Figure 10, we have shown the plot for the first eigenvalue using FEM and
our data-driven model when the mesh parameter $h=0.05$ and $\delta\mu=0.2,0.4$
respectively. Also, we have shown the confidence interval for the GPR of the
eigenvalues. The confidence interval is very narrow, that is the eigenvalues
obtained by FEM and data-driven model will not vary much. That is the
eigenvalues we get using the GPR-based DD model is very accurate.

Also, in Figure 11, we have plotted all the five coefficients of the RB
solutions with 95 percent confidence of interval for the choice $h=0.05$ and
the sample points are $\mu_{tr}=1:0.4:8$ that is the snapshot matrix contains
all the first eigenvectors at these sample points. The test points are
$\mu=0.5:1:8.5$, the end points are outside of the parameter space. So, the
confidence interval at the end points are little larger.

In Figure 12, we displayed the 1st eigenvectors obtained by using FEM and our
DD-model at the test points. The eigenvectors are looks same. To check how
close it is, we have shown the error between them and the error is of order
$10^{-4}$.

h |  | $\delta\mu$
---
| Method
---
| $\mu=1.5$
---
| $\mu=3.5$
---
| $\mu=5.5$
---
| $\mu=7.5$
---
0.05 | 0.4 | FEM | 39.73473805 | 72.62248873 | 108.48892187 | 145.54749238
|  | DD | 39.81667386 | 72.57625001 | 108.50118905 | 145.49823357
0.05 | 0.2 | FEM | 39.73473805 | 72.62248873 | 108.48892187 | 145.54749238
|  | DD | 39.79566310 | 72.57602363 | 108.50603170 | 145.54219385
0.01 | 0.4 | FEM | 39.62169633 | 72.32111463 | 107.85749755 | 144.42856653
|  | DD | 39.70506978 | 72.27383282 | 107.87538316 | 144.37709615
0.01 | 0.2 | FEM | 39.62169633 | 72.32111463 | 107.85749755 | 144.42856653
|  | DD | 39.68283223 | 72.27478421 | 107.87527246 | 144.42228776
Table 5: Comparison of 1st eigenvalues of the EVP (4.14) using FEM and DD
Model.

(a) $\mu_{tr}=1:0.4:8$

(b) $\mu_{tr}=1:0.2:8$

Figure 10: The GPR for the first eigenvalues of the nonlinear problem (4.14)
with different number of snapshots with $h=0.05$.

(a) 1st coefficient

(b) 2nd coefficient

(c) 3rd coefficient

(d) 4th coefficient

(e) 5th coefficient

Figure 11: The GPR plot for the coefficients for the problem (4.14) with
$h=0.05$ and $\delta\mu=0.4$ and test points are $\mu=0.5:1:8.5$.

Figure 12: 1st eigenvectors of problem (4.14) using DD model and FEM with
$h=0.05$ and $\delta\mu=0.4$.

####  4.3.2 Non-affine eigenvalue problem with multiple parameter

Consider the eigenvalue problem: given
$\boldsymbol{\mu}=(\mu_{1},\mu_{2})\in\mathscr{P}$, find
$(\lambda(\boldsymbol{\mu}),u(\boldsymbol{\mu}))$ such that

$\left\\{\begin{array}[]{ll}-\Delta
u&=\lambda(\boldsymbol{\mu})\varepsilon(x;\boldsymbol{\mu})u\quad\text{in}\quad\Omega\\\
u(\mu)=0&\textrm{ on }\partial\Omega\end{array}\right.$ (4.15)

where
$\varepsilon(x;\boldsymbol{\mu})=\epsilon_{1}\chi_{\Omega_{1}(\boldsymbol{\mu})}(x)+\epsilon_{2}\chi_{\Omega_{2}(\boldsymbol{\mu})}(x)$
with $\Omega_{1}(\boldsymbol{\mu})\cup\Omega_{2}(\boldsymbol{\mu})=\emptyset$
and
$\Omega_{1}(\boldsymbol{\mu})\cap\Omega_{2}(\boldsymbol{\mu})=\Omega=(-1,1)^{2}$.
For numerical result $\epsilon_{1}=0.1,\epsilon_{2}=0.2$. The interface
separating the $\Omega_{1}(\boldsymbol{\mu})$ and
$\Omega_{2}(\boldsymbol{\mu})$ is

$x_{1}=\mu_{1}sin(\mu_{2}\pi x_{2})$

with the parameter $\boldsymbol{\mu}=(\mu_{1},\mu_{2})$ ranging in
$\mathscr{P}=[0.1,0.2]\times[1,8].$ Introducing a linear transition between
$\epsilon_{1}$ and $\epsilon_{2}$ near the interface like it is done in
[Fumagallietal16], we get the function presented in Figure 13 for the
parameter $\mu_{1}=0.2$ and $\mu_{2}=1$ and $8$, respectively.

$(\mu_{1},\mu_{2})$ |  | Method
---
| Eigenvalue-I
---
| Eigenvalue-II
---
| Eigenvalue-III
---
| Eigenvalue-IV
---
(0.15,1.5) | FEM | 31.48699165 | 75.22816969 | 89.16775007 | 143.71355581
| DD | 31.49811030 | 75.28078739 | 89.10207780 | 142.38589826
(0.15,3.5) | FEM | 31.54972648 | 75.32810691 | 89.03858404 | 144.08480023
| DD | 31.56475170 | 75.35994556 | 88.96400162 | 143.35827613
(0.15,5.5) | FEM | 31.56473098 | 75.38499108 | 89.05984549 | 144.34148451
| DD | 31.57095675 | 75.40394116 | 89.00349397 | 143.69597428
(0.15,7.5) | FEM | 31.58159550 | 75.42566021 | 89.01447975 | 144.45513191
| DD | 31.58010472 | 75.42981927 | 89.00210188 | 143.80005269
Table 6:  First four eigenvalues using FEM and DD-model at the test points of
the problem (4.15) with mesh size $h=0.1$.

(a) $\mu_{1}=0.2,\mu_{2}=1$

(b) $\mu_{1}=0.2,\mu_{2}=8$

Figure 13: The weight function $\varepsilon(\boldsymbol{\mu})$ in (4.15) with
two different $\boldsymbol{\mu}$.

(a) FEM

(b) DD ROM

(c) Error

Figure 14: First eigenvectors of EVP (4.15) at parameters $\mu=(0.15,3.5)$ and
$\mu=(0.15,7.5)$ respectively, using FEM and our data-driven model.

In this example, we can not write the RHS and hence the corresponding linear
form as affine parametric dependence form because of the nature of the
function $\varepsilon(\boldsymbol{\mu})$. We solve this eigenvalue problem
using GPR based data-driven model.

Numerical result

For this example we have taken 2 uniformly distributed points in the parameter
space $\mu_{1}$ and 36 uniformly distributed points $(\mu_{2}=1:0.2:8)$ in the
parameter space of $\mu_{2}$ and take their cross product and get 72 sample
points in the parametric space. The snapshot matrix is constructed by taking
the first eigenvector at those sample points as columns. As a test set we
choose $\mu_{1}=0.15$ and $\mu_{2}=1,2.5,4,4.5,6.5,8$ respectively which is
different from from those taken to generate snapshots.

In Table 6, we have reported first four eigenvalues of the eigenvalue problem
obtained by using FEM and our data-driven model for the test parameter
mentioned in the last paragraph. The predicted eigenvalues by the data driven
model are matching with the eigenvalues obtained by using FEM. Note that for
construction of GPR we do not use any information of the snapshot matrix, so
we can predict more than one eigenvalues using the data.

In Figure 14, we have shown the first eigenvectors obtained by FEM and data-
driven model for the parameter $\boldsymbol{\mu}=(0.15,3.5)$ and
$\boldsymbol{\mu}=(0.15,7.5)$ and also shown the relative error between them
and it is of order $10^{-4}$.

###  4.4 Results of EVP with intersecting Eigenvalues

Let us consider the following eigenvalue problem

$\left\\{\begin{array}[]{ll}-\nabla\cdot(A\nabla u)&=\lambda
u\quad\text{in}\quad\Omega=[-1,1]^{2}\\\ u(\mu)=0&\textrm{ on
}\partial\Omega\end{array}\right.$ (4.16)

with the matrix $A=\begin{bmatrix}1&0\\\ 0&1+\mu\end{bmatrix}$ where we choose
the parameter space as $\mathcal{P}=[-0.9,0.9]$.The analytic eigenvalues are

$\displaystyle\lambda_{m,n}(\mu)=\frac{\pi^{2}}{4}(m^{2}+(1+\mu)n^{2}),\quad
m,n=1,2,\dots$

and the eigenvectors are

$\displaystyle u_{m,n}(\mu)=cos(\frac{m\pi}{2}x)cos(\frac{n\pi}{2}y).$

In Figure 15, we plotted the first six eigenvalues obtained using FEM. From
the figure, we can see that the eigenvalues other than the first eigenvalues
are intersecting. We will investigate, how our model works for the
intersecting eigenvalues and eigenvectors. We present the numerical results
for the first and the third eigenvalues and eigenvectors for this problem by
data-driven model.To calculate the third eigenvector, we consider the snapshot
matrix consists of the third eigenvectors at the sample parameters
$\mu_{1},\mu_{2},\dots,\mu_{n_{s}}$ as the columns.

Figure 15: First six sorted eigenvalues of the problem (4.16). h |  | Method
---
| $\mu=-0.75$
---
| $\mu=-0.25$
---
| $\mu=0.25$
---
| $\mu=0.75$
---
0.1 | FEM | 3.09172930 | 4.32853369 | 5.56526834 | 6.80197424
| DD | 3.09175362 | 4.32850239 | 5.56525117 | 6.80199994
0.05 | FEM | 3.08606437 | 4.32052203 | 5.55496589 | 6.78940395
| DD | 3.08606913 | 4.32051579 | 5.55496244 | 6.78940910
0.01 | FEM | 3.08432204 | 4.31805168 | 5.55178071 | 6.78550948
| DD | 3.08432225 | 4.31805140 | 5.55178056 | 6.78550971
Table 7: First eigenvalues of (4.16) for different mesh using DD model with sample points $\mu_{tr}=-0.9:0.1:0.9$ and FEM model. h |  | Method
---
| $\mu=-0.75$
---
| $\mu=-0.25$
---
| $\mu=0.25$
---
| $\mu=0.75$
---
0.1 | -FEM | 8.14338931 | 11.78888922 | 14.89196477 | 19.85303433
| DD | 8.21778394 | 11.81042527 | 14.91670069 | 19.85137938
0.05 | FEM | 8.05008647 | 11.73700317 | 14.82575077 | 19.76686171
| DD | 8.12327908 | 11.75458889 | 14.85046418 | 19.76559935
0.01 | FEM | 8.02024667 | 11.72081569 | 14.80523805 | 19.74028492
| DD | 8.09313087 | 11.73711527 | 14.82993019 | 19.73915595
Table 8: 3rd eigenvalues of (4.16) for different mesh using DD model with
sample points $\mu_{tr}=-0.9:0.1:0.9$ and FEM model.

(a) FEM

(b) DD

(c) Error

Figure 16: 1st eigenvectors corresponding to mesh size $h=0.05$ using DD model
with 1 POD basis where the training set is $\mu_{tr}=-0.9:0.1:0.9$ and FEM.

For the experiment, We took $19$ uniformly distributed sample points on the
parameter space $[-0.9,0.9]$, that is the parameter step size is $0.1$. Then
using the eigenvalues and eigenvectors at these points we construct the GPR
following the idea describe in section 3.3. We have presented the eigenvalues
and eigenvectors at four test points $\mu=-0.75,-0.25,0.25,0.75$.

In Table 7, we have reported the first eigenvalues obtained by our DD-model at
these four points for the mesh size $h=0.1,0.05$ and $0.01$ respectively.
Also, for the comparison we have reported the eigenvalues using the FEM. From
the table one can see that the eigenvalues are the DD-based eigenvalues are
matching with the FEM eigenvalues at least up to four decimal places.In Figure
16 we have displayed the first eigenvector at those four parameters obtained
by our model and FEM model. Also, we have displayed the point-wise error
between the eigenvectors obtained by FEM and our DD-model. The error is of
order $10^{-4}-10^{-5}$. Note that, the first eigenvector is independent of
the parameter $\mu$, that is the reason why the criterion (3.13), selects only
one POD basis for this example.

Now to investigate how our model performs when the crossing occurs, we
calculate the third eigenvalues and eigenvectors. The process of constructing
the GPR model is same but instead of first eigenvalues and first RB
eigenvectors we used the third eigenvalues and third reduced eigenvectors as
the training data at the sample points. The training points and the GPR for
the three co-efficients are presented in Figure17 for the mesh size $h=0.05$.
To draw the GPR we have use the test points $\mu=-0.9:0.05:0.9$. Also, we have
shown the confidence interval for all the three reduced coefficients. For,
this case the confidence interval is little large, so the error is higher than
the first eigenvalues and eigenvectors.

(a) First coefficient

(b) Second coefficient

(c) Third coefficient

Figure 17: Three GPR corresponding to the three RB-coefficients for the 3rd
eigenvector of the problem (4.16) with $h=0.05$. The training set is
$\mu_{tr}=-0.9:0.1:0.9$.

(a) $h=0.1$

(b) $h=0.05$

(c) $h=0.01$

Figure 18: The GPR corresponding to 3rd eigenvalues of (4.16) with mesh size
$h=0.1,0.05$ and $0.01$ respectively. The training set is
$\mu_{tr}=-0.9:0.1:0.9$.

In Table 8, we have reported the third eigenvalues obtained by our DD-model at
these four points for the mesh size $h=0.1,0.05$ and $0.01$ respectively.
Also, for the comparison we have reported the eigenvalues using the FEM. From
the table one can see that the eigenvalues of our DD-model are very closed to
the eigenvalues of the FEM method. In Figure 19 we have displayed the third
eigenvector at those four parameters obtained by our model and FEM model.
Also, we have displayed the point-wise error between the eigenvectors obtained
by FEM and our DD-model. The error is of order $10^{-2}-10^{-3}$. Note that,
although the third eigenvector is independent of the parameter $\mu$ but
because of the crossing of the eigenvalues, three different eigenvectors
occurs in the snapshot matrix, that is the reason why the criterion (3.13)
gives three POD basis for calculating third eigenvector.

In Figure 18, we have plotted the 3rd eigenvalues used for the training as
star and also plotted the mean of GPR using the eigenvalues at the test points
$-0.9:0.05:0.9$ in magenta line corresponding to mesh size $h=0.1,0.05$ and
$0.01$. Also, we have shown the 95 percentage confidence interval for the 3rd
eigenvalues.The confidence interval is very narrow, that means the eigenvalues
we obtained using the GPR model is very closed to that of the original one.

(a) FEM

(b) DD ROM

(c) Error

Figure 19: 3rd eigenvectors corresponding to mesh size $h=0.05$ using DD model
with 4 POD basis where the training set is $\mu_{tr}=-0.9:0.1:0.9$ and FEM.

Finally, we have tested our data-driven model for calculating multiple
eigenfunctions simultaneously by using the procedure described in subsection
3.4. We have calculated first three eigenvalues and eigenvectors
simultaneously for the problem (4.16). In Figure 20, we shave shown the plot
for the mean of the GPR and the training set, corresponding to the first three
eigenvalues, with 95 percentage of confidence interval. We observed that the
training data are lying in the curve of the mean GPR for all the three cases.
The confidence interval for the first eigenvalue is very narrow and for the
other two confidence interval is little larger than the fist one. All the
three eigenvalues obtained by DD-model at the test points
$\mu=-0.75,-0.25,0.25,0.75$ corresponding to mesh size $h=0,05$ are reported
in Table 9 and compared them with the FEM result. From the table one can see
that the first eigenvalues obtained by data-driven model are more accurate
than the second and third eigenvalues.

| $\mu$
---
| Method
---
| 1st eigenvalue
---
| 2nd eigenvalue
---
| 3rd eigenvalue
---
-0.75 | FEM | 3.08606437 | 4.94333422 | 8.05008647
| DD | 3.08606913 | 4.93514768 | 8.12327908
-0.25 | FEM | 4.32052203 | 9.88459549 | 11.73700317
| DD | 4.32051579 | 9.91536968 | 11.75458889
0.25 | FEM | 5.55496589 | 12.97341757 | 14.82575077
| DD | 5.55496244 | 13.00419699 | 14.85046418
0.75 | FEM | 6.78940395 | 14.20979040 | 19.76686171
| DD | 6.78940910 | 14.20160725 | 19.76559935
Table 9: First three eigenvalues of (4.16) for different mesh $h=0.05$ using
DD model with sample points $\mu_{tr}=-0.9:0.1:0.9$ and FEM model.

In Figure 21, we have presented the plots of the mean GPR corresponding to the
four reduced coefficients ans the training data when for the simultaneous
case. Also, we have plotted the 95 percentage of confidence interval for each
GPR. The GPR for the first and third coefficients are little large. Note that
for data-driven model the sign for the snapshots are very important. We have
chosen the sign of the snapshot following the following algorithm [Pau07a],
where we used notation of MATLAB:

$S_{s}=[u_{1,h}(\mu_{1});u_{2,h}(\mu_{1});u_{3,h}(\mu_{1})]$

for $j=2:n_{s}$ do

$s_{j}=[]$

for $i=1:3$ do

$e_{1}=\|u_{i,h}(\mu_{j})-u_{i,h}(\mu_{j-1})\|$,

$e_{2}=\|u_{i,h}(\mu_{j})+u_{i,h}(\mu_{j-1})\|$

if $e_{1}\geq e_{2}$ then

$s_{j}=[s_{j};-u_{i,h}(\mu_{j})]$

else

$s_{j}=[s_{j};u_{i,h}(\mu_{j})]$

end if

end for

$S_{s}=[S_{s},s_{j}]$

end for

In Figure 22-Figure 25, we have shown the first three eigenvectors obtained by
simultaneous data-driven model and compared them with that of FEM at the test
points $\mu=-0.75,-0.25,0.25$ and $0.75$ respectively. From Figures one can
see that the DD-based eigenvectors are matching with FEM-based eigenvectors
for all the four points. Note that in this case we only had to train seven
GPR, three corresponding to the three eigenvalues and four corresponding to
the reduced coefficients.

(a) GPR for 1st EV

(b) GPR for 2nd EV

(c) GPR for 3rd EV

Figure 20: The GPR corresponding to first three eigenvalues of (4.16) with
mesh size $h=0.05$. The training set is $\mu_{tr}=-0.9:0.1:0.9$.

(a) GPR for 1st coefficient

(b) GPR for 2nd coefficient

(c) GPR for 3rd coefficient

(d) GPR for 4th coefficient

Figure 21: The GPR corresponding to first three eigenvalues of (4.16) with
mesh size $h=0.05$ for simultaneous case. The training set is
$\mu_{tr}=-0.9:0.1:0.9$.

(a) 1st eigenvector (FEM)

(b) 2nd eigenvector (FEM)

(c) 3rd eigenvector (FEM)

(d) 1st eigenvector (DD)

(e) 2nd eigenvector (DD)

(f) 3rd eigenvector (DD)

Figure 22: Comparison of first three eigenvectors (4.16) at $\mu=-0.75$ with
$h=0.05$.

(a) 1st eigenvector (FEM)

(b) 2nd eigenvector (FEM)

(c) 3rd eigenvector (FEM)

(d) 1st eigenvector (DD)

(e) 2nd eigenvector (DD)

(f) 3rd eigenvector (DD)

Figure 23: Comparison of first three eigenvectors (4.16) at $\mu=-0.25$ with
$h=0.05$.

(a) 1st eigenvector (FEM)

(b) 2nd eigenvector (FEM)

(c) 3rd eigenvector (FEM)

(d) 1st eigenvector (DD)

(e) 2nd eigenvector (DD)

(f) 3rd eigenvector (DD)

Figure 24: Comparison of first three eigenvectors (4.16) at $\mu=0.25$ with
$h=0.05$.

(a) 1st eigenvector (FEM)

(b) 2nd eigenvector (FEM)

(c) 3rd eigenvector (FEM)

(d) 1st eigenvector (DD)

(e) 2nd eigenvector (DD)

(f) 3rd eigenvector (DD)

Figure 25: Comparison of first three eigenvectors (4.16) at $\mu=0.75$ with
$h=0.05$.