Interpolatory rational model order reduction of parametric problems lacking uniform inf-sup stability

TL;DR

This paper introduces a minimal rational interpolation method for efficiently approximating parametric PDE solution maps with meromorphic structure, demonstrating convergence and adaptivity through theoretical analysis and numerical experiments.

Contribution

The paper presents a novel minimal rational interpolation approach for model order reduction of parametric PDEs with meromorphic solutions, including convergence analysis and adaptive error estimation.

Findings

The method accurately captures resonance features of the solution map.

The approach converges under certain structural conditions.

Numerical experiments confirm theoretical predictions and effectiveness.

Abstract

We present a technique for the approximation of a class of Hilbert space-valued maps which arise within the framework of Model Order Reduction for parametric partial differential equations, whose solution map has a meromorphic structure. Our MOR stategy consists in constructing an explicit rational approximation based on few snapshots of the solution, in an interpolatory fashion. Under some restrictions on the structure of the original problem, we describe a priori convergence results for our technique, hereafter called minimal rational interpolation, which show its ability to identify the main features (e.g. resonance locations) of the target solution map. We also investigate some procedures to obtain a posteriori error indicators, which may be employed to adapt the degree and the sampling points of the minimal rational interpolant. Finally, some numerical experiments are carried out to confirm the theoretical results and the effectiveness of our technique.

Equations

$$\text{for }k\in K\subset \mathbb{C},\text{find }u(k)\in V\text{s.t.}-\Delta u(k)-k^{2}u(k)\stackrel{(W)}{=}f\in W\text{,}$$

$${\mathbb{P}}_N^{\star }\left(\mathbb{C}\right)=\left\{Q\in {\mathbb{P}}_N\left(\mathbb{C}\right),{\left|\left|\left|Q\right|\right|\right|}_N=1\right\}\text{,}$$

$${\mathbb{P}}_{S-1}\left(\mathbb{C};V\right)=\left\{\mu \mapsto \sum _{i=0}^{S-1}{p}_i{\mu }^i,p_i\in V\text{ for }i=0,\dots ,S-1\right\}$$

$$u_{N}^{\Xi_{S}}(\mu)=I^{\Xi_{S}}\left(uQ_{N}^{\Xi_{S}}\right)(\mu)\,{\Big{/}}\,Q_{N}^{\Xi_{S}}(\mu),\quad\mu\in\mathbb{C},Q_{N}^{\Xi_{S}}(\mu)\neq 0\text{,}$$

$${\mathbb{P}}_N\left(\mathbb{C}\right)\ni Q\mapsto j^{{\Xi }_S}(Q{)}^2={\Vert \frac{1}{(S-1)!}\frac{d^{S-1}}{d{\mu }^{S-1}}{I}^{{\Xi }_S}(uQ)\Vert }_V^2\text{,}$$

$$\mathbf{w}({\mu }_j)={\left(⟨u({\mu }_j),{\phi }_i{⟩}_V\right)}_{i=1}^S,\text{so that}u({\mu }_j)=\sum _{i=1}^S{w}_i({\mu }_j){\phi }_i,\text{for }j=1,\dots ,S\text{;}$$

$$\Psi ={\left(\frac{1}{(S-1)!}\frac{d^{S-1}}{d{\mu }^{S-1}}{I}^{{\Xi }_S}(w_i{\psi }_l)\right)}_{i=1,\dots ,S,l=0,\dots ,N}\in {\mathbb{C}}^{S\times (N+1)}\text{,}$$

$$j^{{\Xi }_S}(Q{)}^2=$$

$$=$$

$$\text{for }\mu \in K\subset \mathbb{C},\text{find }u(\mu )\in \mathrm{Dom}({\mathcal{F}}_{\mu })\text{s.t.}{\mathcal{F}}_{\mu }(u(\mu ))=0\text{,}$$

$$u:K^{\star }$$

$$\mu$$

$$u(\mu )=\sum _{\lambda \in \Lambda }\frac{v_{\lambda }}{\lambda -\mu },\mu \in \mathbb{C}\text{.}$$

$${\varphi }_K:\overline{\mathbb{C}\setminus K}\to \mathbb{C}\setminus \mathcal{B}(0,\mathrm{Cap}(K))=\{z\in \mathbb{C},\left|z\right|\ge \mathrm{Cap}(K)\}\text{,}$$

$$\sqrt{\frac{\ell (A)}{\pi }}\le \mathrm{Cap}(A)\le \frac{1}{2}\underset{z,w\in A}{\max }\left|z-w\right|\text{,}$$

