This document provides an explanation for the numerical methods provided by the class \code{Rational Interpolant}\footpath{./rrompy/reduction_methods/standard/rational_interpolant.py} and daughters, e.g. \code{Rational Interpolant Greedy}\footpath{./rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py}, as well as most of the pivoted approximants\footpath{./rrompy/reduction_methods/pivoted/{,greedy/}rational_interpolant_*.py}.
Most of the focus will be dedicated to the impact of the (\code{rationalMode},\code{functionalSolve}) parameters, whose allowed values are
\begin{itemize}
\item (\code{MINIMAL},\code{NORM}) (default): see Section~\ref{sec:min:norm}; allows for repeated sample points.
\item (\code{MINIMAL},\code{DOMINANT}): see Section~\ref{sec:min:dominant}; allows for repeated sample points.
\item (\code{BARYCENTRIC},\code{NORM}): see Section~\ref{sec:min:barycentricnorm}; does not allow for a Least Squares approach; undefined for more than one parameter.
\item (\code{BARYCENTRIC},\code{DOMINANT}): see Section~\ref{sec:min:barycentricaverage}; does not allow for a Least Squares approach; undefined for more than one parameter.
\item (\code{STANDARD},\code{NORM}): see Section~\ref{sec:out:norm}; allows for repeated sample points.
\item (\code{STANDARD},\code{DOMINANT}): see Section~\ref{sec:out:dominant}; allows for repeated sample points.
\end{itemize}
We restrict the discussion to the single-parameter case. The main reference throughout the present document is \cite{Pradovera}.
\section{Aim of approximation}
We seek an approximation of $u:\C\to V$, with $(V,\inner{\cdot}{\cdot})$ a complex\footnote{The inner product is linear (resp. conjugate linear) in the first (resp. second) argument: $\inner{\alpha v}{\beta w}=\alpha\overline{\beta}\inner{v}{w}$.} Hilbert space (with induced norm $\norm{\cdot}$), of the form $\widehat{p}/\widehat{q}$, where $\widehat{p}:\C\to V$ and $\widehat{q}:\C\to\C$. The target $u$ might be high-dimensional (for instance, the solution of a PDE after Finite Element discretization) or low(er)-dimensional (for instance, a functional of the above-mentioned PDE solution). For a given denominator $\widehat{q}$, the numerator $\widehat{p}$ is found by interpolation of $\widehat{q}u$. Hence, here we focus on the computation of the denominator $\widehat{q}$.
Other than the choice of target function $u$, the parameters which affect the computation of $\widehat{q}$ are:
\begin{itemize}
\item$\code{mus}\subset\C$ ($\{\mu_j\}_{j=1}^S$ below); for \code{BARYCENTRIC}, the $S$ points must be distinct.
\item$\code{N}\in\N$ ($N$ below); for \code{BARYCENTRIC}, $N$ must equal $S-1$.
To simplify the notation, we set $E=S-1$. For simplicity, we will consider only the case of $S$ distinct sample points. One can deal with the case of confluent sample points by extending the standard (Lagrange) interpolation steps to Hermite-Lagrange ones.
In the following, we will make use of the polynomial interpolation operator $\I^M:(\C\times W)^S\to\mathbb{P}^M(\C;W)$, where $M$ is an integer (either $N$ or $E$), and $W$ a Banach space (either $\C$ or $V$). We define its action on samples $\left((\mu_j,\psi_j)\right)_{j=1}^S\in(\C\times W)^S$ as
In practice, the polynomial basis $\{\phi_i\}_{i=0}^M$ is determined by the value of $\code{polybasis}$:
\begin{itemize}
\item If $\code{polybasis}=\code{"CHEBYSHEV"}$, then $\phi_k(\mu)=\mu^k$ for $k\in\{0,1\}$ and $\phi_k(\mu)=2\mu\phi_{k-1}(\mu)-\phi_{k-2}(\mu)$ for $k\geq2$.
\item If $\code{polybasis}=\code{"LEGENDRE"}$, then $\phi_k(\mu)=\mu^k$ for $k\in\{0,1\}$ and $\phi_k(\mu)=(2-1/k)\mu\phi_{k-1}(\mu)-(1-1/k)\phi_{k-2}(\mu)$ for $k\geq2$.
\item If $\code{polybasis}=\code{"MONOMIAL"}$, then $\phi_k(\mu)=\mu^k$ for $k\geq0$.
\end{itemize}
Moreover, it will prove useful to define the ``snapshot rank'' $r=r(\{u(\mu_j)\}_{j=1}^S)=\text{dim}~\text{span}\{u(\mu_j)\}_{j=1}^S$. This quantity is bounded from above by $S$ and by the dimension of $V$ (e.g., $r=1$ if $V=\C$).
where $(\star)$ is a normalization condition (which changes depending on \code{functionalSolve}) to exclude the trivial minimizer $\widehat{q}\equiv0$. The methods described differ in terms of the constraint $(\star)$, as well as of the degrees of freedom which are chosen to represent the denominator $q$.
For \eqref{eq:glob1} to be well-defined (unique optimizer up to unit scaling), in addition to the condition $N\leq S-1$, we must have (by balance of degrees of freedom vs. constraints) $N\leq r$
\subsection{Polynomial coefficients as degrees of freedom}\label{sec:mri:poly}
If the polynomial basis $\{\phi_i\}_{i=0}^E$ in \eqref{eq:normal} is hierarchical (as the three ones above), then the $E$-th derivative of $\I^E$ is proportional to the coefficient $c_E$, and we have
Using the Kronecker delta ($\delta_{ij}=1$ if $i=j$ and $\delta_{ij}=0$ if $i\neq j$), the last term $[q(\mu_j)u(\mu_j)]_{j=1}^S\in V^S$ can be factored into
where $\widetilde{\Phi}$ is the $S\times(N+1)$ matrix obtained by extracting the first $N+1$ columns of $\Phi$. We remark that we have expanded the polynomial $q$ using the basis\footnote{In theory, nothing prevents us from using different bases for $\I^E$ and $q$, cf. Section~\ref{sec:observations}.}$\{\phi_i\}_{i=0}^N$: $q(\mu)=\sum_{i=0}^Nq_i\phi_i(\mu)$, with coefficients $\mathbf{q}=[q_i]_{i=0}^N\in\C^{N+1}$.
Combining \eqref{eq:poly1} and \eqref{eq:poly2}, it is useful to consider the $(N+1)\times(N+1)$ Hermitian matrix with entries ($0\leq i,i'\leq N$)
If $(\star)$ is quadratic (resp. linear) in $\mathbf{q}$, then we can cast the computation of the denominator as a quadratically (resp. linearly) constrained quadratic program involving $G$.
We constrain $\widehat{\mathbf{q}}$ to have unit (Euclidean) norm. The resulting optimization problem can be cast as a minimal (normalized) eigenvector problem for $G$ in \eqref{eq:poly3}. More explicitly,
We constrain $\widehat{q}_N=1$, thus forcing $q$ to be monic, with degree exactly $N$. Given $G$ in \eqref{eq:poly3}, the resulting optimization problem can be solved rather easily as:
\subsection{Barycentric coefficients as degrees of freedom}\label{sec:bary}
Here we assume that the sample points are distinct, and that $N=E=S-1$, so that, in particular, $\Phi=\widetilde{\Phi}$. Considering the constraint $N\leq r$, we can deduce that the approach presented in this section can only be applied if the rank of the snapshots is either full or defective by 1.
We can choose for convenience a non-hierarchical rational basis, dependent on the sample points, for $q$ and $\I^E$, taking inspiration from barycentric interpolation:
Evaluation of such basis at the sample points leads to unbounded results. To correct this, we symbolically weigh each row by the nodal polynomial $\prod_{j=1}^S(\cdot-\mu_j)$, causing cancellation of the unbounded terms. The resulting ``effective'' polynomial basis is
Considering \eqref{eq:bary3}, it is useful to define the $S\times S$ Hermitian ``snapshot Gramian'' matrix with entries ($0\leq i,i'\leq N$)
\begin{equation}\label{eq:bary4}
G_{ii'}=\inner{u(\mu_{i'+1})}{u(\mu_{i+1})}.
\end{equation}
So, once again, if $(\star)$ is quadratic (resp. linear) in $\mathbf{q}$, then we can cast the computation of the denominator as a quadratically (resp. linearly) constrained quadratic program involving $G$.
Before specifying the kind of normalization enforced, it is important to make a remark on numerical stability. The basis in \eqref{eq:bary1} is actually just a ($i$-dependent) factor away from being the Lagrangian one (for which $\phi_i(\mu_j)$ would equal $\delta_{(i+1)j}$ instead of
as it does in our case). As such, it is generally a bad idea to numerically evaluate $q$ starting from its expansion coefficients with respect to $\{\phi_i\}_{i=0}^N$. We get around this by exploiting the following trick, whose foundation is in \cite[Section 2.3.3]{Klein}: the roots of $\widehat{q}=\sum_{i=0}^N\widehat{q}_i\phi_i$ are the $N$ finite eigenvalues $\lambda$ of the generalized $(N+2)\times(N+2)$ eigenproblem
\begin{equation}\label{eq:bary5}
\textrm{Det}\left(\begin{bmatrix}
0 &\widehat{q}_0 &\cdots&\widehat{q}_N\\
1 &\mu_1 &&\\
\vdots&&\ddots&\\
1 &&&\mu_S
\end{bmatrix}-\begin{bmatrix}
0 &&&\\
& 1 &&\\
&&\ddots&\\
&&& 1
\end{bmatrix}\lambda\right)=0.
\end{equation}
This computation is numerically more stable than most other manipulations of a polynomial in the basis \eqref{eq:bary1}. Once the roots of $\widehat{q}$ have been computed, one can either convert it to nodal form
or forgo using $\widehat{q}$ completely, in favor of a Heaviside-like approximation involving the newly computed roots $\{\widehat{\lambda}_i\}_{i=1}^N$ as poles:
See the final paragraph in Section~\ref{sec:observations} for a slightly more detailed motivation of why the Heaviside form of the approximant might be more useful than the standard rational one $\widehat{p}/\widehat{q}$ in practice.
We constrain $\widehat{\mathbf{q}}$ to have unit (Euclidean) norm. The resulting optimization problem can be cast as a minimal (normalized) eigenvector problem for $G$ in \eqref{eq:bary4}. More explicitly,
We constrain $\sum_{i=0}^N\widehat{q}_i=1$, so that the polynomial $\widehat{q}$ is monic. Given $G$ in \eqref{eq:bary4}, the resulting optimization problem can be solved rather easily as:
\subsection{Minor observations for MRI}\label{sec:observations}
\begin{itemize}
\item If $N=E$, normal equations are not necessary to compute $\I^E$, since $\Phi$ is square and can be inverted directly. However, in practical applications, it may be useful to decrease the degree $E$ of the interpolant (which, in our presentation, we kept fixed to $S-1$ for simplicity) to overcome numerical instabilities which may arise in the (pseudo-)inversion of $\Phi$. If this happens, $\Phi$ becomes non-square, and normal equations are the only option.
\item For \code{BARYCENTRYC}, a specific choice of polynomial basis for $\I^E$ was used to diagonalize the functional. Under the assumptions that the sample points are distinct (and that $E+1=S$), one can employ the quasi-Lagrangian basis \eqref{eq:bary1} to expand $\I^E$ in the other approaches as well, thus simplifying significantly the structure of \eqref{eq:glob1}:
This is independent of the basis used to expand $q$ and, numerically, has repercussions only on the computation of the term $\big(\Phi^H\Phi\big)^{-1}\Phi^H$ in \eqref{eq:poly3}.
\item In general, \code{NORM} can be expected to be more numerically stable than \code{DOMINANT}. This is due to the fact that the normalization is enforced in a more numerically robust fashion.
\item If the snapshots are orthonormalized via \code{POD}\footpath{./rrompy/sampling/engines/sampling_engine_pod.py}, a simple unitary transformation allows to replace $V$ with $\C^r$. As a consequence, all the $V$-inner products (resp. norms) can be recast as Euclidean inner products (resp. norms) involving the $R$ factor of the generalized ($V$-orthonormal) QR decomposition of the snapshots.
\item If a univariate rational surrogate is built in the scope of multivariate pole-matching-based pivoted approximation\footpath{./rrompy/reduction_methods/pivoted/*}, the rational approximant is converted into a Heaviside/nodal representation when different surrogates are combined. As such, the \code{BARYCENTRIC} approach may be preferable to avoid extra computations, as well as additional round-off artifacts.
where $(\star)$ is a normalization condition (which changes depending on \code{functionalSolve}) to exclude the trivial minimizer $\widehat{q}\equiv0$.
For \eqref{eq:glob2} to be well-defined (unique optimizer up to unit scaling), in addition to the conditions $N\leq S-1$ and $M\leq S-2$, we must have (by balance of degrees of freedom vs. constraints) $M+N/r+1\leq S$.
As before, we can expand the target functional by linear algebra
We observe that $I-\Phi\big(\Phi^H\Phi\big)^{-1}\Phi^H$ (with $I$ the $S\times S$ identity matrix) corresponds to the (orthogonal) projection onto the orthogonal complement of the range of $\Phi$, so that we may express it as $\Psi\Psi^H$, with the $S-M-1$ columns of $\Psi$ being an orthonormal basis for the orthogonal complement of the range of $\Phi$ (which can be found by completion of the QR decomposition of $\Phi$). Hence, we have
As above, we can interpret \eqref{eq:poly52} as the minimization of a quadratic form represented by a $(N+1)\times(N+1)$ Hermitian matrix with entries ($0\leq i,i'\leq N$)
The matrix $G$ above can also be expressed as $\widetilde{\Phi}^H\big(\mathring{G}\bullet\left(\Psi\Psi^H\right)\big)\widetilde{\Phi}$, with $\mathring{G}$ the snapshot Gramian \eqref{eq:bary4} and $\bullet$ the Hadamard product. If a Cholesky decomposition $\mathring{G}=\mathring{R}^H\mathring{R}$ is available, then $G=\widetilde{\Phi}^H(\mathring{R}\odot\Psi^H)^H(\mathring{R}\odot\Psi^H)\widetilde{\Phi}$, with $\odot$ the Khatri-Rao product.
If $(\star)$ is quadratic (resp. linear) in $\mathbf{q}$, then we can cast the computation of the denominator as a quadratically (resp. linearly) constrained quadratic program involving $G$.
We constrain $\widehat{\mathbf{q}}$ to have unit (Euclidean) norm. The resulting optimization problem can be cast as a minimal (normalized) eigenvector problem for $G$ in \eqref{eq:poly53}. More explicitly,
We constrain $\widehat{q}_N=1$, thus forcing $q$ to be monic, with degree exactly $N$. Given $G$ in \eqref{eq:poly53}, the resulting optimization problem can be solved rather easily as: