text:
\section{Non-intrusive reduced-basis method}
\label{sec:Non-intrusive reduced-basis method}

After providing a high-level overview on \gls{ROM} in \refsec{sec:Introduction}, we describe the non-intrusive \gls{ROM} using \gls{POD} and regression more formally in the following.
We first describe the method as it is commonly reported in the literature \cite{Hesthaven2018, Guo2018, HesthavenBook2015, Yu2019}. Then in Sections \ref{sec:preprocessing}-\ref{sec:unsteadyRB}, we detail our proposed modifications.

Although the data for the non-intrusive \gls{ROM} can originate from any numerical scheme or even experimental measurements, we illustrate the purpose of a \emph{basis} using \gls{FEM}. The high-fidelity solution to a parametrized \gls{PDE} provided by an \gls{FEM} solver is typically of the form
\begin{equation}
 \label{eq:FEM}
  s (\vekt{x}; \vekt{\mu})
  = \vekt{s}^\top(\vekt{\mu}) \vekt{\phi}(\vekt{x})
  = \sum_{i=1}^{N_h} s_i (\vekt{\mu})  \phi_i(\vekt{x}) \ ,
\end{equation}

where
$ \vekt{\phi} = [ \phi_{1}(\vekt{x}) | \cdots | \phi_{N_h}(\vekt{x}) ]^\top \in \mathbb{R}^{N_h}$
is a collection of the (e.g. Lagrange) basis functions and the solution vector $\vekt{s} (\vekt{\mu}) \in \mathbb{R}^{N_h}$ fixes the coefficient values for all $N_h$ \glspl{DOF}.
The vector $\vekt{\mu} \in \mathbb{R}^{N_d}$ collects all $N_d$ input parameters.
Given a training set with $N_{tr}$ different parameter samples $\mathbb{P}_{tr} = \{ \vekt{\mu}^{(n)} \}_{1 \leq n \leq N_{tr}}$
and corresponding high-fidelity solution coefficients
$\mathbb{S}_{tr} = \{ \vekt{s}(\vekt{\mu}) \}_{\vekt{\mu} \in \mathbb{P}_{tr}}$,
the key idea in \gls{ROM} is to construct a set of $L$ \emph{reduced basis} vectors
$\matr{V} = [ \vekt{v}^{(1)} | ... | \vekt{v}^{(L)}] \in \mathbb{R}^{N_h \times L}$
such that for any solution vector $\vekt{s} (\vekt{\mu})$ an approximation $\vekt{s}_L (\vekt{\mu})$ can be constructed as a linear combination of a small number of reduced basis vectors $L \ll N_h$:
  \begin{equation}
  \label{eq:approx}
    \vekt{s} (\vekt{\mu})
    \approx {\vekt{s}}_L (\vekt{\mu})
    = \sum_{l=1}^{L} y_{l} (\vekt{\mu}) \vekt{v}^{(l)}
    = \matr{V}\vekt{y}(\vekt{\mu}) \ ,
  \end{equation}

with the corresponding \emph{reduced coefficients} $\vekt{y} (\vekt{\mu}) =  [y_l|\cdots|y_L]^\top \in \mathbb{R}^L$.
Inserting \refeqn{eq:approx} in \refeqn{eq:FEM} gives the spatially interpolated approximation:
\begin{equation}
    s(\vekt{x};\vekt{\mu}) \approx s_L(\vekt{x};\vekt{\mu}) = \vekt{y}^\top(\vekt{\mu}) \matr{V}^\top \vekt{\phi}(\vekt{x}) \ .
\end{equation}
We refer to $\matr{V}^\top \vekt{\phi}(\vekt{x})$ as the reduced basis \emph{functions}. However, throughout this work, we mostly operate on the discrete entities $\matr{V}$ and $\vekt{s}(\vekt{\mu})$ and refer to them simply as \emph{reduced basis} and \emph{solution}, respectively.

A commonly used tool to compute the reduced basis is \gls{POD} (see \refsec{sec:POD}), which generates a set of orthonormal reduced basis vectors such that $\matr{V}^\top\matr{V}=\matr{I}$. As a result, the reduced coefficients can be computed by \emph{projecting} the high-fidelity solutions onto the reduced basis:
\begin{equation}
 \label{eq:rbcoeffs}
 \vekt{y} (\vekt{\mu}) = \matr{V}^\top\vekt{s}(\vekt{\mu}) \ .
\end{equation}

The final step of the offline stage is to recover the unknown underlying mapping from parameters to the reduced coefficients $\vekt{\pi}: \vekt{\mu} \mapsto \vekt{y} (\vekt{\mu}) $.
This is done by fitting a regression model $\tilde{\vekt{\pi}}$ on the training data consisting of the input set of parameter samples $\mathbb{P}_{tr}$ and the output set of the corresponding reduced coefficients $\mathbb{Y}_{tr}=\{\vekt{y}(\vekt{\mu})\}_{\vekt{\mu} \in \mathbb{P}_{tr}}$.
 Several choices for regression models are described in Section \ref{sec:Regression models}. During the online stage, the fitted regression model can be evaluated rapidly for new parameter samples to obtain the \emph{predicted reduced coefficients} $\tilde{\vekt{y}}$. These can be transformed back to the full space to reconstruct the \emph{predicted solution} $\tilde{\vekt{s}}_L$:

 \begin{equation}
