# Bayesian Spectral Deconvolution of X-Ray Absorption Near Edge Structure
Discriminating High- and Low-Energy Domains

Shuhei Kashiwamura1, Shun Katakami2, Ryo Yamagami3, Kazunori Iwamitsu4,
Hiroyuki Kumazoe5, Kenji Nagata6, Toshihiro Okajima7, Ichiro Akai5, Masato
Okada1,2<EMAIL_ADDRESS>1Graduate School of Science, The University
of Tokyo, Bunkyo, Tokyo 113-0033, Japan 2Graduate School of Frontier Sciences,
The University of Tokyo, Kashiwa, Chiba 277-8561, Japan
3Graduate School of Science and Technology, Kumamoto University 860-8555,
Japan
4Technical Division, Kumamoto University, Kumamoto 860-8555, Japan
5Institute of Industrial Nanomaterials, Kumamoto University, Kumamoto
860-8555, Japan
6Research and Services Division of Materials Data and Integrated System,
National Institute for Materials Science, Tsukuba, Ibaraki 305-0047, Japan
7Aichi Synchrotron Radiation Center, Seto, Aichi 489-0965, Japan

###### Abstract

In this paper, we propose a Bayesian spectral deconvolution considering the
properties of peaks in different energy domains. Bayesian spectral
deconvolution regresses spectral data into the sum of multiple basis
functions. Conventional methods use a model that treats all peaks equally.
However, in X-ray absorption near edge structure (XANES) spectra, the
properties of the peaks differ depending on the energy domain, and the
specific energy domain of XANES is essential in condensed matter physics. We
propose a model that discriminates between the low- and high-energy domains.
We also propose a prior distribution that reflects the physical properties. We
compare the conventional and proposed models in terms of computational
efficiency, estimation accuracy, and model evidence. We demonstrate that our
method effectively estimates the number of transition components in the
important energy domain, on which the material scientists focus for mapping
the electronic transition analysis by first-principles simulation.

††preprint: APS/123-QED

##  I Introduction

X-ray absorption near edge structure (XANES) appears in the energy region of
$\sim$50 eV around the X-ray absorption edge. It comes from the X-ray
absorption due to electronic transitions from a core level to unoccupied band
states for the element selected by the edge energy of X-ray absorption.[1]
XANES spectra are measured in condensed matter research because they offer
crucial physicochemical information. They markedly change with the atom
valence, the crystal field dependent on microscopic symmetry, and chemical
bonding states. [2] The experimental results of XANES are compared with the
density of states (DoS)[3, 4] obtained by first-principles calculations, and
the electronic states are discussed.

XANES spectra contain multiple peaks and step structures. It is necessary to
decompose the spectrum into these components to extract information about
electronic transitions. However, estimating the number of peaks and each
peak’s width, position, and intensity is generally challenging since the
spectrum is often complex and multimodal. A naive approach is least-squares
fitting by the steepest descent method, but it cannot avoid the local
minima.[5] Local minima show results that depend on the initial values, making
the analysis arbitrary.

An effective approach for the analysis of XANES spectra is Bayesian spectral
deconvolution[6], where the following spectral parameters can be estimated:
position, width, and intensity of the basis function. [7] In previous
studies,[7, 6] the exchange Monte Carlo method (EMC)[8] was used to get out of
the local solution and approach the optimal global solution. EMC makes it
possible to sample from the posterior distribution of spectral parameters. The
number of peaks and noise variance can also be estimated by maximizing a
function called the marginal likelihood.[9]

Previous studies on Bayesian spectral deconvolution used the regression
function treating all peaks equally.[7, 9, 6] However, the properties and
physical background of XANES spectra differ in the high- and low-energy
domains. Conventional methods have not fully utilized these differences.

In this work, we propose the energy-domain-aware regression model for Bayesian
spectral deconvolution. It incorporates the prior knowledge that XANES spectra
have sharp peaks in the low-energy domain and broad peaks in the high-energy
domain. Our proposed model enables us to discuss the number of peaks in each
energy domain separately and shows its superiority in convergence, sampling
efficiency, estimation accuracy, and model evidence compared with the
conventional model. This was confirmed via synthetic data analysis.

This paper is organized as follows. In Sect. 2, we introduce the framework of
the Bayesian spectral deconvolution. We describe our model and show that our
model contains the conventional model depending on the design of prior
distributions. In Sect. 3, we compare the conventional and proposed models and
verify our framework via synthetic data analysis. In Sect. 4, we present our
discussion and conclude this paper.

##  II Framework

Figure 1:  Synthetic model and data used in our analysis. (a) Gray dots show
synthetic data. (b) Graph of spectrum model (black line) consisting of
individual Gaussian peaks (colored lines) and the absorption edge with white
line (gray line).

###  II.1 Model

In XANES spectra, peaks in the low-energy domain reflect discrete electronic
transitions from the atomic core levels to unoccupied states or valence
bands.[1] In contrast, peaks in the high-energy domain reflect transitions to
continuum states, also observed in extended X-ray absorption fine structure
(EXAFS). [10] Thus, peaks in the low-energy domain are focused for mapping to
the electronic transition analysis by first-principles simulation. The shape
of the peak reflects the transition probability and density of states. In the
high-energy domain, the probability of electronic transitions is low because
the nucleus weakly affects the continuous level. Large sharp peaks tend to be
seen in the low-energy domain, and small, broad peaks tend to be seen in the
high-energy domain.[1]

Figure 1 shows the synthetic data and model. Figure 1 shows an example of the
synthetic spectrum generated assuming XANES data measured at the O-K edge
($543.1\,\mathrm{eV}$) of $\rm{Li_{5}La_{3}Ta_{2}O_{12}}$, which has a complex
structure in the low-energy domain. $\rm{Li_{5}La_{3}Ta_{2}O_{12}}$ is
expected to be used as a solid electrolyte material.[11] Figure 1 shows the
spectrum components; the XANES spectrum consists of peaks and a step structure
called the absorption edge. At least one of the peaks appears near the
absorption edge and is called the white line.

We propose a model that discriminates the peaks below and above the absorption
edge. We assume that there are $K_{1}$ peaks below the absorption edge and
$K_{2}$ peaks above the absorption edge. The spectral function
$f(E;\theta,\mathcal{K})$ for the incident energy value of $E\in\mathbb{R}$
can be written as

$\displaystyle f(E;\theta,\mathcal{K})$ $\displaystyle:=$ $\displaystyle
f_{step}(E;\theta_{step})+f_{peak}(E;\theta_{peak},\mathcal{K}),$ (1)

where the first term represents the absorption edge and white line. The second
term represents the $K_{1}+K_{2}$ peaks, where $\mathcal{K}:=(K_{1},K_{2})$ is
the parameter of the number of peaks.
$\theta=\left\\{\theta_{step},\theta_{peak}\right\\}$ is the spectral
parameter. We use the arctangent function for the model of the step structure.
The shape of peaks depends on the instrumental resolution. If the instrumental
resolution is smaller than the intrinsic lifetime broadenings, peaks have the
Lorentzian lineshape. On the other hand, if the instrumental resolution
dominates, peaks have a Gaussian lineshape.

Our framework is envisioned to be applied most recently to XANES data measured
at the O-K edge in Saga Light Source (SAGA-LS), and we use the Gaussian for
the model of peaks because the instrumental resolution dominates. The
following discussion does not lose its generality when the Gaussian is
replaced by the Lorentzian or pseudo-Voigt function.

The function $f_{step}(E;\theta_{step})$ is defined as

$\displaystyle f_{step}(E;\theta_{step})$ $\displaystyle:=$ $\displaystyle
H\left\\{\frac{1}{2}+\frac{1}{\pi}\arctan\left(\frac{E-E_{0}}{\Gamma/2}\right)\right\\}$
(2) $\displaystyle+A\exp\left\\{-4\ln 2\left(\frac{E-(E_{0}+\Delta
E)}{\omega}\right)^{2}\right\\},$

where $\theta_{step}=\\{H,E_{0},\Gamma,A,\Delta E,\omega\\}$ is the parameter
set.

Hartree atomic units are used here. Here, $H$, $E_{0}$, and $\Gamma$ are the
step’s height, position, and width, respectively. $A$, $E_{0}+\Delta E$, and
$\omega$ are the white line’s height, position, and width, respectively.

The function of $K_{1}+K_{2}$ peaks is defined as

$\displaystyle f_{peak}(E;\theta_{peak},\mathcal{K}=(K_{1},K_{2}))$
$\displaystyle:=$ $\displaystyle
f_{1}(E;\theta_{1},K_{1})+f_{2}(E;\theta_{2},K_{2}),$ (3) $\displaystyle
f_{1}(E;\theta_{1},K_{1})$ $\displaystyle:=$
$\displaystyle\sum_{k_{1}=1}^{K_{1}}F_{k_{1}}\exp\left\\{-4\ln
2\left(\frac{E-E_{k_{1}}}{W_{k_{1}}}\right)^{2}\right\\},$ (4) $\displaystyle
f_{2}(E;\theta_{2},K_{2})$ $\displaystyle:=$
$\displaystyle\sum_{k_{2}=1}^{K_{2}}F_{k_{2}}\exp\left\\{-4\ln
2\left(\frac{E-E_{k_{2}}}{W_{k_{2}}}\right)^{2}\right\\},$ (5)

where $\theta_{peak}=\left\\{\theta_{1},\theta_{2}\right\\}$,
$\theta_{1}=\\{F_{k_{1}},E_{k_{1}},W_{k_{1}}\\}_{{k_{1}}=1}^{K_{1}}$, and
$\theta_{2}=\\{F_{k_{2}},E_{k_{2}},W_{k_{2}}\\}_{{k_{2}}=1}^{K_{2}}$ are
parameters and represent each peak’s height, position, and width,
respectively. The function $f_{1}(E;\theta_{1},K_{1})$ represents the $K_{1}$
peaks below the absorption edge and the function $f_{2}(E;\theta_{2},K_{2})$
represents the $K_{2}$ peaks above the absorption edge. The peak positions
must satisfy the conditions $E_{k_{1}}<E_{0}$, $E_{k_{2}}>E_{0}$.

We assume that an observed spectral intensity $I$ is obtained as the sum of
model $f(E;\theta,\mathcal{K})$ and noise $\epsilon$ as

$\displaystyle I$ $\displaystyle=$ $\displaystyle
f(E;\theta,\mathcal{K})+\epsilon.$ (6)

We assume that noise $\epsilon$ is generated from a Gaussian distribution with
zero mean and variance $1/b$. The conditional probability of observing $I$
given model $f(E;\theta,\mathcal{K})$ is then given by

$\displaystyle p(I|E,\theta,\mathcal{K},b)$ $\displaystyle=$
$\displaystyle\left(\frac{b}{2\pi}\right)^{\frac{1}{2}}\exp\left\\{-\frac{b}{2}\left[I-f(E;\theta,\mathcal{K})\right]^{2}\right\\},$
(7)

where $b$ is called precision. Assuming that we have a data set
$D:=\\{E_{i},I_{i}\\}^{N}_{i=1}$ consisting of $N$ observed data, the
conditional probability $p(D|\theta,\mathcal{K},b)$ is given by

$\displaystyle p(D|\theta,\mathcal{K},b)$ $\displaystyle=$
$\displaystyle\prod_{i=1}^{N}p(I_{i}|E_{i},\theta,\mathcal{K},b)$ (8)
$\displaystyle=$
$\displaystyle\left(\frac{b}{2\pi}\right)^{\frac{N}{2}}\exp(-Nb\mathcal{E}_{N}(\theta,\mathcal{K})),$
(9)

where $\mathcal{E}_{N}(\theta,\mathcal{K})$ is the error function of the
fitting function $f(E;\theta,\mathcal{K})$:

$\displaystyle\mathcal{E}_{N}(\theta,\mathcal{K})$ $\displaystyle=$
$\displaystyle\frac{1}{2N}\sum_{i=1}^{N}\left[I_{i}-f(E_{i};\theta,\mathcal{K})\right]^{2}.$
(10)

###  II.2 Bayesian formulation

Bayesian spectral deconvolution treats the data set $D$, spectral parameter
$\theta$, the parameter of the number of peaks $\mathcal{K}$, and the inverse
noise variance $b$ as random variables. We assume that $\mathcal{K}$ and $b$
are generated subject to the probabilities $p(\mathcal{K})$ and $p(b)$,
respectively, $\theta$ is generated subject to the probability
$p(\theta|\mathcal{K})$, and data set $D$ is generated subject to the
probability $p(D|\theta,\mathcal{K},b)$. The joint probability density
$p(D,\theta,\mathcal{K},b)$ is then given by

$\displaystyle
p(D,\theta,\mathcal{K},b)=p(D|\theta,\mathcal{K},b)p(\theta|\mathcal{K})p(\mathcal{K})p(b).$
(11)

The posterior distribution of data set $\theta$ given $D,\mathcal{K}$ and $b$
is represented by Bayes’ theorem as

$\displaystyle p(\theta|D,\mathcal{K},b)$ $\displaystyle=$
$\displaystyle\frac{p(D,\theta,\mathcal{K},b)}{\int d\theta
p(D,\theta,\mathcal{K},b)}$ (12) $\displaystyle=$
$\displaystyle\frac{p(D|\theta,\mathcal{K},b)p(\theta|\mathcal{K})}{\int
d\theta p(D|\theta,\mathcal{K},b)p(\theta|\mathcal{K})}$ (13) $\displaystyle=$
$\displaystyle\frac{1}{Z_{N}(\mathcal{K},b)}\left(\frac{b}{2\pi}\right)^{\frac{N}{2}}\exp\left[-Nb\mathcal{E}_{N}(\theta)\right]p(\theta|\mathcal{K}),$
(14) $\displaystyle Z_{N}(\mathcal{K},b)$ $\displaystyle:=$ $\displaystyle\int
d\theta\left(\frac{b}{2\pi}\right)^{\frac{N}{2}}\exp\left[-Nb\mathcal{E}_{N}(\theta)\right]p(\theta|\mathcal{K}),$
(15)

where $Z_{N}(\mathcal{K},b)$ is marginal likelihood, which is also called
model evidence. The negative logarithm of $Z_{N}(\mathcal{K},b)$ is called as
the Bayesian free energy:

$\displaystyle F_{N}(\mathcal{K},b):=-\log Z_{N}(\mathcal{K},b).$ (16)

###  II.3 Prior distributions

The probability density of the spectral parameter $\theta$ given the parameter
$\mathcal{K}$ is called the prior distribution. We show the discrimination of
the peaks below and above the absorption edge using a prior distribution. We
also show that, depending on the prior distribution, our model includes a
conventional model that treats all peaks equally.[7, 9, 6]

In the following, we define the probability density functions of prior
distributions of each parameter. In this paper, $\text{U}(X;\alpha,\beta)$
denotes the uniform distribution defined in $\alpha\leq x\leq\beta$,
$\text{N}(X;\mu,\sigma^{2})$ denotes the normal distribution with the mean
$\mu$ and the covariance $\sigma$, and $\text{G}(X;\kappa,\vartheta)$ denotes
the gamma distribution with the shape parameter $\kappa$ and the scale
parameter $\vartheta$ as

$\displaystyle\text{U}(X;\alpha,\beta)$ $\displaystyle=$
$\displaystyle\frac{1}{\beta-\alpha}\,\,\,\,\,\,\,(\alpha<X<\beta),$ (17)
$\displaystyle\text{N}(X;\mu,\sigma^{2})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left\\{-\frac{1}{2\sigma^{2}}(X-\mu)^{2}\right\\},$
(18) $\displaystyle\text{G}(X;\kappa,\vartheta)$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma(\kappa)\vartheta^{\kappa}}X^{\kappa-1}e^{-X/\vartheta},$
(19)

where we refer to $\alpha,\beta,\mu,\sigma,\kappa$, and $\vartheta$ as
hyperparameters.

####  II.3.1 Proposed model

We define the relationship between the position of the absorption edge
$E_{0}\in\theta_{step}$ and those of peaks
$E_{k_{1}},E_{k_{2}}\in\theta_{peak}$ as

$\displaystyle\Delta E_{k_{1}}$ $\displaystyle:=$ $\displaystyle
E_{k_{1}}-E_{0},$ (20) $\displaystyle\Delta E_{k_{2}}$ $\displaystyle:=$
$\displaystyle E_{k_{2}}-E_{0},$ (21)

where we can rewrite the conditions $E_{k_{1}}<E_{0}$ and $E_{k_{2}}>E_{0}$ as
$\Delta E_{k_{1}}<0$ and $\Delta E_{k_{2}}>0$, respectively. The parameters
$\theta_{step}$ and $\theta_{peak}$ are not independent variables because
$E_{k_{1}}\in\theta_{peak}$ and $E_{k_{2}}\in\theta_{peak}$ depend on
$E_{0}\in\theta_{step}$ via Eqs. (20) and (21). To treat parameters
independently, we reparametrize $\theta_{peak},\theta_{1},$ and $\theta_{2}$
as
$\theta^{\prime}_{peak}:=\left\\{\theta^{\prime}_{1},\theta^{\prime}_{2}\right\\}$,
$\theta^{\prime}_{1}:=\\{F_{k_{1}},\Delta
E_{k_{1}},W_{k_{1}}\\}_{{k_{1}}=1}^{K_{1}}$, and
$\theta^{\prime}_{2}:=\\{F_{k_{2}},\Delta
E_{k_{2}},W_{k_{2}}\\}_{{k_{2}}=1}^{K_{2}}$. Then, the prior distribution
$p(\theta|\mathcal{K})$ can be rewritten as

$\displaystyle p(\theta|\mathcal{K})$
$\displaystyle=p(\theta_{step},\theta_{peak}|\mathcal{K}),$ (24)
$\displaystyle=p(\theta_{peak}|\theta_{step},\mathcal{K})p(\theta_{step}|\mathcal{K}),$
$\displaystyle=p(\theta^{\prime}_{peak}|\mathcal{K})p(\theta_{step}),$

where $\theta_{step}$ and $\theta^{\prime}_{peak}$ are independent variables.

Table 1: Prior distributions for absorption edge and white line. spectral parameter | prior | hyperparameter |
---|---|---|---
$H$ | $p(H):=\text{U}(H;\alpha_{H},\beta_{H})$ | $\alpha_{H},\,\,\beta_{H}$ |
$E_{0}$ | $p(E_{0}):=\text{N}(E_{0};\mu_{E_{0}},\sigma_{E_{0}}^{2})$ | $\mu_{E_{0}},\,\,\sigma_{E_{0}}$ |
$\Gamma$ | $p(\Gamma):=\text{U}(\Gamma;\alpha_{\Gamma},\beta_{\Gamma})$ | $\alpha_{\Gamma},\,\,\beta_{\Gamma}$ |
$A$ | $p(A):=\text{G}(A;\kappa_{A},\vartheta_{A})$ |  $\kappa_{A}$, $\vartheta_{A}$ |
$\Delta E$ | $p(\Delta E):=\text{N}(\Delta E;\mu_{\Delta E},\sigma_{\Delta E}^{2})$ | $\mu_{\Delta E}=0,\,\,\sigma_{\Delta E}$ |
$\omega$ | $p(\omega):=\text{U}(\omega;\alpha_{\omega},\beta_{\omega})$ | $\alpha_{\omega},\,\,\beta_{\omega}$ |
Table 2: Prior distributions for peaks in proposed model. spectral parameter | prior | hyperparameter |
---|---|---|---
$F_{k_{1}}$ | $p(F_{k_{1}}):=\text{G}(F_{k_{1}};\kappa_{F_{k_{1}}},\vartheta_{F_{k_{1}}})$ | $\kappa_{F_{k_{1}}},\,\,\vartheta_{F_{k_{1}}}$ |
$F_{k_{2}}$ | $p(F_{k_{2}}):=\text{G}(F_{k_{2}};\kappa_{F_{k_{2}}},\vartheta_{F_{k_{2}}})$ | $\kappa_{F_{k_{2}}},\,\,\vartheta_{F_{k_{2}}}$ |
$\Delta E_{k_{1}}$ | $p(\Delta E_{k_{1}}):=\text{U}(\Delta E_{k_{1}};\alpha_{\Delta E_{k_{1}}},\beta_{\Delta E_{k_{1}}})$ | $\alpha_{\Delta E_{k_{1}}},\,\,\beta_{\Delta E_{k_{1}}}$ |
$\Delta E_{k_{2}}$ | $p(\Delta E_{k_{2}}):=\text{U}(\Delta E_{k_{2}};\alpha_{\Delta E_{k_{2}}},\beta_{\Delta E_{k_{2}}})$ | $\alpha_{\Delta E_{k_{2}}},\,\,\beta_{\Delta E_{k_{2}}}$ |
$W_{k_{1}}$ | $p(W_{k_{1}}):=\text{G}(W_{k_{1}};\kappa_{W_{k_{1}}},\vartheta_{W_{k_{1}}})$ | $\kappa_{W_{k_{1}}},\,\,\vartheta_{W_{k_{1}}}$ |
$W_{k_{2}}$ | $p(W_{k_{2}}):=\text{G}(W_{k_{2}};\kappa_{W_{k_{2}}},\vartheta_{W_{k_{2}}})$ | $\kappa_{W_{k_{2}}},\,\,\vartheta_{W_{k_{2}}}$ |

First, we define the prior distributions on the absorption edge and the white
line. The prior distribution of $\theta_{step}$ can be written as the product
of the prior distribution of each parameter:

$\displaystyle p(\theta_{step})$ $\displaystyle=$ $\displaystyle
p(H)p(E_{0})p(\Gamma)p(A)p(\Delta E)p(\omega).$ (25)

We show definitions of the prior distributions
$p(H),p(E_{0}),p(\Gamma),p(A),p(\Delta E)$, and $p(\omega)$ in Table 1. Since
the position of the absorption edge $E_{0}$ appears near the ionization energy
$\mu_{E_{0}}$ of the element of interest, we set up a Gaussian distribution
centered on $\mu_{E_{0}}$ as $p(E_{0})$. Since the white line position
$E_{0}+\Delta E$ is near the absorption edge, we use a Gaussian distribution
centered on $\mu_{\Delta E}:=0$ as $p(\Delta E)$. We use the gamma
distribution for the intensity parameter of white line $A$ to prevent the
intensity becoming too small. We set the uniform distribution for $H,\Gamma$,
and $\omega$ since we have no knowledge except for a rough range.

Next, we define the prior distributions of peaks. We can write the prior
distribution of $\theta^{\prime}_{peak}$ as

$\displaystyle p(\theta^{\prime}_{peak})$ $\displaystyle=$ $\displaystyle
p(\theta^{\prime}_{1}|K_{1})p(\theta^{\prime}_{2}|K_{2}),$ (26) $\displaystyle
p(\theta^{\prime}_{1}|K_{1})$ $\displaystyle=$
$\displaystyle\prod_{k_{1}=1}^{K_{1}}p(F_{k_{1}})p(\Delta
E_{k_{1}})p(W_{k_{1}}),$ (27) $\displaystyle p(\theta^{\prime}_{2}|K_{2})$
$\displaystyle=$ $\displaystyle\prod_{k_{2}=1}^{K_{2}}p(F_{k_{2}})p(\Delta
E_{k_{2}})p(W_{k_{2}}).$ (28)

