Stability of rootfinding for barycentric Lagrange interpolants (Q2454385)

Given a univariate polynomial defined by \(n+1\) values \([f_0,\dots,f_n]\) at distinct nodes \([x_0,\dots,x_n]\), define the nodal polynomial \(\ell(z)\) by \(\ell(z) = \prod_{i=0}^n(z-x_i)\), and the barycentric weights \(w_j\) by \(w_j = \prod_{k \neq j}(x_j-x_k)^{-1}\), for \(0 \leq j \leq n\). The barycentric interpolation formula is then: \[ p(z) = \ell(z) \sum_{j=0}^{n} \frac{w_jf_j}{(z-x_j)}. \] The roots of polynomials expressed in this basis can be found via the eigenvalues of a generalized companion matrix pair \((\mathbf{A},\mathbf{B})\), and it is often advantageous to do so for reasons of numerical stability. While the stability of this method previously has been verified empirically, the present paper provides a detailed analysis of the stability of this approach. The authors show, in particular, that this method is normwise backward stable, and they propose a bound on the backward error which is ``typically only an order of magnitude larger than the actual error''. The analysis is rigorous as well as practical. For example, the computed eigenvalues of \((\mathbf{A},\mathbf{B})\) are assumed to be the exact eigenvalues of the matrix pair \((\mathbf{A+\Delta A}, \mathbf{B+\Delta B})\), where the perturbations satisfy \[ ||(\mathbf{\Delta A},\mathbf{\Delta B})||_F \leq \sigma(n) \varepsilon_M||(\mathbf{A}, \mathbf{B})||_F \] ``(as guaranteed by LAPACK...)'', where \(\varepsilon_M\) is machine precision, and \(||(\mathbf{A}, \mathbf{B})||_F\) is the Frobenius norm of a matrix pair \((\mathbf{A}, \mathbf{B})\). Further details as to how this is calculated can be found in the LAPACK Users' Guide. This leads them to an explicit bound on, essentially, the ratio of the backward error in the coefficients of the polynomial \(p(z)\) to the backward error of the generalized eighenvalue problem. The authors note that for their experiments, they cho0se \(\sigma(n) = \sqrt{n}\) because it seemed to be a ``good experimental fit, but [they] admit that [they] have no theoretical reason for choosing this.'' The authors give a full discussion of balancing and scaling of the companion matrix to bring it closer to a ``normal pair'', which allows them to produce roots with small backward error. They conclude with numerous experiments, comparing their bound to the test problems given in the paper by \textit{A. Edelman} and \textit{H. Murakami} [Math. Comput. 64, No. 210, 763--776 (1995; Zbl 0833.65041)], Chebyshev polynomials of the first and second kinds, the scaled Wilkinson polynomial and a Wilkinson filter example introduced on Page 178 of the article by \textit{J. H. Wilkinson} [Numer. Math. 1, 150--166, 167--180 (1959; Zbl 0202.43701)].