\label{eq:up-projection}
 \tilde{\vekt{s}}_L (\vekt{\mu}) = \matr{V} \tilde{\vekt{y}}(\vekt{\mu}) = \matr{V} \tilde{\vekt{\pi}}(\vekt{\mu}) \ .
\end{equation}

The recovery requires $\mathcal{O}(N_h L)$ operations due to the matrix-vector multiplication, which involves the number of \glspl{DOF} in the full-order problem $N_h$.
Depending on the application, this may, however, not even be required either if the quantity of interest is computed directly from the reduced coefficients or if subsampling of the solutions can be used.
Note that the online evaluation is completely independent of the solver used during the offline phase.

In the following, we introduce the concepts of \gls{POD} and standardization.

\subsection{Proper orthogonal decomposition}
\label{sec:POD}

Given a \emph{snapshot matrix} $ \matr{S} = [ \vekt{s} (\vekt{\mu}^{(1)}) | ... | \vekt{s} (\vekt{\mu}^{(N_{tr})}) ] \in \mathbb{R}^{N_h \times N_{tr}}$ in which the high-fidelity solutions are arranged column-wise,
\gls{POD} makes use of the \gls{SVD} to decompose the normalized snapshot matrix $\matr{S}/\sqrt{N_{tr}}$ into two orthonormal matrices $\matr{W} \in \mathbb{R}^{N_h \times N_h}$, $\matr{Z} \in \mathbb{R}^{N_{tr} \times N_{tr}}$ and a diagonal matrix $\matr{\Sigma} \in \mathbb{R}^{N_h \times N_{tr}}$ such that

\begin{equation}
 \matr{S}/\sqrt{N_{tr}} = \matr{W} \matr{\Sigma} \matr{Z}^\top \ .
\end{equation}

 Columns of $\matr{W} = [\vekt{w}_1 | \cdots | \vekt{w}_{N_h}]$ and $\matr{Z} =  [\vekt{z}_1 | \cdots | \vekt{z}_{N_{tr}}]$ are the \emph{left} and \emph{right singular vectors} of both $\matr{S}$ and $\matr{S}/\sqrt{N_{tr}}$.
The entries in the rectangular diagonal matrix $\matr{\Sigma}=diag({\sigma_1, \cdots, \sigma_N})$ are the \emph{singular values} of $\matr{S}/\sqrt{N_{tr}}$ and are ordered decreasingly $\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_{N} \ge 0 $.
The \emph{closeness} between a snapshot $\vekt{s} (\vekt{\mu})$ and its rank $L$ approximation $ \vekt{s}_L (\vekt{\mu}) $ in the basis $\matr{V}$ can be quantified as their Euclidian distance using Equations (\ref{eq:approx}) and (\ref{eq:rbcoeffs}):

\begin{equation}
\label{eq:delta}
\delta_\text{POD}(\vekt{\mu}; \matr{V}) =
|| \vekt{s}(\vekt{\mu}) - \vekt{s}_L(\vekt{\mu}) ||_2 =
|| \vekt{s}(\vekt{\mu}) - \matr{V}\matr{V}^\top \vekt{s}(\vekt{\mu}) ||_2 \ .
\end{equation}

The \gls{SEY} \cite{Schmidt1907, Eckart1936} theorem states that the first $L$ left singular vectors of $\matr{S}$, i.e., $[\vekt{w}_1 | \cdots | \vekt{w}_{L}]$ are the optimal choice among all orthonormal bases $\hat{\matr{V}} \in \mathbb{R}^{N_h \times L}$

\begin{equation}
\label{eq:delta_train}
 \matr{V} = [\vekt{w}_1 | \cdots | \vekt{w}_L] = \operatorname*{argmin}_{ \hat{\matr{V}} } \delta_\text{POD}(\mathbb{P}_{tr};\hat{\matr{V}})
\end{equation}

with respect to minimizing the root-mean-square of the projection error over all training snapshots
\begin{equation}
\label{eq:Frobenius}
\delta_\text{POD}(\mathbb{P}_{tr}; \hat{\matr{V}})
= \sqrt{\frac{1}{N_{tr}} \sum_{\vekt{\mu} \in \mathbb{P}_{tr}} \delta^2_\text{POD}( \vekt{\mu}; \hat{\matr{V}} ) }
= \frac{1}{\sqrt{N_{tr}}}
|| \matr{S} - \hat{\matr{V}}\hat{\matr{V}}^\top \matr{S} ||_F \ ,
\end{equation}

where $||\cdot||_F$ denotes the Frobenius matrix norm.
The notation $\delta_\text{POD}(\mathbb{P}_{tr})$ is a generalization of $\delta_\text{POD}(\{\vekt{\mu}\})$ = $\delta_\text{POD}(\vekt{\mu})$, where we simply drop the brackets of the singleton set for ease of notation.
The \gls{SEY} theorem also states that the minimal associated projection error can be expressed as the root-squared sum of the left out singular values
\begin{equation}
\label{eq:delta_sumsigma}
  \delta_\text{POD}(\mathbb{P}_{tr} ;\matr{V}) = \sqrt{\sum_{l=L+1}^{N_{tr}} \sigma_l^2} \ .
\end{equation}

Computing the \gls{SVD} directly is prohibitively expensive, so the more efficient \emph{method of snapshots} is used in practice and in this work \cite{Sirovich}.
The idea is to first compute the eigenvalue decomposition of either $\matr{S}\matr{S}^\top /N_{tr} \in \mathbb{R}^{N_h \times N_h}$ or $ \matr{S}^\top\matr{S} /N_{tr} \in \mathbb{R}^{N_{tr} \times N_{tr}} $, depending on which one is smaller. This limits the computational complexity of \gls{POD} to be at worse cubic in the minimum of these dimensions $\mathcal{O}(\min\left\{N_{tr}, N_h\right\}^3)$ \cite{Pan1999}.
 Typically, the number of nodes $N_h$ is much larger than number of snapshots $N_{tr}$, so the eigendecomposition works out as $\matr{S}^\top\matr{S}/N_{tr} = \matr{Z}\matr{\Lambda}\matr{Z}^\top$, with $\matr{Z} \in \mathbb{R}^{N_{tr} \times N_{tr}}$ collecting the orthonormal eigenvectors and $\matr{\Lambda} = diag(\lambda_1, \cdots, \lambda_N)$ containing the real and positive eigenvalues of $\matr{S}^\top\matr{S}/N_{tr}$. With the above definition of \gls{SVD}, it follows that $\matr{\Sigma} = \matr{\Lambda}^{1/2}$ and $\matr{W} = \matr{S}\matr{Z}\matr{\Sigma}^{-1} / \sqrt{N_{tr}}$.

\begin{comment}
The resulting best rank $L$ approximation of the snapshot matrix is expressed in the \emph{thin} \gls{SVD} form \cite{Yu2019} as
$ \matr{S}_L = \matr{V}\matr{D}\matr{Y} $
with
$\matr{D}=diag(\sigma_1, \cdots, \sigma_L)$
and
$\matr{Y} = [\vekt{y}_1 | \cdots | \vekt{y}_L] \in \mathbb{R}^{N \times L}$ with $\vekt{y}_l =\vekt{z}_l $
for
$1 \le l \le L$.
\end{comment}

\subsection{Standardization}
\label{sec:preprocessing}
In the following, we describe the proposed \emph{standardization} steps to the non-intrusive \gls{ROM}. We start by providing theoretical bounds and motivation for the effectiveness of snapshot centering. Next, we show how standardization by singular values -- considered one of 'Tricks of the Trade' in machine learning \cite{LeCun} -- is a natural extension to the \gls{POD}-regression framework, which makes full use of the \gls{POD} and relieves the regression process.

\subsubsection{Snapshot centering}
\label{sec:centering}
Within the \gls{ROM} community, there is no clear consensus on whether to perform \gls{POD} on \emph{non-centered} data as described in \refsec{sec:POD} or on \emph{centered} data. The \emph{centered snapshot matrix}
$\matr{S}^c = [ \vekt{s}^c (\vekt{\mu}^{(1)}) | \cdots | \vekt{s}^c (\vekt{\mu}^{(N_{tr})}) ]$
is obtained by subtracting the mean over all columns $\bar{\vekt{s}} = 1/N_{tr} \sum_{\vekt{\mu} \in \mathbb{P}_{tr}} \vekt{s}(\vekt{\mu})$ from each snapshot:
\begin{equation}
\label{eq:q_h^c}
 \vekt{s}^c (\vekt{\mu}) = \vekt{s} (\vekt{\mu}) - \bar{\vekt{s}} \ .
\end{equation}

Some authors suggest that centering before \gls{POD} is 'customary' in \gls{ROM} \cite{Xiao2015, Ahmed2020, DMDbook2016}, but many counterexamples to this can also be found \cite{Guo2018, Hesthaven2018, Wang2019, Yu2019}.
In \cite{Ahmed2020} it is even argued for and in \cite{Chen2012} against centering based on empirical evidence from specific applications.

We prefer centering based on the following motivation.

According to the \gls{SEY} theorem, the optimal reduced basis for centered data $\matr{V}^c$ still consists of the left-hand singular vectors of $\matr{S}^c$ and the committed projection error can still be expressed as the sum of the truncated modes:

\begin{equation}
\label{eqn:delta_POD_centered}
\delta_\text{POD}(\mathbb{P}_{tr}; {\matr{V}^c}) =
 \sqrt{
  1/N_{tr}
  \sum_{\vekt{\mu} \in \mathbb{P}_{tr}} ||
  \vekt{s}^c(\vekt{\mu}) -
  \matr{V}^c {\matr{V}^c}^\top
  \vekt{s}^c(\vekt{\mu})
  ||_2^2
  }
  =
  \sqrt{\sum_{l=L+1}^{N_{tr}} (\sigma^c_l)^2}
   \ .