We show definitions of the prior distributions
$p(F_{k_{1}}),p(F_{k_{2}}),p(\Delta E_{k_{1}}),p(\Delta
E_{k_{2}}),p(W_{k_{1}})$, and $p(W_{k_{2}})$ in Table 2. We restrict the
positions of peaks by the uniform distribution. Assuming that the range of the
peak positions is $E_{min}\leq E\leq E_{max}$, we set hyperparameters
$\alpha_{\Delta E_{k_{1}}},\beta_{\Delta E_{k_{1}}},\alpha_{\Delta
E_{k_{2}}}$, and $\beta_{\Delta E_{k_{2}}}$ as

$\begin{cases}\;\alpha_{\Delta E_{k_{1}}}=E_{min}-E_{0}-\sigma_{E_{0}}\\\
\;\beta_{\Delta E_{k_{1}}}=0\end{cases}\begin{cases}\;\alpha_{\Delta
E_{k_{2}}}=0\\\ \;\beta_{\Delta
E_{k_{2}}}=E_{max}-E_{0}+\sigma_{E_{0}},\end{cases}$ (29)

where we consider that $p(E_{0})$ is the normal distribution with width
$\sigma_{E_{0}}$. These prior distributions satisfy the condition $\Delta
E_{k_{1}}<0$, $\Delta E_{k_{2}}>0$. We use the gamma distribution for the
prior distribution of the intensity and width. Using the gamma distribution,
we can restrict the parameters to non zero and control the shape of the
distribution by adjusting the hyperparameters. We can incorporate the prior
distributions on the knowledge that $K_{1}$ peaks below $E_{0}$ are sharp and
$K_{2}$ peaks above $E_{0}$ are broad. The examples of concrete
hyperparameters are shown in Sect. 3.

Figure 2: Illustration of conventional and proposed models.

####  II.3.2 Conventional model

We show that our model can represent the conventional model by unifying the
prior distributions below and above the absorption edge as

$\displaystyle p(F_{k_{1}})$ $\displaystyle=$ $\displaystyle p(F_{k_{2}}),$
(30) $\displaystyle p(E_{k_{1}})$ $\displaystyle=$ $\displaystyle
p(E_{k_{2}}),$ (31) $\displaystyle p(W_{k_{1}})$ $\displaystyle=$
$\displaystyle p(W_{k_{2}}).$ (32)

Since the conventional model does not discriminate the peaks below and above
the step, $E_{0}$ is independent of $E_{k_{1}}$ and $E_{k_{2}}$. Thus,
$\theta_{step}$ and $\theta_{peak}$ are independent and the right side of Eq.
(24) can be written as

$\displaystyle
p(\theta_{step},\theta_{peak}|\mathcal{K})=p(\theta_{step})p(\theta_{peak}|\mathcal{K}).$
(33)

Here, all $K_{1}+K_{2}$ peaks have the same conditions. Thus, we do not need
to set $K_{1}$ or $K_{2}$, and only determine $K:=K_{1}+K_{2}$, which
corresponds to their sum. We define the parameter of the number of peaks as
$\mathcal{K}=K$ for the conventional model.

Table 3: Prior distributions for peaks in conventional model. spectral parameter | prior | hyperparameter |
---|---|---|---
$F_{k}$ | $p(F_{k}):=\text{U}(F_{k};\alpha_{F_{k}},\beta_{F_{k}})$ | $\alpha_{F_{k}},\,\,\beta_{F_{k}}$ |
$E_{k}$ | $p(E_{k}):=\text{U}(E_{k};\alpha_{E_{k}},\beta_{E_{k}})$ | $\alpha_{E_{k}},\,\,\beta_{E_{k}}$ |
$W_{k}$ | $p(W_{k}):=\text{U}(W_{k};\alpha_{W_{k}},\beta_{W_{k}})$ | $\alpha_{W_{k}},\,\,\beta_{W_{k}}$ |

Figure 2 shows the relationship between the proposed and conventional models.
As seen in the figure, there are $K_{1}$ peaks in $E_{min}<E<E_{0}$ and
$K_{2}$ peaks in $E_{0}<E<E_{max}$ in the proposed model and $K=K_{1}+K_{2}$
peaks in $E_{min}<E<E_{max}$ in the conventional model. We redefine
$f_{peak}(E;\theta_{peak},\mathcal{K})$ for the conventional model as

$\displaystyle
f_{peak}(E;\theta_{peak},\mathcal{K}=K):=\sum_{k=1}^{K}F_{k}\exp\left\\{-4\ln
2\left(\frac{E-E_{k}}{W_{k}}\right)^{2}\right\\}.$ (34)

The prior distributions of the peaks in the conventional model can be written
as
$p(\theta_{peak}|\mathcal{K})=\prod_{k}^{K=K_{1}+K_{2}}p(F_{k})p(E_{k})p(W_{k})$.
We show definitions of the prior distributions $p(F_{k}),p(E_{k}),$ and
$p(W_{k})$ in Table 3. As shown in Table 3, we use uniform distributions for
all parameters. We set hyperparameters for the prior distributions of the
position of peaks $E_{k}$, as $\alpha_{E_{k}}=E_{min}$ and
$\beta_{E_{k}}=E_{max}$.

Since XANES contains peaks with a wide variety of intensities and widths, the
conventional model has no choice but to use a uniform distribution. The
proposed model distinguishes the peaks below and above the absorption edge by
devising a prior distribution for the energy position. This enabled us to
reflect the physical properties in accordance with the energy domain in the
prior distribution design.

###  II.4 Exchange Monte Carlo method

The calculation of the posterior probability Eq. (14) is generally difficult
since it requires high dimensional integration as Eq. (15). However, we can
obtain an empirical posterior distribution from samples by EMC without
calculating $Z_{N}(\mathcal{K},b)$.

We prepare the different $L$ inverse variances $\\{b_{l}\\}_{l=1}^{L}$ and
posterior distributions $\\{P(\theta_{l}|D,\mathcal{K},b_{l})\\}_{l=1}^{L}$.
We can obtain samples from the following joint density by EMC:

$\displaystyle
p(\theta_{1},\theta_{2},...\theta_{L}|D,\mathcal{K},b_{1},b_{2},...b_{L})$
$\displaystyle=$
$\displaystyle\prod_{l=1}^{L}p(\theta_{l}|D,\mathcal{K},b_{l}),$ (35)

where the inverse variance satisfy the condition $0=b_{1}<b_{2}<...b_{L}$ and
are called the inverse temperatures. In addition to Metropolis-Hasting sweeps
at each $L$ Markov chain called replicas, EMC uses exchange moves between
adjacent pairs of replicas. To satisfy a detailed balance, the swap is
accepted with probability

$\displaystyle u(\theta_{l+1},\theta_{l},b_{l+1},b_{l})$ $\displaystyle:=$
$\displaystyle\min\left[1,\frac{p(\theta_{l+1}|D,\mathcal{K},b_{l})p(\theta_{l}|D,\mathcal{K},b_{l+1})}{p(\theta_{l}|D,\mathcal{K},b_{l})p(\theta_{l+1}|D,\mathcal{K},b_{l+1})}\right]$
(36) $\displaystyle=$
$\displaystyle\min\left[1,\exp\\{N(b_{l+1}-b_{l})(\mathcal{E}_{N}(\theta_{l+1})-\mathcal{E}_{N}(\theta_{l}))\\}\right],$
(37)

where the $\theta_{l}$ is the current state of each replica.

At low temperatures, the posterior density is concentrated around the
parameter $\theta$ that minimizes the error function $\mathcal{E}_{N}(\theta)$
and has complex local minima in the parameter space. As the temperature
increases, the posterior distribution asymptotes to the prior distributions,
enabling the Markov chain to explore the parameter space without becoming
trapped in the local regions. If the replicas at low temperatures are trapped
in the local minima, the exchange move enables it to reach the global minima
through exploration at higher temperatures.

Moreover, we can calculate the model evidence or the free energy from samples
obtained by EMC. Here, we define the auxiliary function: [12, 7, 9]

$\displaystyle\tilde{Z}_{N}(\mathcal{K},b)$ $\displaystyle:=$
$\displaystyle\int
d\theta\exp\left[-Nb\mathcal{E}_{N}(\theta)\right]p(\theta|\mathcal{K}).$ (38)

Equations.(14) and (15) can be rewritten as

$\displaystyle p(\theta|D,\mathcal{K},b)$ $\displaystyle=$
$\displaystyle\frac{1}{\tilde{Z}_{N}(\mathcal{K},b)}\exp\left[-Nb\mathcal{E}_{N}(\theta)\right]p(\theta|\mathcal{K}),$
(39) $\displaystyle Z_{N}(\mathcal{K},b)$ $\displaystyle=$
$\displaystyle\left(\frac{b}{2\pi}\right)^{\frac{b}{2}}\tilde{Z}_{N}(\mathcal{K},b).$
(40)

We can calculate $\tilde{Z}_{N}(\mathcal{K},b_{l})$ as

$\displaystyle\tilde{Z}_{N}(\mathcal{K},b_{l})$ $\displaystyle=$
$\displaystyle\frac{\tilde{Z}_{N}(\mathcal{K},b_{l})}{\tilde{Z}_{N}(\mathcal{K},b_{l-1})}\times\cdots\times\frac{\tilde{Z}_{N}(\mathcal{K},b_{2})}{\tilde{Z}_{N}(\mathcal{K},b_{1})}$
(41) $\displaystyle=$
$\displaystyle\prod_{l^{\prime}=1}^{l-1}\frac{\tilde{Z}_{N}(\mathcal{K},b_{l^{\prime}+1})}{\tilde{Z}_{N}(\mathcal{K},b_{l^{\prime}})}$
(42) $\displaystyle=$
$\displaystyle\prod_{l^{\prime}=1}^{l-1}\left<\exp\left[-N(b_{l^{\prime}+1}-b_{l^{\prime}})\mathcal{E}_{N}(\theta_{l}^{\prime})\right]\right>_{b_{l}^{\prime}},$
(43)

where $<...>_{b_{l}^{\prime}}$ denotes the expectation over the probability
distribution $p(\theta_{l}|D,\mathcal{K},b_{l})$. This is the extended method
of importance sampling. [13] We can then obtain Bayes free energy as

$\displaystyle
F_{N}(\mathcal{K},b_{l})=-\frac{b_{l}}{2}\log\left(\frac{b_{l}}{2\pi}\right)-\log\tilde{Z}_{N}(\mathcal{K},b_{l}).$
(44)

###  II.5 Parameter estimation

Here, we describe the method to estimate the number of peaks, noise variance,
and spectral parameters. The estimation of $\mathcal{K}$ is called model
selection, where the most probable model is chosen.[13] The optimal values of
$\mathcal{K}^{\prime}$ and $b_{l^{\prime}}$ can be obtained by the empirical
Bayes approach as

$\displaystyle(\mathcal{K}^{\prime},l^{\prime})$ $\displaystyle:=$
$\displaystyle\mathop{\rm
arg~{}max}\limits_{\mathcal{K},l}Z_{N}(\mathcal{K},b_{l})$ (45)
$\displaystyle=$ $\displaystyle\mathop{\rm
arg~{}min}\limits_{\mathcal{K},l}F_{N}(\mathcal{K},b_{l}),$ (46)

where we select the noise value from the discrete inverse temperature
$\\{b_{l}\\}_{l=1}^{L}$. We will only mention that we can interpolate the free
energy by the multihistogram method considering the continuity of noise
variance.[9]

Let $\\{\theta_{l^{\prime}}^{(m)}\\}_{m=1}^{M}$ be EMC samples for the $l$-th
replica, where $M$ is the number of EMC samples. The optimal value of spectrum
parameter $\theta_{l^{\prime}}^{(m^{\prime})}$ can be obtained by the
maximization of the posterior probability as

