Finite-sample AAA rational approximation with exact SVD

Input a nonempty ordered list of distinct sample points $Z=(z_1,\dots ,z_N)\subset \mathbb{C}$, data $F=(f_1,\dots ,f_N)\in {\mathbb{C}}^N$, an absolute tolerance $\tau \ge 0$, an integer support budget $B$ with $0\le B\le \lfloor N\mathrm{/}2\rfloor$, and a finite integer bidiagonal-QR sweep budget $S\ge 0$. The procedure is evaluated in exact arithmetic; its finite-precision implementation makes no optimality claim. For $v\in {\mathbb{C}}^k$, use the standard Hermitian inner product and Euclidean norm $\mathrm{\parallel }v{\mathrm{\parallel }}_2=({\sum }_{j=1}^k\mathrm{\mid }{v}_j{\mathrm{\mid }}^2{)}^{1\mathrm{/}2}$. Initialize the accepted support-record list $S_{\mathrm{a}\mathrm{c}\mathrm{c}}=\mathrm{\varnothing }$; each record stores an original index, its sample point, and its data value. Set the remaining original-index list $R=(1,\dots ,N)$. Every return has form $(r_{\mathrm{a}\mathrm{c}\mathrm{c}},S_{\mathrm{a}\mathrm{c}\mathrm{c}},\text{status},\mathcal{R})$, where $\mathcal{R}$ is the residual record maintained below.

1. Set the initial rational candidate $r_0$ to the constant mean $\bar{f}=N^{-1}{\sum }_i{f}_i$. Compute ${\eta }_0={\max }_i\mathrm{\mid }{f}_i-\bar{f}\mathrm{\mid }$, set the last accepted candidate $r_{\mathrm{a}\mathrm{c}\mathrm{c}}=r_0$, and set $\mathcal{R}=((0,{\eta }_0))$. If ${\eta }_0\le \tau$, return $(r_{\mathrm{a}\mathrm{c}\mathrm{c}},S_{\mathrm{a}\mathrm{c}\mathrm{c}},\mathtt{\text{tolerance-met}},\mathcal{R})$.

2. For $m=1,\dots ,B$, choose from $R$ the least original index $s$ that maximizes the absolute residual from $r_{\mathrm{a}\mathrm{c}\mathrm{c}}$. Let $S_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{v}}$ be $S_{\mathrm{a}\mathrm{c}\mathrm{c}}$ with $(s,z_s,f_s)$ appended; remove $s$ from $R$.

3. Form the $(N-m)\times m$ Loewner matrix $L_{ij}=(f_i-f_{s_j})\mathrm{/}(z_i-z_{s_j})$, where $i$ ranges over the indices in $R$ and $s_j$ are the original indices in $S_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{v}}$. The denominator is nonzero because all sample points are distinct.

4. Compute an exact full singular value decomposition (SVD) of $L$ by applying Householder bidiagonalization and its finite precision status followed by Finite budget implicit shifted QR for bidiagonal singular values with zero deflation tolerance and budget $S$. If either call returns failure or the QR budget is exhausted, return $(r_{\mathrm{a}\mathrm{c}\mathrm{c}},S_{\mathrm{a}\mathrm{c}\mathrm{c}},\mathtt{\text{inconclusive}},\mathcal{R})$.

5. Let $E$ be the right singular subspace for the smallest singular value of $L$, and let $P_E$ be its orthogonal projection. Choose the least coordinate index $j$ for which $P_E{e}_j\mathrm{\ne }0$ and set $w=P_E{e}_j\mathrm{/}\mathrm{\parallel }{P}_E{e}_j{\mathrm{\parallel }}_2$. This rule applies at every multiplicity, including dimension one. Multiply by the unique unit scalar that makes the first nonzero coordinate of $w$ positive real. If any coordinate of $w$ is zero, reject this provisional support and return $(r_{\mathrm{a}\mathrm{c}\mathrm{c}},S_{\mathrm{a}\mathrm{c}\mathrm{c}},\mathtt{\text{zero-support-weight}},\mathcal{R})$.

6. Form the barycentric rational quotient using the nodes and data in $S_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{v}}$ and $w$ as its weights. If this rational function has a genuine pole at any remaining sample point, reject this provisional support and return $(r_{\mathrm{a}\mathrm{c}\mathrm{c}},S_{\mathrm{a}\mathrm{c}\mathrm{c}},\mathtt{\text{pole-status}},\mathcal{R})$. Otherwise accept the candidate by setting $r_{\mathrm{a}\mathrm{c}\mathrm{c}}$ to it and $S_{\mathrm{a}\mathrm{c}\mathrm{c}}=S_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{v}}$, and append $(m,{\eta }_m,{\sigma }_m)$, where ${\eta }_m$ is its maximum sampled residual and ${\sigma }_m=\mathrm{\parallel }Lw{\mathrm{\parallel }}_2$, to $\mathcal{R}$. If ${\eta }_m\le \tau$, return $(r_{\mathrm{a}\mathrm{c}\mathrm{c}},S_{\mathrm{a}\mathrm{c}\mathrm{c}},\mathtt{\text{tolerance-met}},\mathcal{R})$.

7. If all $B$ iterations finish without another return, return $(r_{\mathrm{a}\mathrm{c}\mathrm{c}},S_{\mathrm{a}\mathrm{c}\mathrm{c}},\mathtt{\text{support-budget-exhausted}},\mathcal{R})$. In particular $B=0$ returns the initial constant and empty support list with this status when step 1 does not return.

The support tie rule and the projection rule make the exact procedure deterministic, including when the smallest singular value has multiplicity greater than one. Rank deficiency is permitted: if the smallest singular value is zero and all selected weights are nonzero and the quotient is finite at the remaining sample points, the candidate is accepted and its exact sample residual is evaluated. The algorithm asserts no uniform off-grid error bound and no minimax optimality. The Loewner residual ${\sigma }_m$ is a linear least-squares residual, not the nonlinear sampled residual ${\eta }_m$.