\end{equation}

After centering, ${\matr{S}^c}^\top\matr{S}^c/N_{tr}$ becomes the \emph{covariance} matrix and the singular values $\sigma^c$ of $\matr{S}^c/\sqrt{N_{tr}}$ can readily be interpreted as the population standard deviations along the respective reduced basis vectors.

Honeine \cite{Honeine2014} discusses the effects of applying \gls{SVD} to centered data. Importantly, Theorem 3 in \cite{Honeine2014} establishes upper and lower bounds for singular values $ \sigma_{l+1} \le \sigma^c_{l} \le \sigma_{l}$ for all $1 \le l \le N_{tr}$ from which bounds for the projection errors follow:
\begin{equation}
\delta_\text{POD}(\mathbb{P}_{tr}; \matr{V}_{L+1}) \le
\delta_\text{POD}(\mathbb{P}_{tr}; \matr{V}^c_L) \le
\delta_\text{POD}(\mathbb{P}_{tr}; \matr{V}_L).
\end{equation}

This means that with respect to the projection error, centering can effectively save one basis function, but also no more than that. According to Theorem 1 in \cite{Honeine2014}, the difference is diminishing for large $L$, but significant for small $L$:
\begin{equation}
 \delta^2_\text{POD}(\mathbb{P}_{tr}; \matr{V}_{L}) -
\delta^2_\text{POD}(\mathbb{P}_{tr}; \matr{V}^c_{L}) =
\begin{cases}
N_{tr}||\bar{\vekt{s}}||_2  & \text{for } L=0 \\
0                           & \text{for } L=N_{tr} \ .
\end{cases}
\end{equation}

Furthermore, centering also guarantees that the \emph{centered reduced coefficients} $\vekt{y}^c(\vekt{\mu})=\matr{V}^c \vekt{s}^c(\vekt{\mu})$ are also zero mean: $ 1/N_{tr} \sum_{\vekt{\mu} \in \mathbb{P}_{tr}} \vekt{y}^c(\vekt{\mu})  = \vekt{0} $, which is a common assumption in \gls{GPR} (see \refsec{sec:GPR}).

\subsubsection{Coefficient standardization}
\label{sec:scaling}

Another approach adopted in this work is coefficient standardization. The \emph{standardized reduced coefficients} are obtained by scaling the coefficients by the corresponding standard deviations:
\begin{equation}
    \vekt{y}^{s}(\vekt{\mu}) = (\matr{\Sigma}^c)^{-1} \vekt{y}^{c} (\vekt{\mu}) \ .
\end{equation}

$\vekt{y}^{s}(\vekt{\mu})$ have a zero mean and unit variance for each of its $L$ components and can thus be viewed as a multi-variate Z-score of the centered reduced coefficients.

This achieves dataset standardization across different problems and allows to reuse similar regression model architectures and learning processes, which are controlled by their \emph{hyperparameters}. Normally, these must be found using tedious trial-and-error or an extensive and computationally expensive hyperparameter-tuning. Standardization allows to set tight bounds on the search-space or even reuse hyperparameters (see \refsec{sec:Regression models}).
\reffig{fig:preprocessing} helps illustrate the described transformation steps.

\begin{figure}[tbp]
    \centering
    \input{figures/preprocessing}

    \caption[Preprocessing steps.]{
    The figure illustrates the effect of the data preprocessing steps described in Section 2.2.
    For some parameter instance $\vekt{\mu}$, $\vekt{s}(\vekt{\mu})$ denotes the two-dimensional ($N_h=2$) snapshot vector.
    Each element in the dataset is illustrated with a red dot, where $\vekt{e}_i$ is the $i$-th unit vector.
    Data is projected on all available bases ($L=N_h=2$).
    The final standardized data has zero mean and unit standard deviation along each basis vector, i.e., principal direction. The shift and scale invariance serves to improve interpretability and model hyperparameter reusability not only for different principal directions, but also across problems.}
    \label{fig:preprocessing}
\end{figure}

\subsubsection{Parameter standardization}
\label{sec:feature_standardization}

Similar to the standardized reduced coefficients, we also introduce the \emph{standardized parameters} as
\begin{equation}
\vekt{\mu}^s = (\vekt{\mu} - \bar{\vekt{\mu}})/\bar{\bar{\vekt{\mu}}}\end{equation}
with the mean
$\bar{\vekt{\mu}}
=
1/N_{tr} \sum_{\vekt{\mu} \in \mathbb{P}_{tr}} \vekt{\mu}
$
and the standard deviation
$\bar{\bar{\vekt{\mu}}}
=
\sqrt{1/N_{tr} \sum_{\vekt{\mu} \in \mathbb{P}_{tr}} (\vekt{\mu} - \bar{\vekt{\mu}})^2}
$. These again have zero mean and unit variance on any dataset following the same distribution as the training datset.
This standardization step is beneficial for regression tasks in problems, where different parameters are on different scales (see e.g. \refsec{sec:extruder}). Both \gls{RBF} and \gls{GPR} rely on distances in parameter space and would otherwise neglect the smaller parameters (see \refsec{sec:RBF} and \refsec{sec:GPR}). Similarly, \glspl{ANN} show faster convergence with standardized inputs \cite{LeCun98}.

\subsection{Errors}
\label{sec:errors}

Within the \gls{POD}-regression framework, errors are used to describe the quality of the reduced basis and the regression map. In \refsec{sec:errors_absolute}, we show how a careful definition of these errors and their aggregation across multiple samples lead to an intuitive and efficiently computable total error.
In \refsec{sec:errors_standardized}, we propose to standardize these errors, again connecting the results to familiar statistical concepts, achieving better interpretability and comparability across datasets, as well as scale and translation invariances.

\subsubsection{Absolute errors}
\label{sec:errors_absolute}

The \emph{absolute projection error} $\delta_\text{POD}$ has been already introduced in
\refeqn{eq:delta} as the Euclidian distance between the snapshot and its projection onto the reduced basis. The same holds with centering:
\begin{equation}
\delta_{\text{POD}}(\vekt{\mu}; \matr{V}^c) =
|| \vekt{s}(\vekt{\mu}) - \vekt{s}_L(\vekt{\mu})||_2
=
|| \vekt{s}^c(\vekt{\mu}) - \matr{V}^c{\matr{V}^c}^\top \vekt{s}^c(\vekt{\mu}) ||_2 \ .
\end{equation}

Similarly, the \emph{absolute regression error} $\delta_\text{REG}$ between the prediction
$\tilde{\vekt{s}}_L(\vekt{\mu})$
and the projection
$\vekt{s}_L(\vekt{\mu})$ is introduced
to describe how close the learned regression map $\tilde{\vekt{\pi}}^s$ approximates the observed regression map $\vekt{\pi}^s : \vekt{\mu^s} \mapsto \vekt{y}^s$:

\begin{equation}
\label{eq:delta_REG}
\begin{split}
\delta_{\text{REG}}(\vekt{\mu}; \matr{V}^c, \tilde{\vekt{\pi}}^s)
={}&
|| \tilde{\vekt{s}}_L(\vekt{\mu}) - \vekt{s}_L(\vekt{\mu}) ||_2
=
|| \tilde{\vekt{s}}^c_L(\vekt{\mu}) - \vekt{s}_L^c(\vekt{\mu}) ||_2
=
\\
={}&
|| \matr{V}^c\tilde{\vekt{y}}^c(\vekt{\mu}) - \matr{V}^c \vekt{y}^c(\vekt{\mu}) ||_2
=
|| \tilde{\vekt{y}}^c(\vekt{\mu}) - \vekt{y}^c(\vekt{\mu}) ||_2
=
\\
={}&
|| \matr{\Sigma}^c \left(\tilde{\vekt{y}}^s(\vekt{\mu}) - \vekt{y}^s(\vekt{\mu})\right) ||_2
=
|| \matr{\Sigma}^c \left(\tilde{\vekt{\pi}}^s(\vekt{\mu}^s) - \vekt{\pi}^s(\vekt{\mu}^s)\right) ||_2

 \ .
\end{split}
\end{equation}

Lastly, the \emph{absolute total error} $\delta_\text{POD-REG}$ is introduced as the distance between the true and predicted solutions to quantify the performance of the whole non-intrusive \gls{ROM}:

\begin{equation}
\label{eq:delta_PODREG}
\delta_{\text{POD-REG}}(\vekt{\mu}; \matr{V}^c, \tilde{\vekt{\pi}}^s)
=
|| \vekt{s}(\vekt{\mu}) - \tilde{\vekt{s}}_L(\vekt{\mu}) ||_2
=
|| \vekt{s}^c(\vekt{\mu}) - \tilde{\vekt{s}}^c_L(\vekt{\mu}) ||_2 \ .
\end{equation}

\reffig{fig:errors_illustrated} illustrates these three types of errors. Additionally, it can be seen that the total error is composed of two orthogonal components, namely the regression and projection errors:

\begin{equation}
\label{eq:delta_pythagorean}
{\delta_{\text{POD-REG}}(\vekt{\mu})}^2 =
{\delta_{\text{POD}}(\vekt{\mu})}^2 +
{\delta_{\text{REG}}(\vekt{\mu})}^2 \ .
\end{equation}

This can also be shown formally using the definition of the Euclidian vector norm $||\vekt{a}||_2^2=\vekt{a}^\top\vekt{a}$ and the fact that the predictions already belong to the linear span of the reduced basis $\matr{V}^c {\matr{V}^c}^\top \tilde{\vekt{s}}^c_L = \tilde{\vekt{s}}^c_L$:

\begin{equation}
\begin{split}
|| {\vekt{s}^c} - {\matr{V}^c}{\matr{V}^c}^\top{\vekt{s}^c} ||_2^2 +& || {\matr{V}^c}{\matr{V}^c}^\top{\vekt{s}^c} - {\tilde{\vekt{s}^c}_L} ||_2^2 =