$\displaystyle m^{\prime}$ $\displaystyle=$ $\displaystyle\mathop{\rm
arg~{}max}\limits_{m}p(\theta_{l^{\prime}}^{(m)}|D,\mathcal{K}^{\prime},b_{l^{\prime}}),$
(47) $\displaystyle=$ $\displaystyle\mathop{\rm
arg~{}max}\limits_{m}\left[\exp\left[-Nb_{l^{\prime}}\mathcal{E}_{N}(\theta_{l^{\prime}}^{(m)})\right]p(\theta_{l^{\prime}}^{(m)}|\mathcal{K})\right].$
(48)

This technique is called maximum a posteriori (MAP) estimation. From the
aforementioned steps, we can estimate the number of peaks
$\mathcal{K}^{\prime}$, inverse noise variance $b_{l^{\prime}}$, and spectrum
parameter $\theta_{l^{\prime}}^{(m^{\prime})}$.

Table 4: Hyperparameters of priors for absorption edge and the white line. spectral parameter | hyperparameter
---|---
$H$ | $\alpha_{H}=0.8,\,\,\beta_{H}=0.9$
$E_{0}$ | $\mu_{E_{0}}=543.1,\,\,\sigma_{E_{0}}=2.0$
$\Gamma$ | $\alpha_{\Gamma}=0.5,\,\,\beta_{\Gamma}=1.4$
$A$ |  $\kappa_{A}=2.6$, $\vartheta_{A}=0.6$
$\Delta E$ | $\mu_{\Delta E}=0.0,\,\,\sigma_{\Delta E}=2.0$
$\omega$ | $\alpha_{\omega}=2.0,\,\,\beta_{\omega}=4.0$
Table 5: Hyperparameters of priors for peaks in the conventional model. spectral parameter | hyperparameter
---|---
$F_{k}$ | $\alpha_{F_{k}}=0.0,\,\,\beta_{F_{k}}=1.4$
$E_{k}$ | $\alpha_{E_{k}}=530,\,\,\beta_{E_{k}}=590$
$W_{k}$ | $\alpha_{W_{k}}=0.5,\,\,\beta_{W_{k}}=15.0$
Table 6: Hyperparameters of priors for peaks in the proposed model. spectral parameter | hyperparameter
---|---
$F_{k_{1}}$ | $\kappa_{F_{k_{1}}}=2.6,\,\,\vartheta_{F_{k_{1}}}=0.6$
$F_{k_{2}}$ | $\kappa_{F_{k_{2}}}=4.0,\,\,\vartheta_{F_{k_{2}}}=0.1$
$\Delta E_{k_{1}}$ | $\alpha_{E_{k_{1}}}=-15.1,\,\,\beta_{E_{k_{1}}}=0.0$
$\Delta E_{k_{2}}$ | $\alpha_{E_{k_{2}}}=0.0,\,\,\beta_{E_{k_{2}}}=48.9$
$W_{k_{2}}$ | $\kappa_{W_{k_{1}}}=3.0,\,\,\vartheta_{W_{k_{1}}}=1.0$
$W_{k_{2}}$ | $\kappa_{W_{k_{2}}}=11.0,\,\,\vartheta_{W_{k_{2}}}=0.8$

###  II.6 Probability of the number of peaks

As previously described, we estimate the number of peaks and noise variance by
the empirical Bayes approach. In addition, we can calculate the probability
density of the number of peaks and noise variance by the hierarchical Bayes
approach as

$\displaystyle p(\mathcal{K},b|D)$ $\displaystyle=$ $\displaystyle\frac{\int
d\theta\,p(D,\theta,\mathcal{K},b)}{\sum_{\mathcal{K}}\int db\int
d\theta\,p(D,\theta,\mathcal{K},b)}$ (49) $\displaystyle=$
$\displaystyle\frac{\exp\left[-F_{N}(\mathcal{K},b)\right]}{\sum_{\mathcal{K}}\int
db\exp\left[-F_{N}(\mathcal{K},b)\right]},$ (50) $\displaystyle
p(\mathcal{K}|D)$ $\displaystyle=$ $\displaystyle\int dbp(\mathcal{K},b),$
(51)

where we assume that $p(\mathcal{K})$ and $p(b)$ are uniform distributions.
Note that it is also possible to restrict $p(\mathcal{K})$ or $p(b)$ on the
basis of the results of first-principles calculations or the conditions of the
experiment. As seen in Eq. (50), the maximum point of $p(\mathcal{K},b|D)$ is
the equivalent to the minimum point of the free energy $F_{N}(\mathcal{K},b)$.
We can obtain $p(K|D)$ and $p(K_{1},K_{2}|D)$ via Eq. (51) since the parameter
$\mathcal{K}$ corresponds to $(K_{1},K_{2})$ in the proposed model and $K$ in
the conventional model. Moreover, we can calculate $p(K_{1}|D)$ and
$p(K_{2}|D)$ for the proposed model by the principle of marginalization:

$\displaystyle p(K_{1}|D)$ $\displaystyle=$
$\displaystyle\sum_{K_{2}}p(K_{1},K_{2}|D),$ (52) $\displaystyle p(K_{2}|D)$
$\displaystyle=$ $\displaystyle\sum_{K_{1}}p(K_{1},K_{2}|D).$ (53)

Thus, we can calculate the probability densities of the number of peaks below
and above the absorption edge separately. As previously mentioned, the $K_{1}$
peaks are more important than the $K_{2}$ peaks from the viewpoint of physical
properties. Focusing on the crucial domain and discussing the number of peaks
are impossible using the conventional model, but the proposed model makes it
possible.

##  III Numerical Experiments

In this section, we compare the conventional model and proposed models by
numerical experiments on the synthetic data in Fig. 1. The synthetic data was
generated in accordance with Eq. (6) with the number of peaks
$K=(K_{1}+K_{2})=10$, inverse variance of noise $b=3000$, and $N=703$. In the
synthetic data, $K_{1}=5$ peaks are located below the absorption edge, whereas
$K_{2}=5$ peaks are located above. Since the absorption edge appears near the
ionization energy, $E_{0}$ is chosen to be about $543.1$ eV [14].

For the EMC method, the number $L$ of replicas was set using a geometric
progression as[15]