$${\Phi }_K(\mu )=\{\begin{array}{ll}\left|{\varphi }_K(\mu )\right|& \text{if }\mu \notin K,\\ \mathrm{Cap}(K)& \text{if }\mu \in K\end{array}\text{.}$$

$${\Phi }_K({\lambda }_1)\le {\Phi }_K({\lambda }_2)\le \dots \text{.}$$

$$\underset{S\to \infty }{\mathrm{limsup}}{\left|{\omega }^{{\Xi }_S^K}\right|}^{1/S}\le {\Phi }_K\text{uniformly over compact subsets of }\mathbb{C}\text{,}$$

$$\underset{S\to \infty }{\lim }{\left|{\omega }^{{\Xi }_S^K}\right|}^{1/S}={\Phi }_K\text{uniformly over compact subsets of }\mathbb{C}\setminus \partial K\text{.}$$

$$c_{{\mu }_0}\prod _{j=1}^D\left|{\mu }_0-z_j\right|\le {\left|\left|\left|\prod _{j=1}^D\left(\cdot -z_j\right)\right|\right|\right|}_N\le C_{{\mu }_0}\prod _{j=1}^D\left(R_{{\mu }_0}+\left|{\mu }_0-z_j\right|\right)\text{,}$$

$${\left|\left|\left|Q\right|\right|\right|}_{{\mu }_0,N}=\sqrt{\sum _{j=0}^N{\left|q_j\right|}^2},\text{where}Q=\sum _{j=0}^N{q}_j{\left(\cdot -{\mu }_0\right)}^j\text{.}$$

$$c_{{\mu }_0,N}\prod _{j=1}^D\left|{\mu }_0-z_j\right|\le {\left|\left|\left|\prod _{j=1}^D\left(\cdot -z_j\right)\right|\right|\right|}_N\le C_{{\mu }_0,N}\prod _{j=1}^D\left(R_{{\mu }_0}+\left|{\mu }_0-z_j\right|\right)\text{,}$$

$$j^{{\Xi }_S}(Q{)}^2=\sum _{\lambda \in \Lambda }\frac{{\Vert v_{\lambda }\Vert }_V^2}{{\left|{\omega }^{{\Xi }_S}(\lambda )\right|}^2}{\left|Q(\lambda )\right|}^2\text{.}$$

$$j^{{\Xi }_S}(Q_N^{{\Xi }_S})\le \frac{{\left({\sum }_{\lambda \in \Lambda \setminus {\Lambda }_0}{\Vert v_{\lambda }\Vert }_V^2\right)}^{1/2}}{{\left|\left|\left|{\omega }^{{\Lambda }_0}\right|\right|\right|}_N}\underset{\Lambda \setminus {\Lambda }_0}{\sup }\left|\frac{{\omega }^{{\Lambda }_0}}{{\omega }^{{\Xi }_S}}\right|\text{.}$$

$$j^{{\Xi }_S}(Q_N^{{\Xi }_S})\le (1+\epsilon {)}^S\frac{{\left({\sum }_{j>N}{\Vert v_{{\lambda }_j}\Vert }_V^2\right)}^{1/2}}{{\left|\left|\left|{\prod }_{j=1}^N\left(\cdot -{\lambda }_j\right)\right|\right|\right|}_N}\frac{{\prod }_{j=1}^N\left|{\lambda }_{N+1}-{\lambda }_j\right|}{{\Phi }_K({\lambda }_{N+1}{)}^S}\text{for }S>S_{\epsilon ,N}\text{,}$$

$$I^{{\Xi }_S}(\varphi )(\mu )=\sum _{{\mu }^{\prime }\in {\Xi }_S}\frac{\varphi ({\mu }^{\prime }){\omega }^{{\Xi }_S}(\mu )}{(\mu -{\mu }^{\prime }){{\omega }^{{\Xi }_S}}^{\prime }({\mu }^{\prime })}\text{.}$$

$$\frac{1}{(S-1)!}\frac{d^{S-1}}{d{\mu }^{S-1}}{I}^{{\Xi }_S}(uQ)=$$

$$=$$

$$j^{{\Xi }_S}(Q_N^{{\Xi }_S}{)}^2\le j^{{\Xi }_S}{\left(\frac{{\omega }^{{\Lambda }_0}}{{\left|\left|\left|{\omega }^{{\Lambda }_0}\right|\right|\right|}_N}\right)}^2=\sum _{\lambda \in \Lambda \setminus {\Lambda }_0}\frac{{\Vert v_{\lambda }\Vert }_V^2}{{\left|{\omega }^{{\Xi }_S}(\lambda )\right|}^2}\frac{{\left|{\omega }^{{\Lambda }_0}(\lambda )\right|}^2}{{\left|\left|\left|{\omega }^{{\Lambda }_0}\right|\right|\right|}_N^2}\text{,}$$

$$\underset{\Lambda \setminus {\Lambda }_0}{\sup }{r}^{{\Lambda }_0,{\Xi }_S}=\underset{\{{\lambda }_S^{\prime }{\}}_{S>N}}{\sup }{r}^{{\Lambda }_0,{\Xi }_S}\le (1+\epsilon {)}^S\underset{\{{\lambda }_S^{\prime }{\}}_{S>N}}{\sup }\frac{|{\omega }^{{\Lambda }_0}|}{{\Phi }_K^S}\le (1+\epsilon {)}^S\underset{\Lambda \setminus {\Lambda }_0}{\sup }\frac{|{\omega }^{{\Lambda }_0}|}{{\Phi }_K^S}\text{.}$$

$$\left|Q_N^{{\Xi }_S}\left({\lambda }_j\right)\right|\le C_j(1+\epsilon {)}^{4S}{\left(\frac{{\Phi }_K({\lambda }_j)}{{\Phi }_K({\lambda }_{N+1})}\right)}^{2S}\text{for }S>S_{\epsilon ,N}\text{,}$$

$${\mathbb{P}}_N\left(\mathbb{C}\right)\ni Q\leftrightsquigarrow \mathbf{q}\in {\mathbb{C}}^{N+1}\text{implies}{\mathbf{\omega }}_{\lambda }^H\mathbf{q}=Q(\lambda )\text{for }\lambda \in \Lambda \text{.}$$

Full text

Keywords: Model order reduction, rational approximation, Hilbert space-valued meromorphic maps, non-coercive parametric PDEs, Helmholtz equation, frequency response.

AMS Subject Classification: 30D30, 41A20, 41A25, 35J05, 35P15, 65D15.

1 Introduction

Mathematical models based on PDEs are employed to analyze numerically a wide array of physical, financial, and engineering-related phenomena. In many situations, such models depend on one or more parameters, either because of uncertainties or for control/design purposes, and multiple solves of the model have to be performed for different parameter values.

Often the “naive” approach of solving the original problem (which we will refer to as full order model, FOM) as many times as required is not feasible, for instance because the number of evaluations is too large, or because a real-time solution is required. In such cases, model order reduction (MOR) is employed to reduce the computational effort needed for the solution of the FOM. This is achieved by constructing a surrogate model, whose solution is managed to be close to that of the FOM. The construction of the reduced model often requires a considerable computational cost, which, however, can be done offline in a preliminary phase, whereas its evaluation at any parameter value is quite cheap, and can be carried out online in real-time.

Many MOR techniques have been proposed for general FOMs, with the most notable and widely applied being the Reduced Basis (RB) method [15, 16, 17, 21, 26, 35, 36]. The RB method, in its simplest form, assumes that few samples of the solution of the FOM (snapshots) are enough to capture the main features of the solution manifold, i.e. the set of all solutions of the FOM for all parameter values. Exploiting this idea, a surrogate model is built by projecting the FOM onto a subspace generated by few pre-computed snapshots. The effectiveness of the RB method has been verified in many cases. However, it is possible to find some parametric problems for which projection-based MOR does not perform optimally. Such examples are often characterized by some irregularity of the FOM (e.g. nonsmooth dependence on the parameters, non-linearity, lack of stability, etc. of the differential operator therein).

One simple case falling into this category, which is of primary importance in our discussion, is the time-harmonic wave (Helmholtz) equation with parametric wavenumber:

[TABLE]

where $V$ and $W\supset V$ are some suitable Banach function spaces. Depending on the choice of $K$ and of the boundary conditions which complement Eq. 1, the problem above may lack inf-sup stability over $K$ [32], due to the presence of resonant wavenumbers. If one applies carelessly the RB method (in both its main versions, POD and greedy [35, 36]) to Eq. 1, one may observe the appearance of spurious (“non-physical”) resonances in the surrogate model.

Discussions concerning this effect can be found in works related to MOR methods for dynamical systems (e.g. Krylov subspace methods [20, 22]), where problems of similar form are quite common, due to the need for frequency-domain computations. Usually, if the number of parameters is small (often, just 1, the frequency), explicit rather than implicit approximants are considered, with the surrogate model being built by enforcing interpolation or moment matching conditions at the sample points in the parameter space, possibly in a Least Squares (LS) way. The main representatives of this class of methods are the Loewner framework [1, 2, 6, 7, 29] and the vector fitting (VF) algorithm [18, 22, 23, 25, 30, 39], which are tailored to the specific structure of the FOM: rational functions are employed, with the objective of representing each resonance of the FOM by a pole of the approximant. This turns out to be particularly useful when the resonances of the FOM have a physical meaning, since their approximation is implicitly enclosed in the surrogate model, from which it can be obtained with small or no computational effort.

More recently, in the wake of these methods for dynamical systems, and somehow trying to profit from their main advantages, univariate LS Padé approximants have been introduced and studied, both in a standard [9, 11] and a fast [10] version, in the context of a single parameter. Such techniques are based on multiple solves of the FOM at a single parameter value, in the same spirit as Krylov subspace methods, but yield an explicit rational approximant like VF. Actually, LS Padé approximants are actually quite comparable to VF, the main differences between the methods being:

1. (i)

in VF, one usually must achieve a reasonably high “sampling density” [18], namely use at least $2N$ samples to approximate $N$ resonances; conversely, fast LS Padé approximants can get away with half as many samples; 2. (ii)

the fast version of LS Padé approximants relies on the high-dimensionality of the samples, namely on their linear independence, whereas VF can be applied even for scalar outputs (this is not actually a limitation of LS Padé approximants, since it suffices to build a rational approximation of the solution map, and then to apply the desired linear functional to such surrogate solution); 3. (iii)

the error convergence of LS Padé approximants with respect to the number of samples and to the degree of the denominator is well understood for a reasonably broad class of problems [9, 10], while the convergence behavior of the VF error has, to our knowledge, yet to be analyzed thoroughly [18]; 4. (iv)

LS Padé approximants are computed from moments of the solution at a single parameter, the explicit computation of which may feature numerically instabilities, see [12]; instead, VF uses distributed samples in a LS-Lagrange interpolation fashion, for which ill-conditioning issues are much milder.

In this paper we wish to develop and analyze an extension of fast LS Padé approximants which overcomes iv by allowing a distributed parameter sampling, while keeping i and iii intact. In particular, in section 2 we introduce our method, which we name minimal rational interpolation, and we state more precisely to which kinds of FOMs our theory applies to. Then, in section 3 we describe the main convergence results concerning minimal rational interpolation. Afterwards, in section 4 we discuss briefly the topic of a posteriori error/residual estimation for our technique. Possible extensions to multi-parameter problems are described in section 5. In section 6 we show two numerical examples to verify the effectiveness of our method and of the presented error indicator, as well as the theoretical convergence rates. Some concluding remarks are in section 7.

2 Description of the method

Let $\left(V,\langle\cdot,\cdot\rangle_{V}\right)$ be a Hilbert space over $\mathbb{C}$, with induced norm $\left\|\cdot\right\|_{V}$, and $K\subset\mathbb{C}$ be compact. Our task is to find a surrogate for a given map $u:K\to V$. Notably, we seek an approach which builds a surrogate starting only from few evaluations (snapshots) of $u$. Assuming $u$ to be the outcome of a complex computational model, such technique qualifies as non-intrusive, in the sense that it can employ any available solver as a “black box”.

We assume that $S>0$ samples of $u$ at the points $\Xi_{S}=\{\mu_{j}\}_{j=1}^{S}\subset K$ are available. In particular, under some additional regularity assumptions on $u$, we allow confluence of the sample points: if $\mu_{j}$ appears $K_{j}+1$ times in $\Xi_{S}$, we assume to have evaluated $u$, as well as its derivatives with respect to $\mu$ up to order $K_{j}$, at the point $\mu_{j}$.

In order to set up our approximation, it is necessary to specify how many “resonances” of $u$ we wish to approximate. To this aim, we consider the integer $N\in\{0,\ldots,S-1\}$, as well as the set of normalized polynomials of degree at most $N$

[TABLE]

with $\mathbb{P}_{N}\left(\mathbb{C}\right)$ the space of $\mathbb{C}$-valued polynomials of degree up to $N$, and ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{N}$ an arbitrary (but fixed) Hilbertian norm on $\mathbb{P}_{N}\left(\mathbb{C}\right)$. Our aim will be to approximate $N$ resonances of $u$ by the roots of some element of $\mathbb{P}_{N}^{\star}\left(\mathbb{C}\right)$.

Let $I^{\Xi_{S}}$ be the polynomial interpolation operator111For stability reasons it may be wise to replace $I^{\Xi_{S}}$ with a Least Squares-type interpolator. With some effort, our results can be generalized to such case. for $V$-valued functions based on samples at the points $\Xi_{S}$. Namely, given $\phi:\Xi_{S}\to V$, the interpolant $I^{\Xi_{S}}(\phi)$ is the element of the space of $V$-valued polynomials

[TABLE]

such that $\phi(\mu^{\prime})=I^{\Xi_{S}}(\phi)(\mu^{\prime})$ for $\mu^{\prime}\in\Xi_{S}$. In case of confluent sample points, such interpolation condition has to be interpreted in a Hermite sense, with derivatives of the interpolant matching those of the original function.

Definition 1 (Minimal rational interpolants).

A minimal rational interpolant of $u$ of type $[S-1/N]$ with samples at $\Xi_{S}$ is a rational function

[TABLE]

where the denominator $Q_{N}^{\Xi_{S}}$ minimizes the functional

[TABLE]

*over $\mathbb{P}_{N}^{\star}\left(\mathbb{C}\right)$. *

The specific expression of $j^{\Xi_{S}}$ is the natural generalization of the corresponding functional in fast LS Padé approximation [10]. We will see that this form is essential to prove all our convergence results, starting from the central results described in Lemma 6.

We remark that, since $j^{\Xi_{S}}$ is convex, see Eq. 5, the existence of a minimizer $Q_{N}^{\Xi_{S}}$ is guaranteed by the compactness of $\mathbb{P}_{N}^{\star}\left(\mathbb{C}\right)$. From a more practical perspective, the minimization of $j^{\Xi_{S}}$ can be carried out as follows:

1. (i)

build a $V$-orthonormal basis $\left\{\varphi_{i}\right\}_{i=1}^{S}$ of the span of the $S$ samples $\left\{u(\mu_{j})\right\}_{j=1}^{S}$, as well as the corresponding “component map”

[TABLE]

in practice, the basis and the component map can be computed via the Gram-Schmidt procedure or by a generalized QR decomposition of the snapshots; 2. (ii)

fix a basis $\left\{\psi_{l}\right\}_{l=0}^{N}$ of $\mathbb{P}_{N}\left(\mathbb{C}\right)$, orthonormal with respect to the scalar product associated to ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{N}$, and employ it to represent any arbitrary $Q\in\mathbb{P}_{N}^{\star}\left(\mathbb{C}\right)$ by its component vector $\mathbf{q}$ (by construction, $\left\|\mathbf{q}\right\|_{2}=1$); for a more detailed discussion on practical choices of the polynomial basis, we refer to subsection 3.1; 3. (iii)

build the Gramian factor

[TABLE]

which, by linearity, allows to express the target functional as a quadratic form:

[TABLE]

($\overline{\cdot}$ and $\cdot^{H}$ denote conjugation and conjugate transposition, respectively); 4. (iv)

by SVD, obtain the required minimizer, as the minimal right singular vector of the Gramian factor $\Psi$.

We remark that introducing the orthogonal bases $\left\{\varphi_{i}\right\}_{i=1}^{S}$ and $\left\{\psi_{l}\right\}_{l=0}^{N}$ is not strictly necessary for the minimization of $j^{\Xi_{S}}$. However, they allow to state the SVD problem in the usual Euclidean norm, and to reduce it to size $S$, independently of the dimension of $V$.

2.1 Parametric problem framework

While the map $u$ introduced above may be a very general Hilbert-space valued function, we are often interested in approximating maps of a quite specific type. Indeed, assume that we are given a parametric problem:

[TABLE]

with $\mathcal{F}_{\mu}:V\supset\operatorname{Dom}(\mathcal{F}_{\mu})\to V$ a family of densely defined operators, depending continuously on $\mu$. Provided the set $K^{\star}\subset K$ of values of $\mu$, for which a unique $u(\mu)$ solving Eq. 6 exists, is non-empty, we focus on the solution map associated to Eq. 6:

[TABLE]

In the following, we restrict our analysis to a special class of solution maps, namely meromorphic $V$-valued functions over $\mathbb{C}$ with simple poles:

[TABLE]

We denote by $\{v_{\lambda}\}_{\lambda\in\Lambda}\subset V$ the set of residues, which we assume to be $V$-square-summable: $\sum_{\lambda\in\Lambda}\left\|v_{\lambda}\right\|_{V}^{2}<\infty$. Moreover, we assume that the set $\Lambda$ of poles, or resonances, is countable, with no finite limit point. If $\Lambda$ is not finite, Eq. 8 has to be understood in the sense of $V$-convergence, for an arbitrary ordering of $\Lambda$. Without loss of generality, we will assume that $\Lambda$ is unbounded: if this is not the case, one can always extend $\Lambda$ by adding (countably) infinitely many fictitious poles at $\infty$ with arbitrary corresponding residues.

The bulk of our results will be obtained under the additional assumption that the residues form an orthogonal family, i.e. that $\langle v_{\lambda},v_{\lambda^{\prime}}\rangle_{V}=0$ whenever $\lambda\neq\lambda^{\prime}$. This allows considerable simplifications in our derivation: in particular, orthogonal residues have the great advantage of introducing no “correlation” (in the $V$-topology) between components of the solution corresponding to different resonances. Still, many of our conclusions can be generalized to the non-orthogonal case: we defer a discussion of such extensions to subsection 3.4.

To motivate our interest in this class of functions, we give here some examples of operator families $\{\mathcal{F}_{\mu}\}_{\mu\in K}$ which lead to meromorphic solution maps; for their approximation, it is reasonable to employ rational functions, as minimal rational interpolants do. However, we remark that the theory that we develop might not necessarily apply to all the cases which we list here.

• •

Frequency-domain dynamical systems: $\mathcal{F}_{\mu}(u)=(A-\mu E)u+B$, with $A,E\in\mathbb{C}^{n\times n}$, and $B\in V=\mathbb{C}^{n\times n_{I}}$ (one column per input). By employing a determinantal representation [5], one can easily show that the solution map is a rational function; in particular, if the pencil $(A,E)$ admits $n$ (not necessarily finite) eigenvalues counting multiplicity, then $u$ is of the form Eq. 8 [4]. Depending on the spectral properties of $A$, one may be able to ensure residue orthogonality: for instance, it suffices for $E^{\dagger}A$ to be normal, where † denotes the Moore-Penrose pseudo-inverse.

• •

Polynomial eigenvalue-like problems: $\mathcal{F}_{\mu}(u)=A(\mu)u+B$, where $A(\mu)\in\mathbb{C}^{n\times n}$ is a polynomial in $\mu$, and $B\in V=\mathbb{C}^{n\times n_{I}}$. Using augmentation-based approaches [4, 24], one can cast the problem as a frequency-domain dynamical system of size $nM$, where $M$ is the degree of $A$. This implies that the solution map of the original problem (of size $n$) is a rational function; in order to check whether it is actually of the form Eq. 8, one has to investigate the spectral properties of the augmented $(nM)\times(nM)$ problem.

• •

Helmholtz-type problems: $\mathcal{F}_{\mu}(u)=(\mathcal{L}-\mu\mathcal{I})u-f$, where $\mathcal{L}:V\supset\operatorname{Dom}(\mathcal{L})\to V$ is densely defined and bijective, with compact and normal inverse, and $f\in V$ (e.g. $V=L^{2}(\Omega)$, $\mathcal{L}=-\Delta$, and $\operatorname{Dom}(\mathcal{L})=\Delta^{-1}(V)$). The solution map is of the form Eq. 8, and the normality of the resolvent $\mathcal{L}^{-1}$ ensures orthogonal residues [10].

• •

Analytic compact perturbations of bounded bijective operators: $\mathcal{F}_{\mu}(u)=(\mathcal{J}-\mathcal{L}(\mu))u-f$, where $\mathcal{J}:V\to V$ is bounded and bijective, $\mathcal{L}(\mu):V\to V$ is compact for all $\mu\in\mathbb{C}$, and analytic with respect to $\mu$, and $f\in V$; a proof of meromorphy is in [37]. Unfortunately, the form of the solution map is, in general, not (8): for instance, it cannot be guaranteed that the poles are simple, i.e. a $\lambda$-dependent exponent may appear in the denominator of (8).

3 Convergence theory

In this section, we investigate the convergence properties of minimal rational interpolants as the number of samples $S$ increases. To this aim, it becomes necessary to introduce some assumptions on the asymptotic properties of the sample set $\Xi_{S}$. In order to streamline the derivation, we choose to follow a potential theory-inspired approach.

3.1 Preliminaries for convergence theory

Fix $K\subset\mathbb{C}$ either as the closure of the finite region bounded by a Jordan curve or a line segment with positive length. There exist [38] a (unique) positive real number $\operatorname{Cap}(K)$ and a continuous bijective conformal mapping

[TABLE]

holomorphic outside $K$, such that $\left|\phi_{K}\right|=\operatorname{Cap}(K)$ over $\partial K$ and $\lim_{\mu\to\infty}\left|\phi_{K}^{\prime}(\mu)\right|=1$.

The value $\operatorname{Cap}(K)$ is the logarithmic capacity of $K$. It is a powerful concept in complex analysis, algebraic geometry, and approximation theory, and we refer to the monograph [27] for a thorough introduction to the subject. In its more general framework, logarithmic capacity is well defined for any compact set in the complex plane, and coincides with the definition above under our assumptions.

Among its properties we remark some bounds for other metrics over complex sets: for any compact set $A\subset\mathbb{C}$ we have

[TABLE]

with $\ell$ denoting the Lebesgue measure over $\mathbb{C}$. We remark that zero-measure sets are not guaranteed to have capacity zero: for instance, $\operatorname{Cap}([0,4])=1$.

A crucial quantity for stating convergence results in approximation theory is the Green’s potential, which is closely related to the complex magnitude of $\phi_{K}$:

[TABLE]

In particular, we use the Green’s potential to induce an ordering in the set of poles $\Lambda=\{\lambda_{j}\}_{j=1,2,\ldots}$:

[TABLE]

Moreover, we require that $\Phi_{K}(\lambda_{N+1})>\operatorname{Cap}(K)$, i.e. that the number of poles within $K$ is at most $N$, see Eq. 10 and Eq. 11. Indeed, we will show in the following that minimal rational interpolants fail to differentiate between poles of $u$ within $K$; hence, without the assumption $\Phi_{K}(\lambda_{N+1})>\operatorname{Cap}(K)$, the “unidentifiable” poles would distort the approximant and, in general, prevent convergence to the target solution map.

Now, for a given finite point set $\Xi\subset\mathbb{C}$, we define the corresponding nodal polynomial as $\omega^{\Xi}(\mu)=\prod_{\mu^{\prime}\in\Xi}\left(\mu-\mu^{\prime}\right)$. This allows us to state the result which motivates our interest in Green’s potentials.

Theorem 2 (Fejér points [38]).

Let $K$ be either the closure of the finite region bounded by a Jordan curve or a line segment with positive length. There exists a sequence of sample sets \big{(}\Xi_{S}^{K}\big{)}_{S=1,2,\ldots}, whose elements are referred to as the Fejér points for $K$, such that

[TABLE]

This means that the Green’s potential provides an asymptotic upperbound for the magnitude of the nodal polynomial, as long as the sample points are chosen in an optimal way. From here onward we assume that, given $S\in\mathbb{N}$, the sample set $\Xi_{S}$ coincides with some set of Fejér points $\Xi_{S}^{K}$, allowing us to drop the superscript.

Remark 3*.*

Fejér points are explicitly known only for very simple sets $K$; this is the case, for instance, in two situations which are very common in applications: the Fejér points of the unit disk are the roots of unity, and for a line segment they are the Chebyshev nodes [38]. However, our theory applies even for more general sampling strategies. The technical requirement is the following: there exists $\Phi:\mathbb{C}\to\mathbb{R}$ such that

• •

$\sup_{K}\Phi\leq\inf_{\mathbb{C}\setminus K}\Phi$,

• •

$\limsup_{S\to\infty}\left|\omega^{\Xi_{S}}\right|^{1/S}\leq\Phi$ uniformly over $K$,

• •

$\lim_{S\to\infty}\left|\omega^{\Xi_{S}}\right|^{1/S}=\Phi$ uniformly over compact subsets of $\mathbb{C}\setminus K$.

If this is the case, then one can replace the Green’s potential by $\Phi$ in the entirety of our discussion; in particular, $\Phi$ should induce the pole ordering, cf. Eq. 11.

For instance, one could take any sequence which approximates asymptotically the Fejér points [13]. Another interesting case is that of repeated points, employing Lagrange-Hermite interpolation: given $\{\mu_{i}^{0}\}_{i=1}^{L}\subset\mathbb{C}$, one defines the sample points as $\mu_{nL+i}=\mu_{i}^{0}$ for all $n=0,1,\ldots$ and $i=1,\ldots,L$. Then, as long as $K$ is defined as the discrete set $K=\{\mu_{i}^{0}\}_{i=1}^{L}$, the three conditions above hold with $\Phi(\mu)=\prod_{i=1}^{L}\left|\mu-\mu_{i}^{0}\right|^{1/L}$. In particular, this allows to analyze fast LS Padé approximants [10] as a special case of minimal rational interpolants, with $L=1$ and $\Phi(\mu)=\left|\mu-\mu_{1}^{0}\right|$.

A second class of considerations involves the denominator space $\mathbb{P}_{N}^{\star}\left(\mathbb{C}\right)$, whose definition hinges upon the choice of the norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{N}$. Depending on the type of result we wish to prove, we may need to introduce the following assumption on the scaling of ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{N}$ with respect to $N$.

Assumption 4*.*

There exist $\mu_{0}\in\mathbb{C}$ and three positive constants $R_{\mu_{0}}$, $c_{\mu_{0}}$, and $C_{\mu_{0}}$, all independent of $N$, such that

[TABLE]

for all $N\in\mathbb{N}$ and $\{z_{j}\}_{j=1}^{D}\in\mathbb{C}^{D}$, with $0\leq D\leq N$.

In [10] it was shown that Assumption 4 holds (with $R_{\mu_{0}}=c_{\mu_{0}}=C_{\mu_{0}}=1$) for the “monomial” norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{N}={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mu_{0},N}$ defined as

[TABLE]

We describe in Appendix A a broader framework where Assumption 4 can still be shown to hold, with the (often ill-conditioned) monomials in Eq. 14 being replaced by a more general polynomial basis.

Many of our results do not rely on Assumption 4, but exploit the following weaker result.

Lemma 5.

For all $\mu_{0}\in\mathbb{C}$, there exist positive constants $R_{\mu_{0}}$ (independent of $N$), $c_{\mu_{0},N}$, and $C_{\mu_{0},N}$, such that

[TABLE]

*for all $N\in\mathbb{N}$ and $\{z_{j}\}_{j=1}^{D}\in\mathbb{C}^{D}$, with $0\leq D\leq N$. *

Proof*.*

As described above, ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{N}={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mu_{0},N}$ satisfies Lemma 5. The result follows by the equivalence of all norms on a finite-dimensional space.

3.2 Convergence of minimal rational interpolant poles

In the present section we prove that the denominator $Q_{N}^{\Xi_{S}}$ described in Definition 1 performs well in identifying the approximate position of some of the poles $\Lambda$ of the solution map, among which those within the parameter domain $K$. As a preliminary step, it is useful to exploit the definition of $Q_{N}^{\Xi_{S}}$ to derive a bound for the minimal value of $j^{\Xi_{S}}$ over $\mathbb{P}_{N}^{\star}\left(\mathbb{C}\right)$.

Lemma 6.

Let $Q\in\mathbb{P}_{N}\left(\mathbb{C}\right)$. Then

[TABLE]

Moreover, for any subset $\Lambda_{0}$ of $\Lambda$ with cardinality at most $N$, we have

[TABLE]

In particular, for any $\varepsilon>0$, there exists $S_{\varepsilon,N}$ such that

[TABLE]

*provided $\Phi_{K}(\lambda_{N+1})<\Phi_{K}(\lambda_{N+2})$. (If this is not true, there exists a reordering of the poles, still satisfying Eq. 11, for which Eq. 16 holds.) *

Proof*.*

The interpolation operator can be expressed in barycentric coordinates [8] as

[TABLE]

Accordingly, since $Q$ is interpolated exactly,

[TABLE]

The first claim follows by orthonormality of $\{v_{\lambda}\}_{\lambda\in\Lambda}$. From here, we remark that, by construction, $j^{\Xi_{S}}(Q_{N}^{\Xi_{S}})\leq j^{\Xi_{S}}\left(Q\right)$ for all $Q\in\mathbb{P}_{N}^{\star}\left(\mathbb{C}\right)$, so that, by choosing $Q=\omega^{\Lambda_{0}}/{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\omega^{\Lambda_{0}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{N}$, we obtain

[TABLE]

from which Eq. 15 follows.

Now, let $\varepsilon>0$ and $N$ be fixed, and define $\Lambda_{0}=\{\lambda_{j}\}_{j=1}^{N}$. It remains to find a bound for the supremum of $r^{\Lambda_{0},\Xi_{S}}(\lambda)=|\omega^{\Lambda_{0}}(\lambda)/\omega^{\Xi_{S}}(\lambda)|$ over $\lambda\in\Lambda\setminus\Lambda_{0}$. To this aim, let $\lambda^{\prime}_{S}$ maximize $r^{\Lambda_{0},\Xi_{S}}$ over $\Lambda\setminus\Lambda_{0}$ (for all $S$, the supremum must be attained since $r^{\Lambda_{0},\Xi_{S}}(\lambda^{\prime})$ converges to 0 as $\left|\lambda^{\prime}\right|$ diverges to infinity). The set of maximizers $\{\lambda^{\prime}_{S}\}_{S>N}$ must be bounded, due to the polynomial degree of $\omega^{\Lambda_{0}}$ being fixed, and smaller than that of $\omega^{\Xi_{S}}$. Thus, by Eq. 13 there exists $S^{\prime}=S^{\prime}(\varepsilon,N,\Lambda,K)$ such that, for $S>S^{\prime}$,

[TABLE]

Let $\Phi_{K}(\lambda_{N+1})=\Phi_{K}(\lambda_{N+\nu})<\Phi_{K}(\lambda_{N+\nu+1})$ for some $\nu\geq 1$, and assume that $\lambda_{N+1}$ maximizes $|\omega^{\Lambda_{0}}(\lambda)|$ over $\lambda\in\{\lambda_{N+k}\}_{k=1}^{\nu}$. Then, since $|\omega^{\Lambda_{0}}|$ is independent of $S$, the last supremum is attained at $\lambda_{N+1}$, provided $S$ is large enough. This yields Eq. 16.

Before proving convergence of the roots of the approximate denominator to poles of the solution map, we show as a preliminary result that $Q_{N}^{\Xi_{S}}$ takes small values close to some of the elements of $\Lambda$.

Lemma 7.

Let $N$ be fixed. For any $\varepsilon>0$, there exists $S_{\varepsilon,N}$ such that

[TABLE]

*for $j=1,\ldots,N$, with $C_{j}$ independent of $S$ and $\varepsilon$. *

Proof*.*

Through some isometry, we identify $\mathbb{P}_{N}\left(\mathbb{C}\right)$ with $\mathbb{C}^{N+1}$ (endowed with the Euclidean norm), so that, in particular, $\mathbb{P}_{N}^{\star}\left(\mathbb{C}\right)$ corresponds to the boundary of the unit sphere $\partial\mathcal{B}_{N+1}(\bm{0},1)\subset\mathbb{C}^{N+1}$.

By duality, there exists a family $\{\bm{\omega}_{\lambda}\}_{\lambda\in\Lambda}\subset\mathbb{C}^{N+1}$ such that

[TABLE]

Let $\mathbf{q}_{N}^{\Xi_{S}}\in\partial\mathcal{B}_{N+1}(\bm{0},1)$ be the identification of $Q_{N}^{\Xi_{S}}$, and consider the Hermitian matrices $G^{\Xi_{S}},G_{N}^{\Xi_{S}}\in\mathbb{C}^{(N+1)\times(N+1)}$ defined as

[TABLE]

From here onward, the proof resembles quite closely that of Lemma 5.4 in [10], being essentially based on deriving a bound for the smallest singular value of $G^{\Xi_{S}}$ and on finding a rank-$N$ decomposition of $G_{N}^{\Xi_{S}}$. The final result that can be obtained is of the form: for $j=1,\ldots,N$, for large $S$,

[TABLE]

with $C_{j}$ independent of $S$ and $\varepsilon$. From here it suffices to observe that $\left|\omega^{\Xi_{S}}(\lambda_{j})\right|$ is bounded by $\left((1+\varepsilon)\Phi_{K}(\lambda_{j})\right)^{S}$, thanks to Eq. 12.

Finally, we are able to show that the poles of minimal rational approximants converge, as $S$ increases, to some of the elements of $\Lambda$, namely the $N$ poles which are the “closest” according to Eq. 11.

Theorem 8.

Let $N$ be fixed, and denote by $\{\lambda_{j}^{\Xi_{S}}\}_{j=1}^{N}$ the roots of $Q_{N}^{\Xi_{S}}$ (in case of deficient degree, we assume the missing roots to be $\infty$). For any $\varepsilon>0$, there exists $S_{\varepsilon,N}$ such that

[TABLE]

*for $j=1,\ldots,N$, with $C_{j}$ independent of $S$ and $\varepsilon$. *

Proof*.*

We start by observing that, by the normalization of $Q_{N}^{\Xi_{S}}$, it must hold

[TABLE]

for any arbitrary $\mu_{0}\in\mathbb{C}$, see Lemma 5.

From here it suffices to apply the same proof as for Theorem 5.5 in [10], with Eq. 19 replacing the corresponding bound for fast LS-Padé denominators. The main steps are: (i) use Eq. 19 and Eq. 21 to show that, for each $j=1,\ldots,N$, there is a sequence of approximate poles $\{\lambda_{i(S)}^{\Xi_{S}}\}_{S>S_{\varepsilon,N}}$ which converges to $\lambda_{j}$; (ii) exploiting the fact that the poles of $u$ are simple, conclude that $Q_{N}^{\Xi_{S}}$ cannot have, for large $S$, two roots too close to each other; (iii) apply Eq. 19 locally to derive Eq. 20.

Going back to our original problem formulation, let $0\leq\nu\leq N$ be such that $\operatorname{Cap}(K)=\Phi_{K}(\lambda_{\nu})<\Phi_{K}(\lambda_{\nu+1})$. Theorem 8 shows that, for fixed $N$ and large enough $S$, minimal rational interpolants approximate well the $\nu$ poles $\{\lambda_{j}\}_{j=1}^{\nu}$ within $K$, with (exactly) one root of $Q_{N}^{\Xi_{S}}$ converging to each of them at rate

[TABLE]

3.2.1 Pole convergence with respect to denominator degree

The statements of Lemma 7 and Theorem 8 may seem somewhat limited due to the fact that the denominator degree $N$ enters the convergence bounds in a complex way. This makes it impossible to determine how the pole approximation error behaves if $N$ is increased together with $S$, namely if it diverges to infinity. In this section we try to solve this shortcoming. The main result we prove is the following.

Theorem 9.

Take two diverging sequences of integers $(N_{k})_{k\in\mathbb{N}}$ and $(S_{k})_{k\in\mathbb{N}}$, increasing and strictly increasing respectively, such that $N_{k}<S_{k}$ for all $k$. Let $\Xi^{(k)}$ be the shorthand for $\Xi_{S_{k}}$, so that $\smash{Q_{N_{k}}^{\Xi^{(k)}}}$ is the denominator of the $\smash{[S_{k}-1/N_{k}]}$ minimal rational interpolant computed from samples of $u$ at $\Xi_{S_{k}}$. We denote by $\smash{\{\lambda_{j}^{(k)}\}_{j=1}^{N_{k}}}$ the roots of $\smash{Q_{N_{k}}^{\Xi^{(k)}}}$. Also, let Assumption 4 be satisfied. For all $j=1,2,\ldots$, we have

[TABLE]

Proof*.*

Since we are interested in asymptotic properties of the poles, in the present proof we employ a slightly different ordering of the resonances $\Lambda=\{\sigma_{j}\}_{j=1,2,\ldots}$, such that

[TABLE]

Let $j\in\mathbb{N}$ be arbitrary. Thanks to Assumption 4, we can obtain the lower bound

[TABLE]

with $\varphi_{j}=\cdot/(R_{\mu_{0}}+\left|\sigma_{j}-\mu_{0}\right|+\cdot)$ non-negative and strictly increasing over $\mathbb{R}^{+}$. On the other hand, due to Assumption 4 and Lemma 6 with $\Lambda_{0}=\{\sigma_{i}\}_{i=1}^{N_{k}}$, one has

[TABLE]

The two bounds above can be combined with the monotonicity of $\varphi_{j}$ to obtain

[TABLE]

with $C_{j}^{\prime}=C_{\mu_{0}}(\sum_{\lambda\in\Lambda}\left\|v_{\lambda}\right\|_{V}^{2})^{1/2}/(c_{\mu_{0}}\left\|v_{\sigma_{j}}\right\|_{V})$.

Thanks to the bijectivity of $\varphi_{j}$, which maps 0 to itself, the claim follows if we can show that the right hand side of Eq. 24 converges to 0 as $k$ increases. To this aim, a first step is to remove the supremum, whose argument can be simplified by introducing the bounds

[TABLE]

and

[TABLE]

where the term $\left|\mu_{0}-\mu^{\prime}\right|$ is bounded from above by $R_{\mu_{0},K}=\max_{K}\left|\,\cdot-\mu_{0}\right|$.

Now, let $k$ be large enough, so that $\left|\mu_{0}-\sigma_{N_{k}+1}\right|>R_{\mu_{0},K}$. Then the absolute value in the last bound can be omitted, and it must hold

[TABLE]

with

[TABLE]

Going back to the supremum, it is easy to see that all the terms depending on $l$ in the right hand side of Eq. 25 are maximized (over $l>N_{k}$) by the choice $l=N_{k}+1$, thanks to Eq. 23. This, after a suitable rearrangement of the terms, leads to

[TABLE]

All of the factors appearing in Eq. 26, except for the last one, can be easily shown to be bounded as $k$ (and consequently $N_{k}$) goes to $\infty$. Hence, it suffices to show that

[TABLE]

To this aim, the Stolz-Cesàro theorem [3] can be applied to prove that

[TABLE]

yielding the claim.

We have thus shown that, by minimal rational interpolation, we can approximate an arbitrary number of poles of $u$, as long as we allow the denominator degree $N$ to become sufficiently large. However, this more general result does not provide any information on the rate of convergence, at least not a meaningful one for typical applications. Indeed, most of the steps in the proof above introduce gross simplifications, which make it impossible for the intermediate bounds Eq. 24 and Eq. 26 to be reasonably sharp.

3.3 Convergence of minimal rational interpolants

As a complement to the previous results, in this section we investigate the approximation properties of the whole minimal rational interpolant $u_{N}^{\Xi_{S}}$, the reference being the map $u$. Similarly to the analysis of the denominator $Q_{N}^{\Xi_{S}}$, we start by fixing the denominator degree $N$, and by proving convergence with respect to the number of samples $S$.

Theorem 10.

Fix $N$ and $\varepsilon>0$, and let $E$ be an arbitrary compact subset of $\{\mu\in\mathbb{C},\Phi_{K}(\mu)<\Phi_{K}(\lambda_{N+1})\}\setminus\Lambda$. There exists $S_{\varepsilon,N,E}$ such that

[TABLE]

*with $C_{E}$ independent of $S$ and $\varepsilon$. *

Proof*.*

We start by deriving a useful identity for the error $u-u_{N}^{\Xi_{S}}$: since $Q_{N}^{\Xi_{S}}$ is interpolated exactly by $I^{\Xi_{S}}$, we can exploit identity Eq. 17 to obtain

[TABLE]

After dividing by $Q_{N}^{\Xi_{S}}(\mu)$, taking the norm, and applying Lemma 6, we obtain

[TABLE]

for $S$ large enough. In particular, the constant $C^{\prime}$ is independent of $S$ and $\varepsilon$.

Now, due to Eq. 12, the term $\left|\omega^{\Xi_{S}}\right|$ is bounded by $(1+\varepsilon)^{S}\Phi_{K}^{S}$ for $S$ large enough. Thus, it just remains to show that $|Q_{N}^{\Xi_{S}}(\mu)|\min_{\Lambda}\left|\,\cdot-\mu\right|$ is bounded away from 0 for $\mu\in E$. Since $E$ is a compact subset of $\mathbb{C}\setminus\Lambda$, the second factor is not troublesome. The first term is more problematic, since $Q_{N}^{\Xi_{S}}$ may have roots inside $E$. However, Theorem 8 and the triangular inequality ensure that, for $S$ large enough, each root of $Q_{N}^{\Xi_{S}}$ has distance from $E$ bounded away from 0. This yields the claim.

If we restrict our interest to compact subsets $E$ of $K\setminus\Lambda$, we can obtain the simplified result

[TABLE]

Thus, the rate of convergence of the approximation error is the same one that we could expect from polynomial approximants for a function holomorphic over $\{\mu\in\mathbb{C},\Phi_{K}(\mu)<\Phi_{K}(\lambda_{N+1})\}$, having a pole or singularity at $\lambda_{N+1}$.

3.3.1 Error convergence with respect to denominator degree

Just like subsection 3.2.1 describes an extension for variable $N$ of the convergence theory for approximate poles, here we consider a generalization of the error convergence theory in the case of diverging denominator degree. The result we are able to show is of a flavor similar to Theorem 9, and is summarized in the following.

Theorem 11.

Let $(N_{k})_{k\in\mathbb{N}}$, $(S_{k})_{k\in\mathbb{N}}$, and $\Xi^{(k)}$ be as in Theorem 9. Fix $E\subset\mathbb{C}$ compact, and let Assumption 4 be satisfied for some fixed $\mu_{0}\in\mathbb{C}$. Then, for all $\varepsilon>0$, we have

[TABLE]

*with $\operatorname{Cap}$ denoting the logarithmic capacity, see subsection 3.1. *

Proof*.*

We rely on some concepts introduced in the proofs of Theorem 9 and Theorem 10, namely the pole ordering Eq. 23 and the intermediate bound Eq. 29. To start, we follow the same steps as in the proof of Theorem 9, exploiting Eq. 29 and Eq. 15 with $\Lambda_{0}=\{\sigma_{j}\}_{j=1}^{N_{k}}$. The resulting bound is

[TABLE]

with $R=\max_{E}\left|\,\cdot-\mu_{0}\right|$ and

[TABLE]

In order to obtain the desired result, we need to manage carefully the two troublesome $\mu$-dependent terms in the denominator of Eq. 31. To this aim, let $\Lambda^{(k)}=\smash{\{\lambda_{j}^{(k)}\}_{j=1}^{N_{k}}}$ be the set of roots of $Q_{N_{k}}^{\Xi^{(k)}}$. We partition $\Lambda^{(k)}$ into two (potentially empty) sets $\Lambda_{in}^{(k)}$ and $\Lambda_{out}^{(k)}$ according to the characterization

[TABLE]

Given $\#\Lambda_{in}^{(k)}$ and $\overline{N}$ the cardinalities of $\Lambda_{in}^{(k)}$ and $\Lambda\cap E$ respectively, we define the family of lemniscates

[TABLE]

which depends on the index $k$ and on the (small) value $\delta>0$.

We remark that the exponent of $\delta$ in Eq. 32 is equal to the degree of the monic polynomial (in $\mu$) whose magnitude is the left hand side of the inequality. Thus [5, Theorem 6.6.3], the logarithmic capacity of $\mathcal{E}_{k,\delta}$ is equal to $\delta$. (We are assuming without loss of generality that $\smash{\#\Lambda_{in}^{(k)}+\overline{N}}>0$. If this is not the case, $\mathcal{E}_{k,\delta}$ is empty for $\delta<1$, and the claim holds quite trivially.)

The main structure of the remainder of the proof is the following:

1. (i)

we define explicitly a sequence $(\delta_{k})_{k=1,2,\ldots}$ converging to 0; 2. (ii)

we show that $\left\|u(\mu)-u_{N_{k}}^{\Xi^{(k)}}(\mu)\right\|_{V}<\varepsilon$ for all $\mu\in E\setminus\mathcal{E}_{k,\delta_{k}}$, for large $k$.

Then the claim follows, since, by inclusion,

[TABLE]

However, before we can proceed with either step, it is necessary to introduce some additional bounds in Eq. 31. First, since $Q_{N_{k}}^{\Xi^{(k)}}$ is normalized, by Assumption 4 it must hold

[TABLE]

In the first group of factors, $\left|\mu_{0}-\lambda^{\prime}\right|$ is bounded from above by $2R$. In the second group, since $\left|\mu-\mu_{0}\right|\leq R$ for all $\mu\in E$ by definition, we have

[TABLE]

Hence,

[TABLE]

which, provided $\delta\leq R$, for all $\mu\in E\setminus\mathcal{E}_{k,\delta}$, implies

[TABLE]

Accordingly, we can simplify Eq. 31 for all $\mu\in E\setminus\mathcal{E}_{k,\delta}$:

[TABLE]

Since $E$ is compact and $\Lambda$ has no finite limit point, the distance between $E$ and $\Lambda\setminus E$ is strictly positive. Hence, the quantity

[TABLE]

is bounded, for all $\mu\in E$, by a constant $C^{\prime}$ depending only on $E$, $\mu_{0}$, and $\Lambda$.

Now it is trivial to achieve ii by setting

[TABLE]

(We remark that we still need to guarantee that the working assumption $\delta\leq R$ is satisfied. However, since we are planning to prove i, such condition will trivially hold, provided $k$ is large enough.)

At this point, it suffices to apply the same strategy employed in the final part of the proof of Theorem 9. In particular, we have

[TABLE]

while the other factors in Eq. 33 stay bounded as $k$ increases.

Due to Eq. 9, Theorem 11 can be weakened to prove the convergence in probability of $\smash{u_{N_{k}}^{\Xi^{(k)}}}$ to $u$, when both are interpreted as functions from the probability space $(E,\mathcal{B}(E),\ell/\ell(E))$ to the Banach space $(V,\left\|\cdot\right\|_{V})$.

Still, convergence in capacity is somewhat stronger than that in probability: for instance, it ensures that the set of parameter values for which the approximant yields an error above $\varepsilon$ in the limit, namely

[TABLE]

cannot include any curve in $\mathbb{C}$.

3.4 Generalizations to non-orthogonal residues

In the previous sections, we have been working under the assumption that the residues $\{v_{\lambda}\}_{\lambda\in\Lambda}$ form a $V$-orthogonal set. However, most of our results for fixed $N$ can be extended to the general case, in which we impose no condition on the angles between residues. As such, in this section we remove the assumption of residue orthogonality.

We start by extending Lemma 6.

Lemma 12.

Let $Q$, $\Lambda_{0}$, and $\varepsilon$ be the same as in Lemma 6. Then

[TABLE]

Moreover,

[TABLE]

and there exists $S_{\varepsilon,N}$ such that

[TABLE]

*provided $\Phi_{K}(\lambda_{N+1})<\Phi_{K}(\lambda_{N+2})$. *

Proof (sketch)**.

One can follow most of the proof of Lemma 6. However, Eq. 18 does not hold due to the lack of orthogonality. Still, we can replace it by

[TABLE]

which holds by triangular inequality.

We remark the the discrete 2-norm of the residue $V$-norms in Lemma 6 has been replaced by the stronger 1-norm. In order to derive convergence properties, we need to introduce the assumption that the $V$-norms of the residues are summable, so that the inequalities in Lemma 12 are nontrivial.

Theorem 13.

Let $N$, $\{\lambda_{j}^{\Xi_{S}}\}_{j=1}^{N}$, and $\varepsilon$ be the same as in Theorem 8. If $\sum_{\lambda\in\Lambda}\left\|v_{\lambda}\right\|_{V}<\infty$, there exists $S_{\varepsilon,N}$ such that

[TABLE]

and

[TABLE]

*for $j=1,\ldots,N$, with $C_{j}$ independent of $S$ and $\varepsilon$. *

Proof (sketch)**.

The proof is similar to that of Lemma 7 and Theorem 8. However, we can observe that the exponents in the convergence rates are halved. This is caused by “interference” between the modes: the difference between the representative matrices $G^{\Xi_{S}}$ and $G_{N}^{\Xi_{S}}$ is

[TABLE]

which, in particular, contains “off-diagonal” terms, for instance the one corresponding to $\lambda=1$ and $\lambda^{\prime}=N+1$. When computing an upperbound for a term of this form, we cannot hope to obtain more than a single factor of $\Phi_{K}(\lambda_{N+1})^{-S}$, whereas two such factors can be gained in the orthogonal case. Similarly, the correlation between the first $N$ resonances causes one less factor of $\left|\omega^{\Xi_{S}}(\lambda_{j})\right|$ to be obtainable from $G_{N}^{\Xi_{S}}$.

Theorem 14.

Let $N$, $\varepsilon$, and $E$ be as in Theorem 10. If $\sum_{\lambda\in\Lambda}\left\|v_{\lambda}\right\|_{V}<\infty$, there exists $S_{\varepsilon,N,E}$ such that

[TABLE]

*with $C_{E}$ independent of $S$ and $\varepsilon$. *

Proof (sketch)**.

Since Lemma 6 does not hold, it is impossible to apply the middle part of the proof of Theorem 10. Instead, we can split the sum in Eq. 28 into two parts: for the first $N$ poles $\lambda_{1},\ldots,\lambda_{N}$, Theorem 13 can be applied to bound $|Q_{N}^{\Xi_{S}}/\omega^{\Xi_{S}}|$ from above; for the remaining resonances, it suffices to use the argument from the last portion of Lemma 6, with some modifications to account for the presence of $Q_{N}^{\Xi_{S}}$.

Unfortunately, convergence results for variable $N$ do not appear to be possible in the non-orthogonal case, without some strong hypotheses on the “amount of correlation” between resonances: for instance, we believe that one could work under the assumption that the worst-case angle between residues

[TABLE]

is bounded away from 0 for all poles $\lambda$ in the region of interest.

4 A posteriori error indicators

All convergence results presented above are a priori estimates. In particular, they contain quantities which are often not available, namely the number and locations of the poles inside (and of some of the ones outside) the parameter domain $K$. Moreover, the results in Theorems 8 and 11 hold only (with some extensions) if the solution map is of the form Eq. 8, whereas minimal rational interpolants can, in principle, be applied to more general parametric problems and sample points sets.

A typical MOR problem can be cast in the following form: let a certain parameter set $K$ and a parametric problem as in Eq. 6 be given, along with some fixed tolerance $\varepsilon>0$. The task is to obtain an approximate solution map $\widetilde{u}$ such that

[TABLE]

In our framework, the presence of singularities in the true and, possibly, the approximate solution map within $K$ may make the conditions above meaningless. In such cases, it may be preferable [17] to consider a residual-based accuracy requirement

[TABLE]

with the $\left\|\cdot\right\|_{W}$-norm replacing $\left\|\cdot\right\|_{V}$ to account for regularity differences between residual and solution.

Now, let us take as approximate map $\widetilde{u}$ the minimal rational interpolant $u_{N}^{Z_{S}}$ with a certain denominator degree $N$ and samples at the $S$ points $Z_{S}$, which need not be Fejér points for $K$. We face the task of understanding whether the required tolerance $\varepsilon$ is achieved, in the sense of Eq. 34.

In particular, let us consider the problem of computing a posteriori the norm of the residual $\left\|\mathcal{F}_{\mu}\left(\widetilde{u}(\mu)\right)\right\|_{W}$ at some given point $\mu\in K$. In order to proceed, we will assume that the operator $\mathcal{F}_{\mu}$ is linear and has a separable form in $\mu$, i.e. that there exist two families of complex-valued functions $\{\theta_{i}^{F}\}_{i=1}^{n_{F}}$ and $\{\theta_{i}^{f}\}_{i=1}^{n_{f}}$, a family of operators $\{F_{i}\}_{i=1}^{n_{F}}$ over (suitable subsets of) $V$, and a family $\{f_{i}\}_{i=1}^{n_{f}}$ of elements of $V$, such that

[TABLE]

This is the ideal situation to employ RB methods [16, 26, 35, 36], and many strategies have been devised to approximate general parametric problems into the form Eq. 35, e.g. the (D)EIM [15, 21] and hyper-reduction techniques [14].

Now, $F(\mu)$ can be applied to both sides of Eq. 27 to obtain

[TABLE]

with $\Delta^{Z_{S}}(\phi)=\phi-I^{Z_{S}}(\phi)$ the interpolation error for a function $\phi$.

This, after dividing by $\smash{Q_{N}^{Z_{S}}}(\mu)$, provides an affine decomposition of the residual with $n_{F}+n_{f}$ terms. Hence, its norm (squared) admits an affine decomposition with $(n_{F}+n_{f})^{2}$ terms. In particular, assuming an offline-online framework [26, 35, 36], the separable terms of the residual can be precomputed offline, i.e. together with the $S$ expensive evaluations of the target solution map $u$, without increasing the complexity of the overall MOR procedure.

Of course, the discussion above can intrinsically be applied only in an intrusive fashion, since we require access to the operators defining the parametric problem. In general, obtaining an a posteriori non-intrusive error indicator can be quite tricky. Here we propose a simplified approach, which is rigorous only for very specific problem structures, namely

[TABLE]

For such problems, thanks to Eq. 17 and Eq. 36, it must hold

[TABLE]

Now it suffices to take the $\left\|\cdot\right\|_{W}$-norm and exploit Eq. 4 to obtain

[TABLE]

One obvious limitation of Eq. 38 and Eq. 39 is the need to evaluate or bound expensive norms (either the $W$ one or the $V\to W$ operator one) involving $F_{1}$. We remark that this drawback can be interpreted as an issue in identifying the exact scaling of the error estimator, a common problem even for stable parametric problems, where the inf-sup (or coercivity) constant must be estimated to link residual and error [28, 35, 36]. If one can overcome such limitation, see subsection 6.2, Eq. 38 and Eq. 39 are certainly quite appealing from a computational standpoint, since their evaluation consists only of few scalar computations.

If Eq. 37 does not hold, we can replace $F(\mu)$ by a suitable linear approximation (e.g. its first-order Taylor series around some parameter) in a spirit similar to that of the Empirical Interpolation Method. Then Eq. 38 and Eq. 39 can be applied heuristically, with an accuracy which may depend quite sharply on the smoothness of $\mathcal{F}_{\mu}$, in particular on its second-order variations over the parameter domain.

As a supplement to our discussion, we would like to note that a posteriori error indicators are often used in RB procedures not only to determine whether the required accuracy has been achieved, but also to drive the selection of the sample points. For instance, in the weak-greedy RB approach [35, 36], one keeps adding a new snapshot at the parameter value which maximizes the a posteriori error indicator of choice, thus updating the surrogate model (and the error indicator as well), until convergence. We envision that a similar adaptive procedure could be devised also for minimal rational interpolation, with new parameter values being added to the sample set according to the same logic.

5 Extensions to multiple parameters

Throughout our discussion, we have only dealt with the approximation of univariate meromorphic solution maps. However, it is often of interest to approximate the solution map of multi-parameter problems: for example, in the analysis of parametric dynamical systems, one may need to approximate the solution of

[TABLE]

where $\mathbf{p}$ is a collection of physical or geometric parameters. For such problems, it is often quite complicated to prove properties of $u$, jointly in $\mu$ and $\mathbf{p}$; however, the analysis is much simpler if the extra parameters $\mathbf{p}$ are kept fixed (“marginalized out”) and one looks at the univariate solution map $u_{\mathbf{p}}(\mu):\mu\mapsto u(\mu,\mathbf{p})$.

This has lead to the development of parametric MOR (pMOR) [6, 19, 33, 40, 41], whose core idea is the following: take a set of parameter values $\{\mathbf{p}_{1},\ldots,\mathbf{p}_{T}\}\subset\mathbb{C}^{d-1}$ and, for each $j=1,\ldots,T$, build a surrogate $\widetilde{u}_{\mathbf{p}_{j}}$ of $u_{\mathbf{p}_{j}}$ (this involves computing the snapshots $u(\mu_{i},\mathbf{p}_{j})$ for $i=1,\ldots,S_{j}$, and applying a single-parameter MOR strategy, with respect to $\mu$ only); then, the joint reduced model $\widetilde{u}$ of $u$ is constructed as a suitable combination of $\widetilde{u}_{\mathbf{p}_{1}},\ldots,\widetilde{u}_{\mathbf{p}_{T}}$.

In a pMOR framework, one could apply minimal rational approximation to obtain the surrogates in frequency only $\{\widetilde{u}_{\mathbf{p}_{j}}\}_{j=1}^{T}$. In particular, the results presented in the previous sections apply: for instance, the residual estimator in section 4 is still valid, allowing for a greedy approach with respect to $\mu$. This line of approximation for multi-parameter problems is currently under investigation.

There are several reasons to prefer pMOR over global MOR (in $d$ parameters) for problems like Eq. 40. Many of them are related to the curse of dimensionality, which forces to take a large number of snapshots, and to build a surrogate model of considerable size, undermining the cost-effectiveness of MOR. We refer to the general pMOR literature cited above for a broader discussion.

Here, we wish to remark that these disadvantages affect also the natural extension of minimal rational interpolation to $d$ parameters, which we could define as follows: assume that the sets of polynomials $\mathbb{P}_{N}\left(\mathbb{C}^{d};\mathbb{C}\right)$ and $\mathbb{P}_{N}^{\star}\left(\mathbb{C}^{d};\mathbb{C}\right)$, as well as the interpolation operator $I^{\Xi_{S}}$, are suitably defined (for instance, $\mathbb{P}_{N}\left(\mathbb{C}^{d};\mathbb{C}\right)$ could denote the space of polynomials of total degree at most $N$); the optimization problem which yields the surrogate denominator $Q_{N}^{\Xi_{S}}$ is now: find $Q_{N}^{\Xi_{S}}\in\mathbb{P}_{N}^{\star}\left(\mathbb{C}^{d};\mathbb{C}\right)$ which minimizes over $\mathbb{P}_{N}^{\star}\left(\mathbb{C}^{d};\mathbb{C}\right)$ the convex functional

[TABLE]

where $M$ is the maximal total degree of the interpolant $I^{\Xi_{S}}(uQ)$, so that, in general, $S\sim M^{d}$. In practice, an algorithm for minimizing this functional can be found as a generalization of the strategy proposed in section 2; while we omit the details here, we remark that the Gramian matrix $\Psi$ representing $j^{\Xi_{S}}$, see Eq. 5 for the univariate case, becomes of size approximately $S^{2-1/d}\times N^{d}$.

In addition to the other limitations of multivariate minimal rational interpolation, we also observe that it is impossible to extend directly our single-parameter theory to the multivariate case, since many of the core ideas (e.g. the barycentric expansion Eq. 17) do not generalize nicely to more than one dimension. Overall, also considering the numerical issues intrinsic to high-dimensional interpolation, we believe minimal rational approximants to be an appropriate MOR technique only in low dimension: extensions to many parameters should rely on pMOR or similar approaches.

6 Numerical examples

In this section, we perform some numerical experiments to show the usefulness of minimal rational interpolation, and to verify the theoretical convergence rates and error indicator. The first example satisfies the assumption under which our results have been derived; in particular, the residues are orthogonal. Instead, the solution map in the second example can be proven meromorphic, but not necessarily of the form Eq. 8. For the sake of reproducibility, the code used to run the simulations is available at [34].

6.1 Normal eigenvalue problem

Let $n=100$ and take $A\in\mathbb{C}^{n\times n}$ a normal matrix whose eigenvalues are randomly generated according to a uniform distribution over $\{x+\iota y,(x,y)\in[-5,5]^{2}\}$, and whose eigenvectors are prescribed by orthonormalizing a matrix with random (complex) standard gaussian entries. In particular, we order the eigenvalues $\{\lambda_{j}\}_{j=1}^{n}$ according to their distance from 0.

Our task is to estimate the eigenvalues of $A$ which lie within the unit disk $\smash{K=\overline{\mathcal{B}(0,1)}}$. To this aim, we fix $\mathbf{v}\in\mathbb{C}^{n}$ a random gaussian vector, and we approximate by minimal rational interpolants the solution map of the parametric problem

[TABLE]

with $I\in\mathbb{C}^{n\times n}$ the identity matrix. We endow $\mathbb{C}^{n}$ with the Euclidean scalar product. As sample points $\Xi_{S}$, we take the $S$-th roots of unity, and we choose to employ the polynomial norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{0,N}$, see Eq. 14. This norm is induced by monomials, which are $L^{2}(\partial K)$-orthogonal. Moreover, the action of the interpolation operator $\smash{I^{\Xi_{S}}}$ can be evaluated quite efficiently as a Fast Fourier Transform.

Let us assume that our interest lies in approximating just the eigenvalues of $A$ inside $K$. As their exact number is unknown, we can choose the denominator degree $N$ according to different strategies:

1. •

we fix $N$ (in our case, we set $N=10$) and increase progressively the number of samples $S$; to verify if the starting guess on $N$ was large enough, we can simply check if at least one of the approximate poles has converged (for large $S$) to some point outside $K$, see Theorem 8. 2. •

we can increase $N$ along with the number of samples $S$, for instance by choosing $N=S-1$ for all $S$; in this case we are sure to approximate all poles, see Theorem 9.

We compare visually the two approaches in Figure 2, where we show the error norm obtained by approximants of type $[20/10]$ and $[20/20]$, both based on $S=21$ samples at roots of unity. We can see that the error increases much quicker near the boundary of the sampled square in the case $N=10$. This is to be expected, as the region of convergence increases together with $N$, see Theorem 10. A (possibly) more interesting observation is that the error appears to be globally smaller when we choose $N=20$. Overall, we can conclude that the choice $N=S-1$ leads to a better approximant.

Now it just remains to check how well each approach manages to capture the location of the eigenvalues of $A$. To this aim, we increase $S$ from $11$ to $30$, and compute the pole approximation error, defined in the left hand side of Eq. 20, for the 5 poles within $K$. The results are shown in Figure 2. For fixed $N=10$, the rate predicted by Eq. 22 can be observed (Green’s potential for a centered disk coincides with the complex magnitude outside the disk). Remarkably, a similar, if not slightly better, rate of convergence is obtained for increasing $N=S-1$, despite the absence of theoretical results in this regard.

In Figure 2 we also show the pole approximation error obtained through a RB-like projection technique: given the same sample set $\Xi_{S}$, we perform a Proper Orthogonal Decomposition (POD) [35] of the samples, identifying the $N$ dominant “sample directions”; then we project orthogonally the original problem on the subspace spanned by these $N$ directions, and we use the solution of the resulting $N\times N$ eigenvalue problem to approximate the original one. We can observe that the projection technique performs much worse than minimal rational interpolants for fixed $N=10$, whereas for $N=S-1$ the two methods achieve very similar results.

Lastly, we wish to test two other approximation strategies: a fast LS Padé approximant [10] centered at the origin, and a minimal rational interpolant based on quasi-random (Halton) points in the unit disk. We show the error norms and the approximated poles in Figure 4 for $N=S-1=20$. Qualitatively, it appears that, for both choices, the accuracy in the approximation of the resonances is similar to that achieved with samples at the roots of unity. With respect to the $[20/20]$ minimal rational interpolant in Figure 2, the errors for the two new approximants appear smaller within $K$, but larger outside $K$; this is to be expected, since the new sample points are closer to the origin.

More quantitatively, we show in Figure 4 the convergence of the pole approximation error with respect to $S$ (here as well, we keep $N=S-1$). For fast LS Padé approximants, when $S$ is small, the error converges quite quickly (the error is uniformly smaller than the one obtained with samples at the roots of unity), but then it stagnates, and even increases starting from $S=20$; this is a symptom of the instability which affects the construction of this centered approximant, namely the computation of high-order derivatives of a meromorphic function.

The results for quasi-random sample points appear relatively similar, with the error flattening out starting from $S=20$. In this case, the non-monotone behavior of the error for large $S$ is due to the interaction between two conflicting effects: the addition of sample points close to the resonances (which would, on its own, improve the accuracy of the approximation) and the ill-conditioning of the interpolation problem (for $S=30$, the Lebesgue constant is $\sim 10^{7}$, as opposed to $\sim 3$ for the $30$-th roots of unity).

6.2 Time-harmonic vibrations of a tuning fork

Let $\Omega\subset\mathbb{R}^{3}$ be the region occupied by a tuning fork. Assume that the tuning fork is clamped on a portion $\Gamma_{D}$ of the handle. Let $\smash{\widetilde{\mathbf{g}}}(t,\mathbf{x})=e^{\iota 2\pi\overline{\nu}t}\mathbf{g}(\mathbf{x})$ be a pressure pulse (sinusoidal in time and gaussian in space) impinging on a portion $\Gamma_{\mathbf{g}}$ of the exterior of the tuning fork, which we assume to be sound-hard. See Figure 5 for a graphical representation of the setup.

We are interested in computing the time-harmonic deformation $\mathbf{u}(\nu)$ with frequency $\nu$ which the tuning fork undergoes. Assuming the tuning fork to have density $\rho$ and linear stress constitutive relation $\sigma(\cdot)$, the displacement $\mathbf{u}(\nu)$ can be found by solving a (damped) time-harmonic linear elasticity problem

[TABLE]

As $V$-inner product, we choose the elastic energy inner product:

[TABLE]

In particular, we choose to approximate problem Eq. 42 using $\mathbb{P}_{1}$ FE on a sufficiently fine tetrahedral discretization of $\Omega$. The FE discretization of Eq. 42 defines a parametric problem of the form Eq. 6, whose solution map $\mathbf{u}$ can be proven meromorphic using compactness arguments [11, 31]. Thus, we decide to approximate it for $\nu\in K=[10,40]$ kHz with minimal rational interpolants relying on samples at the Chebyshev nodes of $K$. The chosen polynomial norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{N}$ is the one induced by the family of Chebyshev polynomials over $K$, see subsection 3.1.

First, we investigate the behavior of the relative approximation error $\frac{\left\|\widetilde{\mathbf{u}}(\nu)-\mathbf{u}(\nu)\right\|_{V}}{\left\|\mathbf{u}(\nu)\right\|_{V}}$ at $\nu=13$ kHz and $\nu=31$ kHz, for $\eta=0$ Hz. We apply different approximation strategies: for $S=1,\ldots,51$, we consider $[S-1/S-1]$ fast LS Padé approximants [10] with $S$ samples at the center of $K$, $[S-1/N]$ minimal rational approximants for fixed $N=10,20,30$, and diagonal $[S-1/S-1]$ minimal rational approximants.

The results are shown in Figure 8. As expected from Theorem 10, minimal rational approximants of type $[S-1/10]$ do not converge ($K$ contains 16 resonances), whereas the ones of type $[S-1/20]$ and $[S-1/30]$ do; we can observe that the diagonal approximant $[S-1/S-1]$ yields the best results overall.

The error achieved by fast LS Padé approximants converges extremely quickly at $\nu=31$ kHz for small $S$, but no convergence can be observed at $\nu=13$ kHz; this is reasonable, since this type of approximant, relying on information at a single point, can deliver an extremely accurate surrogate close to the center, but the quality of the approximation degrades as we move further away. Also, we observe that the convergence plots for fast LS Padé approximants stop decreasing already at $S=11$, and terminate abruptly at $S=35$: for $S\geq 36$, the construction of the surrogate model fails, due to conditioning issues in the evaluation of the derivatives of the solution map, and in the application of the (Taylor) interpolation operator [10, 12].

The norm of the solution map is shown in Figure 8 for $\eta\in\{0,500\}$ Hz, along with the approximated norm profiles, obtained by minimal rational interpolants of type $[20/20]$ and $[40/40]$. If $S=21$, for both values of $\eta$, we can observe that the approximation of $\left\|\mathbf{u}\right\|_{V}$ is quite accurate except for the low-frequency region, where the approximant for $\eta=0$ Hz actually appears quite unstable. This behavior can be improved by increasing the number of samples to $S=41$.

In a (more) realistic application, in order to determine whether there is actual need for more than 21 samples, one could rely on residual-based estimators such as Eq. 34. To this aim, we can apply Eq. 36 to obtain an expression for the residual in the $L^{2}(\Omega)^{3}$-norm. The evaluation of Eq. 36 on a fine grid is shown in Figure 8.

In the same Figure, we also show the values obtained through the simplified residual estimator Eq. 39. In particular, in order to evaluate the estimator, we need to compute or approximate the operator norm $\left\|F_{1}\right\|_{\mathcal{L}(V;L^{2}(\Omega)^{3})}$, with $F_{1}$ the derivative of the operator in Eq. 42 at some point $\mu^{\prime}$, see section 4. We manage to avoid this, thus keeping the estimator as non-intrusive as possible, by the following argument: if we assume that Eq. 38 holds approximately (the problem is quadratic and not linear in the parameter), we can identify a constant $C$, independent of $\mu$, for which, employing the notation of section 4,

[TABLE]

In practice, the value of $C$ can be easily found by evaluating (offline and at a relatively low computational cost) the exact value of the residual at some point $\mu^{\prime}\in K\setminus Z_{S}$, so that the value of $C$ can be set to

[TABLE]

In our example we choose $\mu^{\prime}=25.5$ kHz, and we can observe the corresponding residual estimator to be extremely accurate. Indeed, the values obtained with the two estimators never differ by more than 1%, and the corresponding curves cannot be distinguished in the plot.

Finally, in Figure 8 we also perform a comparison between minimal rational interpolants and RBs, by showing the norm of the residual obtained with the RB approximant of Eq. 42 computed from samples at $\Xi_{21}$. We can observe that the two residuals are quite similar in terms of both order of magnitude and behavior with respect to $\nu$. In particular, simplifying grossly, minimal rational interpolants seem to have a very slight edge on RBs at higher frequencies, while the opposite appears to be true at lower frequencies.

7 Conclusions

In this paper we have presented minimal rational interpolants as the natural extension of fast LS-Padé approximants [10], and we have shown that most of the theoretical properties of this latter technique generalize nicely to the former. The new method can be considered superior to the original one for two main reasons:

1. •

the stability is greatly improved, since computing derivatives of the solution map accumulates numerical errors; the Arnoldi-like reorthogonalization technique described in [10] can somehow limit but not totally eliminate this issue; 2. •

the distributed sampling helps to achieve a globally (e.g. using the $L^{\infty}(K)$-metric) small approximation error, which is the target in most applications; of course, an interpolatory approach cannot be as accurate close to the samples as a Taylor-type approximation, assuming the number of samples to be the same.

Our numerical example has shown minimal rational interpolants to preform well in the MOR of a normal eigenproblem and of a time-harmonic problem. In particular, we have shown the effectiveness of a heuristic and cheap to compute a posteriori residual estimator. Accordingly, we deem of interest to investigate a possible greedy-type algorithm, in the same spirit as the one for the RB method, but with minimal rational interpolants replacing Galerkin projection. In particular, it is still unclear how well the approximation properties of the method are preserved if the samples are not computed at optimal (Fejér) points.

Acknowledgments

The author gratefully thanks Fabio Nobile for the many fruitful discussions and for his helpful advice at several stages of this work.

Appendix A Polynomial norm bounds

Let $\Gamma\subset\mathbb{C}$ be a Jordan curve, and consider $\{\varphi_{j}\}_{j=0,1,\ldots}$ a hierarchical (degree$(\varphi_{j})=j$ for all $j$) family of polynomials, orthonormal over $L^{2}(\Gamma)$:

[TABLE]

For instance, if $\Gamma=\partial\mathcal{B}(\mu_{0},1)$, we can choose $\varphi_{j}(\mu)=\mu^{j}$ for all $j$.

This, for all $N\in\mathbb{N}$, induces a norm over $\mathbb{P}_{N}\left(\mathbb{C}\right)$ in the following sense:

[TABLE]

The following holds.

Lemma 15.

*The norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{N}$ in Eq. 44 satisfies Assumption 4 for any $\mu_{0}$ chosen within the interior of $\Gamma$. *

Proof*.*

For a given $N\in\mathbb{N}$, let $0\leq D\leq N$, $\{z_{j}\}_{j=1}^{D}\in\mathbb{C}^{D}$, and $Q=\prod_{j=1}^{D}\left(\,\cdot-z_{j}\right)$. Due to orthonormality, it holds

[TABLE]

On one side, let $R_{\mu_{0}}=\max_{\Gamma}\left|\,\cdot-\mu_{0}\right|$. We can use the triangular inequality to obtain

[TABLE]

On the other hand, let $c_{\mu_{0}}=\left(2\pi\min_{\Gamma}\left|\,\cdot-\mu_{0}\right|/\int_{\Gamma}\left|\textup{d}l\right|\right)^{2}$. We can exploit the Cauchy-Schwarz inequality to show that

[TABLE]

so that, due to the Cauchy integral formula,

[TABLE]

Dividing Eq. 46 and Eq. 47 by $\int_{\Gamma}\left|\text{d}l\right|$ yields the proof.