\\ =
  {{\vekt{s}^c}^\top}{\vekt{s}^c}
- {{\vekt{s}^c}^\top}{\matr{V}^c}&{{\matr{V}^c}^\top}{\tilde{\vekt{s}}^c_L}
+ {{\tilde{\vekt{s}}^c_L}{}^\top} {\tilde{\vekt{s}}^c_L} =
\\ =
  {{\vekt{s}^c}^\top}{\vekt{s}^c}
- {{\vekt{s}^c}^\top}&{\tilde{\vekt{s}}^c_L}
+ {{\tilde{\vekt{s}}^c_L}{}^\top} {\tilde{\vekt{s}}^c_L} =
\\ =
|| \vekt{s}^c -& {\tilde{\vekt{s}}^c_L} ||_2^2

\end{split}
\end{equation}

\begin{figure}[tbp]
\centering
\input{figures/errors_illustrated}

\caption[Relation between errors.]{The total error $\delta_\text{POD-REG}$ can be decomposed into two orthogonal components -- the projection error $\delta_\text{POD}$ perpendicular to the span of centered basis $V^c$ and the regression error $\delta_\text{REG}$ belonging to the span. Illustration uses three degrees of freedom $N_h=3$ and two reduced basis vectors $L=2$.}
\label{fig:errors_illustrated}
\end{figure}

Lastly, we note that the equivalence in \refeqn{eq:delta_REG} paired with with \refeqn{eq:delta_pythagorean} allows to compute both $\delta_{\text{REG}}$ and $\delta_{\text{POD-REG}}$ in $\mathcal{O}(L)$ operations given that $\delta_{\text{POD}}$ is regression independent and can be pre-computed once.
In contrast to the naive approach via computing $\tilde{\vekt{s}}^c_L=\vekt{V}^c\tilde{\vekt{s}}^c$ requiring $\mathcal{O}(N_h L)$ operations due to matrix-vector multiplication, the presented approach allows to efficiently use and monitor all errors during the training of the regression models.

The corresponding \emph{aggregate} projection, regression and total errors $\delta\left(\mathbb{P}\right)_*, *\in\left\{\text{POD},\text{REG},\text{POD-REG}\right\}$ on any dataset $\mathbb{P}$ of cardinality $N$ are defined via the \gls{RMSE}:
\begin{equation}
    \delta_*\left( \mathbb{P} \right)
    :=
    \sqrt{
    1/N \textstyle\sum_{\vekt{\mu} \in \mathbb{P}}
    {\delta_*\left( \vekt{\mu} \right)}^2
    } \ .
\end{equation}
This choice is motivated by the form of the \gls{SEY} theorem and treats both physical and parameter dimensions equally.

Hence, \refeqn{eq:delta_pythagorean} holds for both individual and aggregate errors.

\subsubsection{Standardized errors}
\label{sec:errors_standardized}

The use of absolute errors makes it difficult to assess and interpret the performance of \gls{ROM}, especially across different datasets.

Thus we introduce the corresponding \emph{standardized projection, regression} and \emph{total errors} in the following way:
\begin{equation}
    \varepsilon_*\left(\mathbb{P}\right) := \frac
    { \delta_*\left(\mathbb{P}\right) }
    {\sqrt{
        1/N_{tr} \sum_{\vekt{\mu} \in \mathbb{P}_{tr}}
        || \vekt{s}\left(\vekt{\mu}\right) - \bar{\vekt{s}} ||^2
    }} \ .
\end{equation}

The choice of the denominator is motivated by its relation to the standard deviation of the training set. Hence, $\varepsilon$ can be interpreted as a multivariate extension of \emph{Z-score} and describes the ratio between the Euclidian distance of two vectors and the standard deviation of the particular dataset.
Always naively predicting the mean vector $ \bar{\vekt{s}} $ leads to 100\% aggregated total error $\varepsilon_{\text{POD-REG}}(\mathbb{P}_{tr})=1$.
Furthermore, the denominator can be related to \gls{POD}, since
$ 1/N_{tr} \sum_{\vekt{\mu} \in \mathbb{P}_{tr}}
|| \vekt{s}\left(\vekt{\mu}\right) - \bar{\vekt{s}} ||_2^2
= 1/N_{tr} \Tr({\matr{S}^c}^\top \matr{S}^c)
= \sum_{n=1}^{N_{tr}}{(\sigma^c_n)^2}
$ .
Consequentially, it is easy to relate the standardized projection error to \gls{SEY} and \refeqn{eq:delta_sumsigma}:
\begin{equation}
{\varepsilon^2_\text{POD}(\mathbb{P}_{tr}; \matr{V}^c_L)} =
\frac
{\sum_{n=L+1}^{N_{tr}}{(\sigma^c_n)^2}}
{\sum_{n=1}^{N_{tr}}{(\sigma^c_n)^2}}  \ .
\end{equation}
This corresponds to a commonly used measure
known as the \emph{relative information loss}, which can be used to guide the choice of number of bases $L$.