$b_{l}=\left\\{\begin{aligned} &\,3000\cdot\xi^{(l-L+2)}&(l\neq 1)\\\
&\,0&(l=1),\end{aligned}\right.$ (54)

where we set $b_{L-2}=3000$, which is the true value of the noise variance.
The hyperparameters of the prior distributions are shown in Tables 4-6. Here,
we define $50$ iterations of the metropolis update as one Monte Carlo step
(MCS). We performed 60,000 MCSs and discarded the first 30,000 of them as the
burn-in period.

Figure 3: Transition and autocorrelation of error function
$\mathcal{E}_{N}(\theta)$. (a) Error function at each MCS. The red and blue
lines show $\mathcal{E}_{N}(\theta)$ of the conventional and proposed models,
respectively. The shade represents the standard deviation for 10 independent
runs. (b) The red and blue lines show the autocorrelation functions of the
proposed and conventional models, respectively. The shade represents the
standard error for 10 independent runs.

Figure 4: Bayesian free energy as a function of inverse temperature $b$ of (a)
conventional model $F_{N}(K,b)$ and (b) proposed model $F_{N}(K_{1},K_{2},b)$.
The vertical black dashed lines show the true value $b=3000$.

Figure 5: Bayesian free energy and probability of number of peaks for
conventional and proposed models. (a) Free energy $F_{N}(K,b=3000)$. (b)
Probability $p(K|D)$. (c) Free energy $F_{N}(K_{1},K_{2},b=3000)$. (d)
Probability $p(K_{1},K_{2}|D)$.

First, we compare the proposed and conventional models in terms of convergence
speed and sampling efficiency. We set the number of peaks as $K=10$ for the
conventional model and $(K_{1},K_{2})=(5,5)$ for the proposed model. The
number of replicas and inverse temperatures were set as $\xi=1.18$ and $L=92$
for both models, respectively. Figure 3 shows the error function
$\mathcal{E}_{N}(\theta)$ against MCSs. The initial value of each state
$\theta_{l}$ for the EMC method was randomly chosen from the prior density
$p(\theta|\mathcal{K})$. As seen from these results, the simulation using the
proposed model can converge faster than that using the conventional model.
Figure 3 shows the autocorrelation function[16] of $\mathcal{E}_{N}(\theta)$.
We can see that the simulation of the proposed model has a lower
autocorrelation than that of the conventional model. From these results, we
can say that the simulation using the proposed model can converge faster and
obtain low-correlation samples in a small iteration compared with the
conventional model.

In the following, we show the results of the estimation of noise variance,
number of peaks, and spectral parameters. For inverse temperatures, we set
$\xi=1.18$ and $L=92$ for the proposed model and $\xi=1.132$ and $L=120$ for
the conventional model. Figure 4 shows the free energies $F_{N}(K,b)$ and
$F_{N}(K_{1},K_{2},b)$ calculated as described in Sect. 2.4. As seen from
these results, the minimum point of the free energy as a function of $b$ is
$\hat{b}=3000$. According to these results, the noise estimation via Eq.
(2.10) seems to obtain the correct noise variance.

Figure 6: Probability of the number of peaks below and above the absorption
edge. (a)$p(K_{1}|D)$ and (b) $p(K_{2}|D)$.

Figure 7: True synthetic spectra and results of MAP fitting. (a) True spectral
model and synthetic data. (b) Fitting result of conventional model for the
estimated number of $K=12$ peaks. (c) Fitting result of proposed model for the
estimated number of $(K_{1},K_{2})=(5,5)$ peaks.

Figure 5 shows the results of the model selection. Figures 5 and 5 show the
results of $F_{N}(K,b=3000)$ and $p(K|D)$, respectively. From these results,
the estimation by the conventional method yields an incorrect number $K$ of
peaks, $K=12$. Figures 5 and 5 show the results of $F_{N}(K_{1},K_{2},b=3000)$
and $p(K_{1},K_{2}|D)$ respectively. From these results, the estimation by the
proposed method yields an appropriate number $(K_{1},K_{2})$ of peaks,
$(K_{1},K_{2})=(5,5)$. As seen from Fig. 5, both the free energy
$F_{N}(K_{1},K_{2},b=3000)$ and the probability density $p(K_{1},K_{2}|D)$
have similar values between the model of $(K_{1},K_{2})=(5,5)$ and
$(K_{1},K_{2})=(5,4)$.

For a detailed discussion, we show the probability densities $p(K_{1}|D)$ and
$p(K_{2}|D)$ in Fig. 6. Although the probability $p(K_{1}=5|D)$ has a
significantly large value, the difference between $p(K_{2}=5|D)$ and
$p(K_{2}=4|D)$ is small. Since the low-energy domain is the primary
consideration in correspondence with first-principles calculations, it is an
advantage of the proposed model to obtain more significant results for
$K_{1}$.

Figure 7 shows the results of the MAP fitting with the estimated number of
peaks, $K=12$ and $(K_{1},K_{2})=(5,5)$. As seen from Fig. 7, the result of
the conventional model is inappropriate although the peak positions in the
low-energy region seem to be estimated correctly. However, the proposed model
can estimate the overall structure of the true model correctly, as seen from
Fig. 7. Thus, we can estimate the values of the number of peaks and the
position, intensity, and width of each peak by using the proposed model.

In addition, we see that the values of free energy of the proposed model are
lower than those of the conventional model. In other words, the proposed model
has higher model evidence than the conventional model.

##  IV Discussion and Conclusion

In this paper, we proposed a model that discriminates between the low- and
high-energy domains. Through numerical experiments, we obtained the following
observations. First, the proposed model is superior to the conventional one in
terms of convergence speed and sampling efficiency. Second, the proposed model
showed better accuracy than the conventional one in the model selection and
the MAP estimation. Third, the proposed model has a lower free energy than the
conventional one. Fourth, the proposed model enables us to estimate the number
of peaks focusing on the important energy domain in terms of material science.

In this paper, we validated our framework via synthetic data analysis. Future
work will consist of real data analysis and comparison of the analysis results
with first-principles calculations. Our method will be used for this as well
and will contribute to materials science.

Futhermore, superior performances of the proposed model result from the design
of prior distributions. We also showed the importance of using physical
knowledge to design prior distributions. Our framework is applicable to not
only the analyses of XANES but those of any other spectra with energy-
dependent structures.

###### Acknowledgements.

The authors would like to express their sincere gratitude to Koki Okajima for
his insightful comment. Their sincere thanks also go to Koki Obinata for
valuable advice and discussion. This work was supported by JST CREST (Grant
Nos. JPMJCR1761 and JPMJCR1861) from the Japan Science and Technology Agency
(JST).

## References

  * [1] J. Stohr, NEXAFS Spectroscopy (Springer-Verlag, Heidelberg, 2003) 2nd ed., Vol. 25.
  * [2] F. de Groot and A. Kotani, in Core Level Spectroscopy of Solids, eds. D. D. Sarma, G. Kotliar, and Y. Tokura (CRC Press, London, 2008) Vol. 6.
  * [3] M. Kotobuki and M. Nogami, Phys. Chem. Chem. Phys. 14, 1008 (2012).
  * [4] Y. Chen, C. Chen, C. Zheng, S. Dwaraknath, M. K. Horton, J. Cabana, J. Rehr, J. Vinson, A. Dozier, J. J. Kas, K. A. Persson, and S. P. Ong, Sci. Data 8, 153 (2021).
  * [5] G. H. Major, N. Fairley, P. M. Sherwood, M. R. Linford, J. Terry, V. Fernandez, and K. Artyushkova, J. Vac. Sci. Technol. A, 38, 061203 (2021).
  * [6] K. Iwamitsu, T. Yokota, K. Murata, M. Kamezaki, M. Mizumaki, T. Uruga, and I. Akai, Phys. Status Solidi B 257, 2000107 (2020).
  * [7] K. Nagata, S. Sugita, and M. Okada, Neural Networks 28, 82 (2012).
  * [8] K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 (1996).
  * [9] S. Tokuda, K. Nagata, and M. Okada, J. Phys. Soc. Jpn. 86, 024001 (2017).
  * [10] B. K. Teo, EXAFS : Basic Principles and Data Analysis (Springer, Berlin, 1986).
  * [11] R. Jalem, M. Nakayama, W. Manalastas, J. A. Kilner, R. W. Grimes, T. Kasuga, and K. Kanamura, J. Phys. Chem. C. 119, 20783 (2015).
  * [12] Y. Ogata, Ann. Inst. Stat. Math. 42, 403 (1990).
  * [13] C. M. Bishop, Pattern Recognition and Machine Learning (Springer, New York, 2006).
  * [14] B. L. Henke, E. M. Gullikson, and J. C. Davis, At. Data Nucl. Data Tables. 54, 181 (1993).
  * [15] K. Nagata and S. Watanabe, Neural Networks 21, 980 (2008).
  * [16] M. Newman and G. Barkema, Monte Carlo Methods in Statistical Physics (Oxford University Press, Oxford, U.K., 1999).


11institutetext: Department of Mathematics, University of Wisconsin–Madison,
Madison, WI 53706, USA 11email<EMAIL_ADDRESS>22institutetext:
Curtin Institute for Computation & School of Electrical Engineering, Computing
and Mathematical Sciences, Curtin University, P.O. Box U1987, Perth, WA 6845,
Australia 22email<EMAIL_ADDRESS>
# Outlier removal for isogeometric spectral approximation with the optimally-
blended quadratures ††thanks: This work and visit of Quanling Deng in Krakow
was partially supported by National Science Centre, Poland grant no.
017/26/M/ST1/00281. This publication was made possible in part by the CSIRO
Professorial Chair in Computational Geoscience at Curtin University and the
Deep Earth Imaging Enterprise Future Science Platforms of the Commonwealth
Scientific Industrial Research Organisation, CSIRO, of Australia. Additional
support was provided by the European Union’s Horizon 2020 Research and
Innovation Program of the Marie Sklodowska-Curie grant agreement No. 777778
and the Curtin Institute for Computation.

Quanling Deng 11 0000-0002-6159-1233 Victor M. Calo 22 0000-0002-1805-4045

###### Abstract

It is well-known that outliers appear in the high-frequency region in the
approximate spectrum of isogeometric analysis of the second-order elliptic
operator. Recently, the outliers have been eliminated by a boundary penalty
technique. The essential idea is to impose extra conditions arising from the
differential equation at the domain boundary. In this paper, we extend the
idea to remove outliers in the superconvergent approximate spectrum of
isogeometric analysis with optimally-blended quadrature rules. We show
numerically that the eigenvalue errors are of superconvergence rate $h^{2p+2}$
and the overall spectrum is outlier-free. The condition number and stiffness
of the resulting algebraic system are reduced significantly. Various numerical
examples demonstrate the performance of the proposed method.

###### Keywords:

isogeometric analysis boundary penalty spectrum eigenvalue superconvergence
optimally-blended quadrature

##  1 Introduction

Isogeometric analysis (IGA) is a widely-used analysis tool that combines the
classical finite element analysis with computer-aided design and analysis
tools. It was introduced in 2005 [17, 6]. There is a rich literature since its
first development; see an overview paper [21] and the references therein. In
particular, a rich literature on IGA has been shown that the method
outperforms the classical finite element method (FEM) on the spectral
approximations of the second-order elliptic operators. The initial work [7]
showed that the spectral errors of IGA were significantly smaller when
compared with FEM approximations. In [18], the authors further explored the
advantages of IGA on the spectral approximations.

To further reduce the spectral errors on IGA spectral approximations, on one
hand, the recent work [5, 24] introduced Gauss-Legendre and Gauss-Lobatto
optimally-blended quadrature rules. By invoking the dispersion analysis which
was unified with the spectral analysis in [19], the spectral errors were shown
to be superconvergent with two extra orders. The work [11] generalized the
blended rules to arbitrary $p$-th order IGA with maximal continuity. Along the
line, the work [10] further studied the computational efficiency and the work
[1, 4, 23, 15, 14] studied its applications.

On the other hand, the spectral errors in the highest-frequency regions are
much larger than in the lower-frequency region. There is a thin layer in the
highest-frequency region which is referred to as “outliers”. The outliers in
isogeometric spectral approximations were first observed in [7] in 2006. The
question of how to efficiently remove the outliers remained open until
recently. In [12], the authors removed the outliers by a boundary penalty
technique. The main idea is to impose higher-order consistency conditions on
the boundaries to the isogeometric spectral approximations. The work proposed
to impose these conditions weakly. Outliers are eliminated and the condition
numbers of the systems are reduced significantly.

In this paper, we propose to further reduce spectral errors by combining the
optimally-blended quadrature rules and the boundary penalty technique. To
illustrate the idea, we focus on tensor-product meshes on rectangular domains.
We first develop the method in 1D and obtain the generalized matrix eigenvalue
problems. We then apply the tensor-product structure to generate the matrix
problems in multiple dimensions. By using the optimally-blended quadrature
rules, we retain the eigenvalue superconvergence rate. By applying the
boundary penalty technique, we remove the outliers in the superconvergent
spectrum. The method also reduces the condition numbers.

The rest of this paper is organized as follows. Section 2 presents the problem
and its discretization by the standard isogeometric analysis. Section 3
concerns the Gauss-Legendre and Gauss-Lobatto quadrature rules. We then
present the optimally-blended rules. In Section 4, we apply the boundary
penalty technique developed in [12] to the IGA setting with blending rules.
Section 5 collects numerical results that demonstrate the performance of the
proposed method. In particular, we perform the numerical study on the
condition numbers. Concluding remarks are presented in Section 6.

##  2 Problem setting

Let $\Omega=[0,1]^{d}\subset\mathbb{R}^{d},d=1,2,3$ be a bounded open domain
with Lipschitz boundary $\partial\Omega$. We use the standard notation for the
Hilbert and Sobolev spaces. For a measurable subset $S\subseteq\Omega$, we
denote by $(\cdot,\cdot)_{S}$ and $\|\cdot\|_{S}$ the $L^{2}$-inner product
and its norm, respectively. We omit the subscripts when it is clear in the
context. For any integer $m\geq 1$, we denote the $H^{m}$-norm and
$H^{m}$-seminorm as $\|\cdot\|_{H^{m}(S)}$ and $|\cdot|_{H^{m}(S)}$,
respectively. In particular, we denote by $H^{1}_{0}(\Omega)$ the Sobolev
space with functions in $H^{1}(\Omega)$ that are vanishing at the boundary. We
consider the classical second-order elliptic eigenvalue problem: Find the
eigenpairs $(\lambda,u)\in\mathbb{R}^{+}\times H^{1}_{0}(\Omega)$ with
$\|u\|_{\Omega}=1$ such that

$\displaystyle-\Delta u$ $\displaystyle=\lambda u\quad$
$\displaystyle\text{in}\quad\Omega,$ (1) $\displaystyle u$
$\displaystyle=0\quad$ $\displaystyle\text{on}\quad\partial\Omega,$

where $\Delta=\nabla^{2}$ is the Laplacian. The variational formulation of (1)
is to find $\lambda\in\mathbb{R}^{+}$ and $u\in H^{1}_{0}(\Omega)$ with
$\|u\|_{\Omega}=1$ such that

$a(w,u)=\lambda b(w,u),\quad\forall\ w\in H^{1}_{0}(\Omega),$ (2)

where the bilinear forms

$a(v,w):=(\nabla v,\nabla w)_{\Omega},\qquad b(v,w):=(v,w)_{\Omega}.$ (3)

It is well-known that the eigenvalue problem (2) has a countable set of
positive eigenvalues (see, for example, [2, Sec. 9.8])

$0<\lambda_{1}<\lambda_{2}\leq\lambda_{3}\leq\cdots$

and an associated set of orthonormal eigenfunctions
$\\{u_{j}\\}_{j=1}^{\infty}$, that is, $(u_{j},u_{k})=\delta_{jk},$ where
$\delta_{jk}=1$ is the Kronecker delta. Consequently, the eigenfunctions are
also orthogonal in the energy inner product since there holds
$a(u_{j},u_{k})=\lambda_{j}b(u_{j},u_{k})=\lambda_{j}\delta_{jk}.$

At the discretize level, we first discretize the domain $\Omega$ with a
uniform tensor-product mesh. We denote a general element as $\tau$ and its
collection as $\mathcal{T}_{h}$ such that
$\overline{\Omega}=\cup_{\tau\in\mathcal{T}_{h}}\tau$. Let
$h=\max_{\tau\in\mathcal{T}_{h}}\text{diameter}(\tau)$. In the IGA setting,
for simplicity, we use the B-splines. The B-spline basis functions in 1D are
defined by using the Cox-de Boor recursion formula; we refer to [8, 22] for
details. Let $X=\\{x_{0},x_{1},\cdots,x_{m}\\}$ be a knot vector with knots
$x_{j}$, that is, a nondecreasing sequence of real numbers. The $j$-th
B-spline basis function of degree $p$, denoted as $\phi^{j}_{p}(x)$, is
defined recursively as

$\displaystyle\phi^{j}_{0}(x)$ $\displaystyle=\begin{cases}1,\quad\text{if}\
x_{j}\leq x<x_{j+1},\\\ 0,\quad\text{otherwise},\\\ \end{cases}$ (4)
$\displaystyle\phi^{j}_{p}(x)$
$\displaystyle=\frac{x-x_{j}}{x_{j+p}-x_{j}}\phi^{j}_{p-1}(x)+\frac{x_{j+p+1}-x}{x_{j+p+1}-x_{j+1}}\phi^{j+1}_{p-1}(x).$

A tensor-product of these 1D B-splines produces the B-spline basis functions
in multiple dimensions. We define the multi-dimensional approximation space as
$V^{h}_{p}\subset H^{1}_{0}(\Omega)$ with (see [3, 16] for details):

$V^{h}_{p}=\text{span}\\{\phi_{j}^{p}\\}_{j=1}^{N_{h}}=\begin{cases}\text{span}\\{\phi^{j_{x}}_{p_{x}}(x)\\}_{j_{x}=1}^{N_{x}},&\text{in
1D},\\\
\text{span}\\{\phi^{j_{x}}_{p_{x}}(x)\phi^{j_{y}}_{p_{y}}(y)\\}_{j_{x},j_{y}=1}^{N_{x},N_{y}},&\text{in
2D},\\\
\text{span}\\{\phi^{j_{x}}_{p_{x}}(x)\phi^{j_{y}}_{p_{y}}(y)\phi^{j_{z}}_{p_{z}}(z)\\}_{j_{x},j_{y},j_{z}=1}^{N_{x},N_{y},N_{z}},&\text{in
3D},\\\ \end{cases}$

where $p_{x},p_{y},p_{z}$ specify the approximation order in each dimension.
$N_{x},N_{y},N_{z}$ is the total number of basis functions in each dimension
and $N_{h}$ is the total number of degrees of freedom. The isogeometric
analysis of (1) in variational formulation seeks $\lambda^{h}\in\mathbb{R}$
and $u^{h}\in V^{h}_{p}$ with $\|u^{h}\|_{\Omega}=1$ such that

$a(w^{h},u^{h})=\lambda^{h}b(w^{h},u^{h}),\quad\forall\ w^{h}\in V^{h}_{p}.$
(5)

At the algebraic level, we approximate the eigenfunctions as a linear
combination of the B-spline basis functions and substitute all the B-spline
basis functions for $w^{h}$ in (5). This leads to the generalized matrix
eigenvalue problem

$\mathbf{K}\mathbf{U}=\lambda^{h}\mathbf{M}\mathbf{U},$ (6)

where
$\mathbf{K}_{kl}=a(\phi_{p}^{k},\phi_{p}^{l}),\mathbf{M}_{kl}=b(\phi_{p}^{k},\phi_{p}^{l}),$
and $\mathbf{U}$ is the corresponding representation of the eigenvector as the
coefficients of the B-spline basis functions. In practice, we evaluate the
integrals involved in the bilinear forms $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$
numerically using quadrature rules. In the next section, we present the
Gauss–Legendre and Gauss–Lobatto quadrature rules, followed by their
optimally-blended rules.

##  3 Quadrature rules and optimal blending

In this section, we first present the classic Gauss-type quadrature rules and
then present the optimally-blended rules developed recently in [5]. While the
optimally-blended rules have been developed in [11] for arbitrary order
isogeometric elements, we focus on the lower-order cases for simplicity.

###  3.1 Gaussian quadrature rules

Gaussian quadrature rules are well-known and we present these rules by
following the book [20]. On a reference interval $\hat{\tau}$, a quadrature
rule is of the form

$\int_{\hat{\tau}}\hat{f}(\hat{\boldsymbol{x}})\
\text{d}\hat{\boldsymbol{x}}\approx\sum_{l=1}^{m}\hat{\varpi}_{l}\hat{f}(\hat{n_{l}}),$
(7)

where $\hat{\varpi}_{l}$ are the weights, $\hat{n_{l}}$ are the nodes, and $m$
is the number of quadrature points. We list the lower-order Gauss-Legendre and
Gauss-Lobatto quadrature rules in 1D below; we refer to [20] for rules with
more points. The Gauss-Legendre quadrature rules for $m=1,2,3,4$ in the
reference interval $[-1,1]$ are as follows:

$\displaystyle m=1:$ $\displaystyle\hat{n}_{1}=0,\qquad\hat{\varpi}_{1}=2;$
(8) $\displaystyle m=2:$
$\displaystyle\hat{n}_{1}=\pm\frac{\sqrt{3}}{3},\qquad\hat{\varpi}_{1}=1;$
$\displaystyle m=3:$
$\displaystyle\hat{n}_{1}=0,\quad\hat{n}_{2,3}=\pm\sqrt{\frac{3}{5}},\qquad\hat{\varpi}_{1}=\frac{8}{9},\quad\hat{\varpi}_{2,3}=\frac{5}{9};$
$\displaystyle m=4:$
$\displaystyle\hat{n}_{1,2,3,4}=\pm\sqrt{\frac{3}{7}\mp\frac{2}{7}\sqrt{\frac{6}{5}}},\qquad\hat{\varpi}_{1,2,3,4}=\frac{18\pm\sqrt{30}}{36}.$

A Gauss-Legendre quadrature rule with $m$ points integrates exactly a
polynomial of degree $2m-1$ or less. The Gauss-Lobatto quadrature rules with
$m=2,3,4,5$ in the reference interval $[-1,1]$ are as follows:

$\displaystyle m=2:$ $\displaystyle\hat{n}_{1,2}=\pm
1,\qquad\hat{\varpi}_{1}=1;$ (9) $\displaystyle m=3:$
$\displaystyle\hat{n}_{1}=0,\quad\hat{n}_{2,3}=\pm
1,\qquad\hat{\varpi}_{1}=\frac{4}{3},\quad\hat{\varpi}_{2,3}=\frac{1}{3};$
$\displaystyle m=4:$
$\displaystyle\hat{n}_{1,2}=\pm\sqrt{\frac{1}{5}},\quad\hat{n}_{3,4}=\pm
1,\quad\qquad\hat{\varpi}_{1,2}=\frac{5}{6},\quad\hat{\varpi}_{3,4}=\frac{1}{6};$
$\displaystyle m=5:$
$\displaystyle\hat{n}_{1}=0,\quad\hat{n}_{2,3}=\pm\sqrt{\frac{3}{7}},\quad\hat{n}_{4,5}=\pm
1,\quad\hat{\varpi}_{1}=\frac{32}{45},\quad\hat{\varpi}_{2,3}=\frac{49}{90},\hat{\varpi}_{4,5}=\frac{1}{10}.$

A Gauss-Lobatto quadrature rule with $m$ points integrates exactly a
polynomial of degree $2m-3$ or less.

For each element $\tau$, there is a one-to-one and onto mapping $\sigma$ such
that $\tau=\sigma(\hat{\tau})$, which leads to the correspondence between the
functions on $\tau$ and $\hat{\tau}$. Let $J_{\tau}$ be the corresponding
Jacobian of the mapping. Using the mapping, (7) induces a quadrature rule over
the element $\tau$ given by

$\int_{\tau}f(\boldsymbol{x})\
\text{d}\boldsymbol{x}\approx\sum_{l=1}^{N_{q}}\varpi_{l,\tau}f(n_{l,\tau}),$
(10)

where $\varpi_{l,\tau}=\text{det}(J_{\tau})\hat{\varpi}_{l}$ and
$n_{l,\tau}=\sigma(\hat{n}_{l})$. For blending rules, we denote by $Q$ a
general quadrature rule, by $G_{m}$ the $m-$point Gauss-Legendre quadrature
rule, by $L_{m}$ the $m-$point Gauss-Lobatto quadrature rule, by $Q_{\eta}$ a
blended rule, and by $O_{p}$ the optimally-blended rule for the $p$-th order
isogeometric analysis.

###  3.2 Optimal blending quadrature rules

Let $Q_{1}=\\{\varpi_{l,\tau}^{(1)},n_{l,\tau}^{(1)}\\}_{l=1}^{m_{1}}$ and
$Q_{2}=\\{\varpi_{l,\tau}^{(2)},n_{l,\tau}^{(2)}\\}_{l=1}^{m_{2}}$ be two
quadrature rules. We define the blended quadrature rule as

$Q_{\eta}=\eta Q_{1}+(1-\eta)Q_{2},$ (11)

where $\eta\in\mathbb{R}$ is a blending parameter. We note that the blending
parameter can be both positive and negative. The blended rule $Q_{\eta}$ for
the integration of a function $f$ is understood as

$\displaystyle\int_{\tau}f(\boldsymbol{x})\ \text{d}\boldsymbol{x}$
$\displaystyle\approx\tau\sum_{l=1}^{m_{1}}\varpi_{l,\tau}^{(1)}f(n_{l,\tau}^{(1)})+(1-\tau)\sum_{l=1}^{m_{2}}\varpi_{l,\tau}^{(2)}f(n_{l,\tau}^{(2)}).$
(12)

The blended rules for isogeometric analysis of the eigenvalue problem can
reduce significantly the spectral errors. In particular, the optimally-blended
rules deliver two extra orders of convergence for the eigenvalue errors; see
[5, 11] for the developments. We present the optimally-blended rule for the
$p$-th order isogeometric elements

$O_{p}=Q_{\eta}=\eta G_{p+1}+(1-\eta)L_{p+1},$ (13)

where the optimal blending parameters are given in the Table 1 below.

Table 1: Optimal blending parameters for isogeometric elements. $p$ | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---
$\eta$ | $\frac{1}{2}$ | $\frac{1}{3}$ | $-\frac{3}{2}$ | $-\frac{79}{5}$ | $-174$ | $-\frac{91177}{35}$ | $-\frac{105013}{2}$

We remark that there are other optimally-blended quadrature rules for
isogeometric elements developed in [5, 11]. Non-standard quadrature rules are
developed in [10] and they are shown to be equivalent with the optimally-
blended rules. They all lead to the same stiffness and mass matrices. Herein,
for simplicity, we only adopt the blended rules in the form of (13).

##  4 The boundary penalty technique

Outliers appear in the isogeometric spectral approximations when using $C^{2}$
cubic B-splines and higher-order approximations. These outliers can be removed
by a boundary penalty technique. In this section, we first recall the boundary
penalty technique introduced recently in [12]. We present the idea for 1D
problem with $\Omega=[0,1]$ and $\partial\Omega=\\{0,1\\}$. We then generalize
it at the algebraic level using the tensor-product structure for the problem
in multiple dimensions. We denote

$\alpha=\lfloor\frac{p-1}{2}\rfloor=\begin{cases}\frac{p-1}{2},&p\quad\text{is
odd},\\\ \frac{p-2}{2},&p\quad\text{is even}.\\\ \end{cases}$ (14)

The isogeometric analysis of (1) in 1D with the boundary penalty technique is:
Find $\tilde{\lambda}^{h}\in\mathbb{R}$ and $\tilde{u}^{h}\in V^{h}_{p}$ such
that

$\tilde{a}(w^{h},\tilde{u}^{h})=\tilde{\lambda}^{h}\tilde{b}(w^{h},\tilde{u}^{h}),\quad\forall\
w^{h}\in V^{h}_{p}.$ (15)

Herein, for $w,v\in V^{h}_{p}$

$\displaystyle\tilde{a}(w,v)$ $\displaystyle=\int_{0}^{1}w^{\prime}v^{\prime}\
dx+\sum_{\ell=1}^{\alpha}\eta_{a,\ell}\pi^{2}h^{6\ell-3}\Big{(}w^{(2\ell)}(0)v^{(2\ell)}(0)+w^{(2\ell)}(1)v^{(2\ell)}(1)\Big{)},$
(16a) $\displaystyle\tilde{b}(w,v)$ $\displaystyle=\int_{0}^{1}wv\
dx+\sum_{\ell=1}^{\alpha}\eta_{b,\ell}h^{6\ell-1}\Big{(}w^{(2\ell)}(0)v^{(2\ell)}(0)+w^{(2\ell)}(1)v^{(2\ell)}(1)\Big{)},$
(16b)

where $\eta_{a,\ell},\eta_{b,\ell}$ are penalty parameters set to
$\eta_{a,\ell}=\eta_{b,\ell}=1$ in default. The superscript $(2\ell)$ denotes
the $2\ell$-th derivative. We further approximate the inner-products by the
quadrature rules discussed above. For $p$-th order element and
$\tau\in\mathcal{T}_{h}$, we denote
$G_{p+1}=\\{\varpi_{l,\tau}^{G},n_{l,\tau}^{G}\\}_{l=1}^{p+1}$ and
$L_{p+1}=\\{\varpi_{l,\tau}^{L},n_{l,\tau}^{L}\\}_{l=1}^{p+1}$. Applying the
optimally-blended quadrature rules in the form of (13) to (15), we obtain the
approximated form

$\tilde{a}_{h}(w^{h},\tilde{u}^{h})=\tilde{\lambda}^{h}\tilde{b}_{h}(w^{h},\tilde{u}^{h}),\quad\forall\
w^{h}\in V^{h}_{p},$ (17)

where for $w,v\in V^{h}_{p}$

$\displaystyle\tilde{a}_{h}(w,v)$
$\displaystyle=\sum_{\tau\in\mathcal{T}_{h}}\sum_{l=1}^{p+1}\big{(}\eta\varpi_{l,\tau}^{G}\nabla
w(n_{l,\tau}^{G})\cdot\nabla
v(n_{l,\tau}^{G})+(1-\eta)\varpi_{l,\tau}^{L}\nabla
w(n_{l,\tau}^{L})\cdot\nabla v(n_{l,\tau}^{L})\big{)}$ (18)
$\displaystyle\quad+\sum_{\ell=1}^{\alpha}\eta_{a,\ell}\pi^{2}h^{6\ell-3}\Big{(}w^{(2\ell)}(0)v^{(2\ell)}(0)+w^{(2\ell)}(1)v^{(2\ell)}(1)\Big{)},$

and

$\displaystyle\tilde{b}_{h}(w,v)$
$\displaystyle=\sum_{\tau\in\mathcal{T}_{h}}\sum_{l=1}^{p+1}\big{(}\eta\varpi_{l,\tau}^{G}w(n_{l,\tau}^{G})\cdot
v(n_{l,\tau}^{G})+(1-\eta)\varpi_{l,\tau}^{L}w(n_{l,\tau}^{L})\cdot
v(n_{l,\tau}^{L})\big{)}$ (19)
$\displaystyle\quad+\sum_{\ell=1}^{\alpha}\eta_{b,\ell}h^{6\ell-1}\Big{(}w^{(2\ell)}(0)v^{(2\ell)}(0)+w^{(2\ell)}(1)v^{(2\ell)}(1)\Big{)}.$

With this in mind, we arrive at the matrix eigenvalue problem

$\tilde{\mathbf{K}}\mathbf{U}=\tilde{\lambda}^{h}\tilde{\mathbf{M}}\mathbf{U},$
(20)

where
$\tilde{\mathbf{K}}_{kl}=\tilde{\mathbf{K}}^{1D}_{kl}=\tilde{a}_{h}(\phi_{p}^{k},\phi_{p}^{l}),\tilde{\mathbf{M}}_{kl}=\tilde{\mathbf{M}}^{1D}_{kl}=\tilde{b}_{h}(\phi_{p}^{k},\phi_{p}^{l})$
in 1D. Using the tensor-product structure and introducing the outer-product
$\otimes$ (also known as the Kronecker product), the corresponding 2D matrices
(see [5] for details) are

$\displaystyle\tilde{\mathbf{K}}$
$\displaystyle=\tilde{\mathbf{K}}^{2D}=\tilde{\mathbf{K}}^{1D}_{x}\otimes\tilde{\mathbf{M}}^{1D}_{y}+\tilde{\mathbf{M}}^{1D}_{x}\otimes\tilde{\mathbf{K}}^{1D}_{y},$
(21) $\displaystyle\tilde{\mathbf{M}}$
$\displaystyle=\tilde{\mathbf{M}}^{2D}=\tilde{\mathbf{M}}^{1D}_{x}\otimes\tilde{\mathbf{M}}^{1D}_{y},$

and the 3D matrices are

$\displaystyle\tilde{\mathbf{K}}$
$\displaystyle=\tilde{\mathbf{K}}^{3D}=\tilde{\mathbf{K}}^{1D}_{x}\otimes\tilde{\mathbf{M}}^{1D}_{y}\otimes\tilde{\mathbf{M}}^{1D}_{z}+\tilde{\mathbf{M}}^{1D}_{x}\otimes\tilde{\mathbf{K}}^{1D}_{y}\otimes\tilde{\mathbf{M}}^{1D}_{z}+\tilde{\mathbf{M}}^{1D}_{x}\otimes\tilde{\mathbf{M}}^{1D}_{y}\otimes\tilde{\mathbf{K}}^{1D}_{z},$
(22) $\displaystyle\tilde{\mathbf{M}}$
$\displaystyle=\tilde{\mathbf{M}}^{3D}=\tilde{\mathbf{M}}^{1D}_{x}\otimes\tilde{\mathbf{M}}^{1D}_{y}\otimes\tilde{\mathbf{M}}^{1D}_{z},$

where $\tilde{\mathbf{K}}^{1D}_{q},\tilde{\mathbf{M}}^{1D}_{q},q=x,y,z$ are 1D
matrices generated from the modified bilinear forms in (17). We remark that in
the limiting case with uniform meshes, analytical eigenpairs to the matrix
eigenvalue problem 20 can be derived based on [9]; see [12] for IGA with Gauss
quadrature rule. We omit it here for brevity.

##  5 Numerical examples

In this section, we present numerical tests to demonstrate the performance of
the method. We consider the problem (1) with $d=1,2,3.$ The 1D problem has
true eigenpairs $\big{(}\lambda_{j}=j^{2}\pi^{2},u_{j}=\sin(j\pi
x)\big{)},j=1,2,\cdots,$ the 2D problem has true eigenpairs
$\big{(}\lambda_{jk}=(j^{2}+k^{2})\pi^{2},u_{jk}=\sin(j\pi x)\sin(k\pi
y)\big{)},j,k=1,2,\cdots,$ and the 3D problem has true eigenpairs
$\big{(}\lambda_{jkl}=(j^{2}+k^{2}+l^{2})\pi^{2},u_{jkl}=\sin(j\pi x)\sin(k\pi
y)\sin(l\pi z)\big{)},j,k,l=1,2,\cdots.$ We sort both the exact and
approximate eigenvalues in ascending order. Since outliers appear in the
spectrum for cubic ($p=3$) and higher order isogeometric elements, we focus on
$p$-th order elements with $p\geq 3$.

###  5.1 Numerical study on error convergence rates and outliers

To study the errors, we consider both the $H^{1}$-seminorm and $L^{2}$-norm
for the eigenfunctions. The optimal convergence rates in $H^{1}$-seminorm and
$L^{2}$-norm are $h^{p}$ and $h^{p+1}$ for $p$-th order elements,
respectively. For the eigenvalues, we consider the relative eigenvalue errors
defined as $\frac{|\tilde{\lambda}^{h}_{j}-\lambda_{j}|}{\lambda_{j}}.$ The
optimal convergence rate is $h^{2p}$ for $p$-th order elements. We denote by
$\rho_{p}$ the convergence rate of $p$-th order isogeometric elements.

Table 2 shows the eigenvalue and eigenfunction errors for $p=3,4,5$ in 1D
while Table 3 shows the eigenvalue errors for the problems in 2D and 3D. The
first eigenvalue error approaches the machine precision fast. Thus, we
calculate the convergence rates with coarser meshes. In all these scenarios,
the eigenfunction errors are convergent optimally, while the eigenvalue errors
are superconvergent with two extra orders, i.e., $h^{2p+2}$. These results are
in good agreement with the theoretical predictions.

Table 2: Errors and convergence rates for the first and sixth eigenpairs in 1D when using IGA with optimally-blended quadratures and the boundary penalty technique. $p$ | $N$ | $\frac{|\tilde{\lambda}^{h}_{1}-\lambda_{1}|}{\lambda_{1}}$ | $|u_{1}-\tilde{u}^{h}_{1}|_{H^{1}}$ | $\|u_{1}-\tilde{u}_{1}^{h}\|_{L^{2}}$ | $\frac{|\tilde{\lambda}^{h}_{6}-\lambda_{6}|}{\lambda_{6}}$ | $|u_{6}-\tilde{u}^{h}_{6}|_{H^{1}}$ | $\|u_{6}-\tilde{u}^{h}_{6}\|_{L^{2}}$
---|---|---|---|---|---|---|---
| 5 | 3.52e-7 | 4.95e-3 | 1.67e-4 | 3.05e-1 | 1.83e1 | 9.71e-1
| 10 | 1.32e-9 | 5.75e-4 | 9.25e-6 | 3.21e-3 | 1.50 | 3.38e-2
3 | 20 | 5.09e-12 | 7.05e-5 | 5.60e-7 | 9.57e-6 | 1.12e-1 | 1.01e-3
| 40 | 1.12e-13 | 8.77e-6 | 3.47e-8 | 3.45e-8 | 1.20e-2 | 4.94e-5
| $\rho_{3}$ | 8.04 | 3.05 | 4.08 | 7.76 | 3.54 | 4.79
| 5 | 6.90e-9 | 5.26e-4 | 1.80e-5 | 3.05e-1 | 1.87e1 | 9.91e-1
4 | 10 | 6.31e-12 | 2.91e-5 | 4.72e-7 | 7.59e-4 | 6.38e-1 | 1.45e-2
| 20 | 5.75e-14 | 1.76e-6 | 1.41e-8 | 4.42e-7 | 1.91e-2 | 1.75e-4
| $\rho_{4}$ | 10.09 | 4.11 | 5.16 | 9.70 | 4.97 | 6.23
| 5 | 1.13e-10 | 5.65e-5 | 1.97e-6 | 3.06e-1 | 1.89e1 | 1.00
5 | 10 | 1.15e-14 | 1.48e-6 | 2.43e-8 | 1.41e-4 | 2.63e-1 | 6.17e-3
| 20 | 1.27e-14 | 4.42e-8 | 3.54e-10 | 1.72e-8 | 3.30e-3 | 3.06e-5
| $\rho_{5}$ | 13.26 | 5.16 | 6.22 | 12.04 | 6.24 | 7.50
Table 3: Errors for the first and sixth eigenvalues in 2D and 3D when using IGA with optimally-blended quadratures and the boundary penalty technique. | $p$ | 3 | 4 | 5
---|---|---|---|---
$d$ | $N$ | $\frac{|\tilde{\lambda}^{h}_{1}-\lambda_{1}|}{\lambda_{1}}$ | $\frac{|\tilde{\lambda}^{h}_{6}-\lambda_{6}|}{\lambda_{6}}$ | $\frac{|\tilde{\lambda}^{h}_{1}-\lambda_{1}|}{\lambda_{1}}$ | $\frac{|\tilde{\lambda}^{h}_{6}-\lambda_{6}|}{\lambda_{6}}$ | $\frac{|\tilde{\lambda}^{h}_{1}-\lambda_{1}|}{\lambda_{1}}$ | $\frac{|\tilde{\lambda}^{h}_{6}-\lambda_{6}|}{\lambda_{6}}$
| 3 | 2.28e-5 | 1.63e-4 | 1.32e-6 | 8.52e-3 | 6.55e-8 | 1.30e-7
| 6 | 8.05e-8 | 6.10e-4 | 1.09e-9 | 9.03e-5 | 1.29e-11 | 1.15e-5
2 | 12 | 3.07e-10 | 1.94e-6 | 1.07e-12 | 6.06e-8 | 5.21e-14 | 1.61e-9
| 24 | 1.31e-12 | 7.15e-9 | 1.35e-14 | 5.37e-11 | 9.45e-14 | 2.86e-13
| $\rho_{p}$ | 8.02 | 8.19 | 10.12 | 10.54 | 7.95 | 12.8
| 2 | 6.78e-4 | 3.33e-1 | 1.00e-4 | 3.33e-1 | 1.28e-5 | 3.33e-1
| 4 | 2.15e-6 | 6.03e-4 | 6.74e-8 | 8.92e-5 | 1.79e-9 | 1.14e-5
3 | 8 | 7.94e-9 | 1.92e-6 | 5.96e-11 | 5.99e-8 | 5.30e-13 | 1.59e-9
| 16 | 3.06e-11 | 7.07e-9 | 3.00e-15 | 5.30e-11 | 4.66e-15 | 3.66e-13
| $\rho_{p}$ | 8.13 | 8.48 | 11.5 | 10.82 | 12.26 | 13.19

Figures 1–3 show the overall spectral errors when using the standard IGA and
IGA with optimally-blended rules and the boundary penalty technique. The
polynomial degrees are $p\in\\{3,4,5\\}$. There are 100 elements in 1D,
$20\times 20$ elements in 2D, and $20\times 20\times 20$ elements in 3D,
respectively. In all the scenarios, we observe that there are outliers in the
IGA spectra. These outliers are eliminated by the proposed method. Moreover,
we observe that these spectral errors are reduced significantly, especially in
the high-frequency regions.

Figure 1: Outliers in IGA spectra and their eliminations when using the IGA
with optimally-blended rules and the boundary penalty technique. There are 100
elements in 1D with polynomial degrees $p\in\\{3,4,5\\}$.

Figure 2: Outliers in IGA spectra and their eliminations when using the IGA
with optimally-blended rules and the boundary penalty technique. There are
$20\times 20$ elements in 2D with polynomial degrees $p\in\\{3,4,5\\}$.

Figure 3: Outliers in IGA spectra and their eliminations when using the IGA
with optimally-blended rules and the boundary penalty technique. There are
$20\times 20\times 20$ elements in 3D with polynomial degrees
$p\in\\{3,4,5\\}$.

###  5.2 Numerical study on condition numbers

We study the condition numbers to show further advantages of the method. Since
the stiffness and mass matrices are symmetric, the condition numbers of the
generalized matrix eigenvalue problems (6) and (20) are given by

$\gamma:=\frac{\lambda^{h}_{\max}}{\lambda^{h}_{\min}},\qquad\tilde{\gamma}:=\frac{\tilde{\lambda}^{h}_{\max}}{\tilde{\lambda}^{h}_{\min}},$
(23)

where $\lambda^{h}_{\max},\tilde{\lambda}^{h}_{\max}$ are the largest
eigenvalues and $\lambda^{h}_{\min},\tilde{\lambda}^{h}_{\min}$ are the
smallest eigenvalues of IGA and the proposed method, respectively. The
condition number characterizes the stiffness of the system. We follow the
recent work of soft-finite element method (SoftFEM) [13] and define the
condition number reduction ratio of the method with respect to IGA as

$\rho:=\frac{\gamma}{\tilde{\gamma}}=\frac{\lambda^{h}_{\max}}{\tilde{\lambda}^{h}_{\max}}\cdot\frac{\tilde{\lambda}^{h}_{\min}}{\lambda^{h}_{\min}}.$
(24)

In general, one has $\lambda^{h}_{\min}\approx\tilde{\lambda}^{h}_{\min}$ for
IGA and the proposed method with sufficient number of elements (in practice,
these methods with only a few elements already lead to good approximations to
the smallest eigenvalues). Thus, the condition number reduction ratio is
mainly characterized by the ratio of the largest eigenvalues. Finally, we
define the condition number reduction percentage as

$\varrho=100\frac{\gamma-\tilde{\gamma}}{\gamma}\,\%=100(1-\rho^{-1})\,\%.$
(25)

Table 4 shows the smallest and largest eigenvalues, condition numbers and
their reduction ratios and percentages for 1D, 2D, and 3D problems. We observe
that the condition numbers of the proposed method are significantly smaller.
For higher-order isogeometric elements, there are more outliers and these
outliers pollute a larger high-frequency region. Consequently, they lead to
larger errors in the high-frequency region. The proposed method reduces more
errors for higher-order elements. This leads to small errors in the high-
frequency region also for high-order elements. The condition number of the
proposed method reduces by about 32% for $C^{2}$ cubic, 60% for $C^{3}$
quartic, and 75% for $C^{4}$ quintic elements. This holds valid for both 2D
and 3D problems.

Table 4: Minimal and maximal eigenvalues, condition numbers, reduction ratios and percentages when using IGA and IGA with optimally-blended quadratures and the boundary penalty technique. The polynomial degrees are $p\in\\{3,4,5\\}$. There are 100, $48\times 48$ elements, and $16\times 16\times 16$ elements in 1D, 2D, and 3D, respectively. $d$ | $p$ | $\lambda_{\min}^{h}$ | $\lambda_{\max}^{h}$ | $\tilde{\lambda}_{\max}^{h}$ | $\gamma$ | $\tilde{\gamma}$ | $\rho$ | $\varrho$
---|---|---|---|---|---|---|---|---
| 3 | 9.87 | 1.46e5 | 9.87e4 | 1.47e4 | 1.00e4 | 1.47 | 32.17%
1 | 4 | 9.87 | 2.45e5 | 9.87e4 | 2.48e4 | 1.00e4 | 2.48 | 59.69%
| 5 | 9.87 | 3.93e5 | 1.00e5 | 3.98e4 | 1.02e4 | 3.92 | 74.47%
| 3 | 1.97e1 | 6.71e4 | 4.55e4 | 3.40e3 | 2.30e3 | 1.47 | 32.17%
2 | 4 | 1.97e1 | 1.13e5 | 4.55e4 | 5.72e3 | 2.30e3 | 2.48 | 59.69%
| 5 | 1.97e1 | 1.81e5 | 4.57e4 | 9.17e3 | 2.31e3 | 3.96 | 74.77%
| 3 | 2.96e1 | 1.12e4 | 7.58e3 | 3.78e2 | 2.56e2 | 1.48 | 32.23%
3 | 4 | 2.96e1 | 1.88e4 | 7.58e3 | 6.36e2 | 2.56e2 | 2.48 | 59.72%
| 5 | 2.96e1 | 3.02e4 | 7.59e3 | 1.02e3 | 2.56e2 | 3.98 | 74.89%

##  6 Concluding remarks

We improve the isogeometric spectral approximations by combining the two
ideas: optimally-blended quadratures and a boundary penalty technique. As a
result, we obtained a superconvergence of rate $h^{2p+2}$ for the eigenvalue
errors and eliminated the outliers in the spectra. The technique can be also
used to improve the spectral approximations of the Neumann eigenvalue
problems. These improvements lead to a better spatial discretization for the
time-dependent partial differential equations, which in return, improve the
overall performance of numerical methods. As future work, it would be
interesting to study the method for higher-order differential operators and
nonlinear application problems.

## References

  * [1] Bartoň, M., Calo, V., Deng, Q., Puzyrev, V.: Generalization of the Pythagorean eigenvalue error theorem and its application to isogeometric analysis. In: Numerical Methods for PDEs, pp. 147–170. Springer (2018)
  * [2] Brezis, H.: Functional analysis, Sobolev spaces and partial differential equations. Universitext, Springer, New York (2011)
  * [3] Buffa, A., De Falco, C., Sangalli, G.: Isogeometric analysis: new stable elements for the Stokes equation. International Journal for Numerical Methods in Fluids (2010)
  * [4] Calo, V., Deng, Q., Puzyrev, V.: Quadrature blending for isogeometric analysis. Procedia Computer Science 108, 798–807 (2017)
  * [5] Calo, V., Deng, Q., Puzyrev, V.: Dispersion optimized quadratures for isogeometric analysis. Journal of Computational and Applied Mathematics 355, 283–300 (2019)
  * [6] Cottrell, J.A., Hughes, T.J.R., Bazilevs, Y.: Isogeometric analysis: toward integration of CAD and FEA. John Wiley & Sons (2009)
  * [7] Cottrell, J.A., Reali, A., Bazilevs, Y., Hughes, T.J.R.: Isogeometric analysis of structural vibrations. Computer methods in applied mechanics and engineering 195(41), 5257–5296 (2006)
  * [8] De Boor, C.: A practical guide to splines, vol. 27. Springer-Verlag New York (1978)
  * [9] Deng, Q.: Analytical solutions to some generalized and polynomial eigenvalue problems. Special Matrices 9(1), 240–256 (2021). https://doi.org/doi:10.1515/spma-2020-0135, https://doi.org/10.1515/spma-2020-0135
  * [10] Deng, Q., Bartoň, M., Puzyrev, V., Calo, V.: Dispersion-minimizing quadrature rules for C1 quadratic isogeometric analysis. Computer Methods in Applied Mechanics and Engineering 328, 554–564 (2018)
  * [11] Deng, Q., Calo, V.: Dispersion-minimized mass for isogeometric analysis. Computer Methods in Applied Mechanics and Engineering 341, 71–92 (2018)
  * [12] Deng, Q., Calo, V.: A boundary penalization technique to remove outliers from isogeometric analysis on tensor-product meshes. Computer Methods in Applied Mechanics and Engineering, accepted (2021)
  * [13] Deng, Q., Ern, A.: SoftFEM: revisiting the spectral finite element approximation of elliptic operators. arXiv preprint arXiv:2011.06953 (2020)
  * [14] Deng, Q., Puzyrev, V., Calo, V.: Isogeometric spectral approximation for elliptic differential operators. Journal of Computational Science (2018)
  * [15] Deng, Q., Puzyrev, V., Calo, V.: Optimal spectral approximation of 2n-order differential operators by mixed isogeometric analysis. Computer Methods in Applied Mechanics and Engineering 343, 297–313 (2019)
  * [16] Evans, J.A., Hughes, T.J.: Isogeometric divergence-conforming b-splines for the Darcy–Stokes–Brinkman equations. Mathematical Models and Methods in Applied Sciences 23(04), 671–741 (2013)
  * [17] Hughes, T.J.R., Cottrell, J.A., Bazilevs, Y.: Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer methods in applied mechanics and engineering 194(39), 4135–4195 (2005)
  * [18] Hughes, T.J.R., Evans, J.A., Reali, A.: Finite element and NURBS approximations of eigenvalue, boundary-value, and initial-value problems. Computer Methods in Applied Mechanics and Engineering 272, 290–320 (2014)
  * [19] Hughes, T.J.R., Reali, A., Sangalli, G.: Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: comparison of p-method finite elements with k-method NURBS. Computer methods in applied mechanics and engineering 197(49), 4104–4124 (2008)
  * [20] Kythe, P.K., Schäferkotter, M.R.: Handbook of computational methods for integration. CRC Press (2004)
  * [21] Nguyen, V.P., Anitescu, C., Bordas, S.P., Rabczuk, T.: Isogeometric analysis: an overview and computer implementation aspects. Mathematics and Computers in Simulation 117, 89–116 (2015)
  * [22] Piegl, L., Tiller, W.: The NURBS book. Springer Science & Business Media (1997)
  * [23] Puzyrev, V., Deng, Q., Calo, V.: Spectral approximation properties of isogeometric analysis with variable continuity. Computer Methods in Applied Mechanics and Engineering 334, 22–39 (2018)
  * [24] Puzyrev, V., Deng, Q., Calo, V.M.: Dispersion-optimized quadrature rules for isogeometric analysis: modified inner products, their dispersion properties, and optimally blended schemes. Computer Methods in Applied Mechanics and Engineering 320, 421–443 (2017)