Furthermore, $\varepsilon$ is invariant to scaling and translation of snapshots and thus independent of the choice of units. This alleviates the need for designing dimensionless or normalized problems in order to interpret the errors and altogether allows for a more data-driven approach to non-intrusive \gls{ROM}.
Finally, within the broader classification of existing error measures in machine learning, our definition of $\varepsilon$ can be viewed as a multivariate extension of the root-relative-squared error, option 2 as classified by Botchkarev \cite{Botchkarev2018}.

\subsection{Time-resolved problems}
\label{sec:unsteadyRB}

In general, authors within \gls{ROM} community approach time-resolved problems by treating time as a continuous parameter \cite{Yu2019}. Wang et al. \cite{Wang2019} point out two issues that arise using this approach. Firstly, for problems with many time steps $N_t$, the snapshot matrix becomes wide $\matr{S} \in \mathbb{R} ^ {N_h {\times} N_{tr} N_t}$ and the cost of the \gls{POD} $\mathcal{O}\left( \min \left\{ N_h , N_{tr} N_t \right\}^3 \right)$ becomes high even using the method of snapshots because the smallest dimension is increased significantly.
Secondly, the regression model must be able to predict the solution at arbitrary time and parameter values. Thus it needs to learn both time and parameter dynamics.
Several solutions to these problems have been proposed, mainly using two-level \gls{POD} \cite{Audouze2013, Chen2018, Xiao2015, Guo2019, Wang2019} or operator inference \cite{Peherstorfer2016}.

We propose an alternative method by relaxing the second requirement of needing to predict solutions at any arbitrary time instance. Instead, we treat time in discrete levels, similar to physical space.
This is inspired by typical numerical schemes, which also provide the solution at discrete spatial and temporal points, with any intermediate values being interpolated.
This approach ensures that the snapshot matrix becomes taller not wider $\matr{S} \in \mathbb{R} ^ {N_h N_t {\times} N_{tr} }$, preserving the smaller dimension $N_{tr}$ and ensuring efficiency of the \gls{POD}.
Each resulting basis vector spans space and time, and as a consequence the regression model must only account for parameter dynamics. This enables much simpler models and significantly reduces the effort of the fitting process.


TITLE: Is this an unknown pattern in prime numbers?
QUESTION [4 upvotes]: I am trying to figure out if the pattern I've found concerning twin primes is a known pattern or not. It turns out that with every set of twin primes, if the higher of the two numbers is converted to radix 7, and then the individual digits of the number are added together and continually added together until a 1 or 2 digit number is leftover, the number is always equal to 6 mod +1.
The lower of the two numbers is always 6 mod -1 with same calculation.
Examples:
Lower twin (radix 10)/ Lower twin (radix 7)/ [sum of digits to 2 digits]/ MOD 6

59 / 113 / 5 / 5

71 / 131 / 5 / 5

101 / 203 / 5 / 5

107 / 212 / 5 / 5

137 / 254 / 11 / 5

149 / 302 / 5 / 5

179 / 344 / 11 / 5

191 / 362 / 11 / 5

197 / 401 / 5 / 5

227 / 443 / 11 / 5

239 / 461 / 11 / 5

269 / 533 / 11 / 5

281 / 551 / 11 / 5

311 / 623 / 11 / 5

347 / 1004 / 5 / 5

419 / 1136 / 11 / 5

431 / 1154 / 11 / 5

461 / 1226 / 11 / 5

521 / 1343 / 11 / 5

569 / 1442 / 11 / 5

599 / 1514 / 11 / 5

617 / 1541 / 11 / 5

641 / 1604 / 11 / 5

659 / 1631 / 11 / 5

809 / 2234 / 11 / 5

821 / 2252 / 11 / 5

827 / 2261 / 11 / 5

857 / 2333 / 11 / 5

881 / 2366 / 17 / 5

1019 / 2654 / 17 / 5

1031 / 3002 / 5 / 5

1049 / 3026 / 11 / 5

1061 / 3044 / 11 / 5

1091 / 3116 / 11 / 5

1151 / 3233 / 11 / 5

1229 / 3404 / 11 / 5

1277 / 3503 / 11 / 5

1289 / 3521 / 11 / 5

1301 / 3536 / 17 / 5

1319 / 3563 / 17 / 5

1427 / 4106 / 11 / 5

1451 / 4142 / 11 / 5

1481 / 4214 / 11 / 5

1487 / 4223 / 11 / 5

1607 / 4454 / 17 / 5

1619 / 4502 / 11 / 5

963426767 / 32605664252 / 41 / 5

963427259 / 32605665554 / 47 / 5

963427301 / 32605665644 / 47 / 5

963427559 / 32605666463 / 47 / 5

963427919 / 32606000516 / 29 / 5

963428021 / 32606001023 / 23 / 5

963428099 / 32606001164 / 29 / 5

963428561 / 32606002424 / 29 / 5

963428861 / 32606003333 / 29 / 5

963428957 / 32606003531 / 29 / 5

963429167 / 32606004251 / 29 / 5

963430019 / 32606006606 / 35 / 5

963430079 / 32606010023 / 23 / 5

963430289 / 32606010443 / 29 / 5

963431177 / 32606013152 / 29 / 5

963431321 / 32606013446 / 35 / 5

963431477 / 32606014061 / 29 / 5

963431717 / 32606014553 / 35 / 5

963432131 / 32606016014 / 29 / 5

963432917 / 32606021216 / 29 / 5

963432989 / 32606021351 / 29 / 5

963433319 / 32606022332 / 29 / 5

963433439 / 32606022563 / 35 / 5

963433697 / 32606023412 / 29 / 5

963434411 / 32606025452 / 35 / 5

963434579 / 32606026112 / 29 / 5

963434609 / 32606026154 / 35 / 5

963434891 / 32606030036 / 29 / 5

963435227 / 32606031026 / 29 / 5

963435491 / 32606031554 / 35 / 5

963436037 / 32606033264 / 35 / 5

963436601 / 32606035031 / 29 / 5

963437261 / 32606036663 / 41 / 5

963437399 / 32606040251 / 29 / 5

963437927 / 32606041634 / 35 / 5

963437939 / 32606041652 / 35 / 5

963438017 / 32606042123 / 29 / 5

963438041 / 32606042156 / 35 / 5

Higher twin (radix 10)/ Higher twin (radix 7)/ [sum of digits to 2 digits]/ MOD 6

571 / 1444 / 13 / 1

601 / 1516 / 13 / 1

619 / 1543 / 13 / 1

643 / 1606 / 13 / 1

661 / 1633 / 13 / 1

811 / 2236 / 13 / 1

823 / 2254 / 13 / 1

829 / 2263 / 13 / 1

859 / 2335 / 13 / 1

883 / 2401 / 7 / 1

1021 / 2656 / 19 / 1

1033 / 3004 / 7 / 1

1051 / 3031 / 7 / 1

1063 / 3046 / 13 / 1

1093 / 3121 / 7 / 1

1153 / 3235 / 13 / 1

1231 / 3406 / 13 / 1

1279 / 3505 / 13 / 1

1291 / 3523 / 13 / 1

1303 / 3541 / 13 / 1

1321 / 3565 / 19 / 1

1429 / 4111 / 7 / 1

1453 / 4144 / 13 / 1

1483 / 4216 / 13 / 1

961750903 / 32555514331 / 37 / 1

961751209 / 32555515246 / 43 / 1

961752301 / 32555521366 / 43 / 1

961752349 / 32555521465 / 43 / 1

961752553 / 32555522206 / 37 / 1

961753789 / 32555525623 / 43 / 1

961753831 / 32555526013 / 37 / 1

961754011 / 32555526361 / 43 / 1

961754071 / 32555526505 / 43 / 1

961754461 / 32555530603 / 37 / 1

961755019 / 32555532331 / 37 / 1

961757059 / 32555541304 / 37 / 1

961757311 / 32555542114 / 37 / 1

961757431 / 32555542345 / 43 / 1

961757683 / 32555543155 / 43 / 1

961758673 / 32555546101 / 37 / 1

961759111 / 32555550265 / 43 / 1

961759483 / 32555551336 / 43 / 1

961759831 / 32555552344 / 43 / 1

961759861 / 32555552416 / 43 / 1

961760119 / 32555553235 / 43 / 1

961760719 / 32555555053 / 43 / 1

961761013 / 32555555653 / 49 / 1

961761139 / 32555556223 / 43 / 1

961761343 / 32555556634 / 49 / 1

961761403 / 32555560051 / 37 / 1

961761571 / 32555560411 / 37 / 1

961762033 / 32555561641 / 43 / 1

961762591 / 32555563366 / 49 / 1

I have other questions related to prime numbers but first want to see how valid or known this part is before I continue. I am not a mathematician.

REPLY [5 votes]: \begin{align}
d_0 + 7d_1 + 7^2 d_2 + 7^3 d_3 + \cdots & \equiv d_0 + 1d_1 + 1^2 d_2 + 1^3 d_3+\cdots & &\mod 6 \\[10pt]
& \equiv d_0 + d_1 + d_2 + d_3 + \cdots & & \mod 6
\end{align}
What is at work here is something that says if $a\equiv A\bmod 6$ and $b\equiv B\bmod 6$ then $ab\equiv AB\bmod6$. Proving that takes a bit of elementary algebra. Applying it here we have $7\equiv 1;$ therefore $7\times7\equiv 1\times 1,$ etc.
The fact that twin primes are always of the form $6n\pm1,$ plus the facts above lead to the conclusion that the pattern you've identified will persist.


TITLE: Good book on integral calculus (improper integrals, integrals with parameters, special functions)
QUESTION [10 upvotes]: Can you recommend a good book (with theoretical results with proofs, and with plenty of solved problems and examples) on the topics of improper integrals, (improper) integrals with parameters, special functions (Beta, Gamma, ...)?

REPLY [2 votes]: You may have a look at the affordable and accessible Schaum's Outline of Advanced Calculus. Not too much advanced :)