Confluent Vandermonde with Arnoldi

The recent work [1] presents a linear algebra method for the polynomial interpolation and fitting on arbitrary grids. It does not need the prior knowledge of a well-conditioned polynomial basis, and works easily on two intervals and subsets of complex plane, for example. The method starts with the ill-conditioned Vandermonde system

$$\sum_{k=0}^{n}x_{j}^{k}c_{k}=f(x_{j}),\quad j=1,\ldots,m$$

such that $f(x)\approx\sum_{k=0}^{n}x^{k}c_{k}$, and transforms it into a much more well-conditioned system

$$\sum_{k=0}^{n}q_{jk}d_{k}=f(x_{j}),\quad j=1,\ldots,m$$

such that $f(s_{j})\approx\sum_{k=0}^{n}w_{jk}d_{k}$. In matrix form, the Vandermonde system $V\mathbf{c}=\mathbf{f}$, i.e.,

$$\begin{pmatrix}1&x_{1}&x_{1}^{2}&\cdots&x_{1}^{n}\\ 1&x_{2}&x_{2}^{2}&\cdots&x_{2}^{n}\\ \vdots&\vdots&\vdots&\cdots&\vdots\\ 1&x_{m}&x_{m}^{2}&\cdots&x_{m}^{n}\end{pmatrix}\begin{pmatrix}c_{0}\\ c_{1}\\ \vdots\\ c_{n}\end{pmatrix}=\begin{pmatrix}f(x_{1})\\ f(x_{2})\\ \vdots\\ f(x_{m})\end{pmatrix}$$

is transformed into $Q\mathbf{d}=\mathbf{f}$ where $Q=(q_{jk})$ is a matrix and $\mathbf{d}=(d_{k})$ is a column vector, and the approximation at $s_{j}$, $j=1,\ldots,M$ is evaluated as $(f(s_{j}))\approx W\mathbf{d}$ where $W=(w_{jk})$. The columns of $Q$ form an orthogonal basis of the column space of $V$. To construct $Q$, it is noted that $\mathrm{Col}\,V$ is just the Krylov subspace $\mathrm{span}\{\mathbf{e},X\mathbf{e},X^{2}\mathbf{e},\ldots,X^{n}\mathbf{e}\}$ with $X=\mathrm{diag}(x_{1},x_{2},\ldots,x_{m})$, $\mathbf{e}=(1,\ldots,1)^{T}$. So, the Arnoldi process for the Krylov subspace is exploited to generate the columns of $Q$ as well as a Hessenberg matrix $H$ such that $XQ_{-}=QH$ (here $Q_{-}$ is obtained by deleting the last column of $Q$). The evaluation matrix $W$ has the first column all one’s and is constructed column by column such that $SW_{-}=WH$, where $S=\mathrm{diag}(s_{1},\ldots,s_{M})$ and $W_{-}$ is $W$ without the $(n+1)$-th column.

Now, we would like to employ the same idea for the Hermite interpolation and fitting. We start with the confluent Vandermonde system

$$\sum_{k=0}^{n}x_{j}^{k}c_{k}=f(x_{j}),\quad\sum_{k=0}^{n}kx_{j}^{k-1}c_{k}=f^{\prime}(x_{j}),\quad j=1,\ldots,m$$

or in matrix form (with $V$ as before)

$$\begin{pmatrix}V\\ V_{{}^{\prime}}\end{pmatrix}\mathbf{c}=\begin{pmatrix}\mathbf{f_{\phantom{{}^{\prime}}}}\\ \mathbf{f}_{{}^{\prime}}\end{pmatrix},\quad V_{{}^{\prime}}=\begin{pmatrix}0&1&2x_{1}&3x_{1}^{2}&\ldots&nx_{1}^{n-1}\\ 0&1&2x_{2}&3x_{2}^{2}&\ldots&nx_{2}^{n-1}\\ \vdots&\vdots&\vdots&\cdots&\vdots\\ 0&1&2x_{m}&3x_{m}^{2}&\ldots&nx_{m}^{n-1}\end{pmatrix},\quad\mathbf{f}_{{}^{\prime}}=\begin{pmatrix}f^{\prime}(x_{1})\\ f^{\prime}(x_{2})\\ \vdots\\ f^{\prime}(x_{m})\end{pmatrix}.$$

We notice that the column space of $[V;V_{{}^{\prime}}]$ (semicolon means as in MATLAB) is a Krylov subspace:

$$\begin{pmatrix}V\\ V_{{}^{\prime}}\end{pmatrix}=\left(\begin{pmatrix}\mathbf{e}\\ \mathbf{0}\end{pmatrix},\;\;\begin{pmatrix}X&O\\ I&X\end{pmatrix}\begin{pmatrix}\mathbf{e}\\ \mathbf{0}\end{pmatrix},\;\;\begin{pmatrix}X&O\\ I&X\end{pmatrix}^{2}\begin{pmatrix}\mathbf{e}\\ \mathbf{0}\end{pmatrix},\;\;\ldots,\begin{pmatrix}X&O\\ I&X\end{pmatrix}^{n}\begin{pmatrix}\mathbf{e}\\ \mathbf{0}\end{pmatrix}\right).$$

So, the Arnoldi process applies and can be used to generate a matrix $Q$ with orthogonal columns and a Hessenberg matrix $H$. Then, the evaluation matrix $W$ can be constructed with

$$\text{the first column }\begin{pmatrix}\mathbf{e}\\ \mathbf{0}\end{pmatrix},\quad\begin{pmatrix}S&O\\ I&S\end{pmatrix}W_{-}=WH.$$

Following [1], we give the MATLAB codes for the method of confluent Vandermonde with Arnoldi.

     function [d,H] = polyfitAh(x,f,fp,n)
     m = size(x,1); Q = zeros(2*m,n+1); Q(1:m,1) = 1; H = zeros(n+1,n);
     for k = 1:n
         q = [x.*Q(1:m,k); Q(1:m,k)+x.*Q(m+1:2*m,k)];
         for j = 1:k
             H(j,k) = Q(:,j)’*q/m; q = q - H(j,k)*Q(:,j);
         end
         H(k+1,k) = norm(q)/sqrt(m); Q(:,k+1) = q/H(k+1,k);
     end
     d = Q\[f;fp];

     function [y,yp] = polyvalAh(d, H, s)
     M = size(s,1);  n = size(d,1)-1; W = zeros(2*M,n+1); W(1:M,1) = 1;
     for k = 1:n
         w = [s.*W(1:M,k); W(1:M,k)+s.*W(M+1:2*M,k)];
         for j = 1:k
             w = w - H(j,k)*W(:,j);
         end
         W(:,k+1) = w/H(k+1,k);
     end
     y = W*d;  yp = y(M+1:2*M,:);  y = y(1:M,:);

To compare, we give also the naive method of confluent Vandermonde without Arnoldi.

     function c = polyfith(x,f,fp,n)
     A = x.^(0:n); o = zeros(size(x)); D = diag(1:n); A = [A; o A(:,1:n)*D]; c = A\[f;fp];

     function [y,yp] = polyvalh(c,s)
     n = length(c)-1;  B = s.^(0:n);  y = B*c;
     o = zeros(size(s));  D = diag(1:n);  yp = [o B(:,1:n)*D]*c;

Other than the Hermite scenario for both fitting and evaluation, we can also fit only the function values using polyfitA (or polyfit) from [1] but evaluates the function values and derivatives using polyvalAh (or polyvalh) here. This scenario is more interesting if the data of derivatives are not available for fitting. In other words, we construct first a discrete orthogonal basis for the function values by polyfitA and then the derivative matrix W(M+1:2*M,:) under the same basis by polyvalAh. We notice that an equivalent formula for doing this has been given in [2], which we should detail in the Discussion section. There is a third possibility: one can still construct a discrete orthogonal basis for the Hermite data, but fit only the function values, and then construct the derivative evaluation matrix under that basis; however we found no advantage of it. The fourth possibility is to use the function values orthogonalized basis for Hermite fitting and evaluation, which we found performs equally well as using the Hermite data orthogonalized basis for Examples 1–2 in the following sections. We will see also for the other examples the comparison results of Hermite and function values orthogonalized bases.

Both the fitting and evaluation can be generalized to high-order derivatives. For example, let $V_{{}^{\prime\prime}}$ be the second-order derivative matrix under the standard basis. We have

$$\begin{pmatrix}V\\ V_{{}^{\prime}}\\ V_{{}^{\prime\prime}}\end{pmatrix}=\left(\begin{pmatrix}\mathbf{e}\\ \mathbf{0}\\ \mathbf{0}\end{pmatrix},\;\;\begin{pmatrix}X&O&O\\ I&X&O\\ O&2I&X\end{pmatrix}\begin{pmatrix}\mathbf{e}\\ \mathbf{0}\\ \mathbf{0}\end{pmatrix},\;\;\begin{pmatrix}X&O&O\\ I&X&O\\ O&2I&X\end{pmatrix}^{2}\begin{pmatrix}\mathbf{e}\\ \mathbf{0}\\ \mathbf{0}\end{pmatrix},\;\;\ldots,\begin{pmatrix}X&O&O\\ I&X&O\\ O&2I&X\end{pmatrix}^{n}\begin{pmatrix}\mathbf{e}\\ \mathbf{0}\\ \mathbf{0}\end{pmatrix}\right).$$

It may be worthwhile to state the general case as follows.

This example is used in [1] for the Lagrange interpolation. Now we interpolate the function $f(x)=1/(1+25x^{2})$ and its derivative simultaneously (polyfitAh+polyvalAh) by a degree $n$ (odd) polynomial in $m=(n+1)/2$ Chebyshev points $x_{j}=\cos((m-j)\pi/(m-1))$, $1\leq j\leq m$. We also do the Lagrange interpolation in $n+1$ points as in [1] and evaluate the derivative as proposed here (polyfitA+polyvalAh). The results are shown in Figure 1. We see that with Arnoldi both $f$ and $f^{\prime}$ can be approximated to high accuracy with the exponential convergence rate, while without Arnoldi the method is prone to round-off errors as $n>50$. The derivative evaluation enabled by the confluent Vandermonde with Arnoldi (polyvalAh) works equally well for the Hermite and Lagrange interpolation. Throughout the paper, the function norm is the maximum norm computed on a very fine grid.

For the function $f(x)=\sqrt{|x|}$ on $[-1,-1/3]\cup[1/5,1]$, we sample the function values (and its derivatives for the Hermite version) on every interval with equally spaced points, and fit by least-squares these data with a degree $n$ polynomials. The number of sampling points $m$ on each interval is chosen as $m=5(n+1)$ for the Hermite version and $m=10(n+1)$ for the Lagrange version. The results are shown in Figure 2.

This example is used in [1] for approximation of $f(x)=1/(10-9x)$ on $[-1,1]$ by the Fourier series on the larger interval $[-2,2]$. Now we consider both the function values and its derivatives:

$$f(x)\approx\mathrm{Re}\left(\sum_{k=0}^{n}\mathrm{e}^{\mathrm{i}k\pi x/2}c_{k}\right),\quad-\frac{2}{\pi}f^{\prime}(x)\approx\mathrm{Im}\left(\sum_{k=0}^{n}k\mathrm{e}^{\mathrm{i}k\pi x/2}c_{k}\right).$$

Following [1], this is realized with $z:=\mathrm{e}^{\mathrm{i}\pi x/2}$, $c_{k}=a_{k}+\mathrm{i}b_{k}$ ($a_{k},b_{k}\in\mathbb{R}$) and

$$f(x)\approx\sum_{k=0}^{n}(\mathrm{Re}z^{k}a_{k}-\mathrm{Im}z^{k}b_{k}),\quad-\frac{2}{\pi}f^{\prime}(x)\approx\sum_{k=0}^{n}\left(\mathrm{Im}(kz^{k})a_{k}+\mathrm{Re}(kz^{k})b_{k}\right).$$

The last line of polyfith is changed to

     m = length(x);  c = [real(A(1:m,:))        imag(A(1:m,2:n+1)); ...
                          imag(A(m+1:2*m,:))   -real(A(m+1:2*m,2:n+1))]\[f;fp];
     c = c(1:n+1) - 1i*[0; c(n+2:2*n+1)];

where fp stores the scaled derivatives $-\frac{2}{\pi}f^{\prime}(x_{j})$. Similar changes should be made also in polyfitAh but with $Q$ instead of $A$. In polyvalh and polyvalAh, now we take the real part for y and imaginary part for yp. Note again that the output yp stores the scaled derivatives $-\frac{2}{\pi}p^{\prime}(s_{j})$. Following [1], 500 Chebyshev points on $[-1,1]$ are used for the Hermite fitting. The results are shown in Figure 3. The Lagrange fitting is not present in the figure since the derivative evaluation for it gave totally wrong results.

Given a domain $\Omega\subset\mathbb{R}^{2}$ with boundary $\partial\Omega$ and a function $f:\partial\Omega\to\mathbb{R}$, the Dirichlet-to-Neumann map $T$ involves first solving the Laplace equation $-\Delta u=0$ in $\Omega$, $u=f$ on $\partial\Omega$ and then evaluating the outer normal derivative $\partial_{\boldsymbol{\nu}}u=:Tf$ on $\partial\Omega$. By identifying $(x,y)\in\mathbb{R}^{2}$ and $z=x+\mathrm{i}y\in\mathbb{C}$, $u$ can be viewed as the real part of an analytic function $h$ in the complex plane, so that $-\Delta u=0$ is automatically satisfied. Following [1], we approximate $h(z)$ by a polynomial $\sum_{k=0}^{n}c_{k}z^{k}$. Let $c_{k}=a_{k}+\mathrm{i}b_{k}$ ($a_{k}$, $b_{k}$ $\in\mathbb{R}$). Then, $u(z)\approx\sum_{k=0}^{n}(\mathrm{Re}z^{k}a_{k}-\mathrm{Im}z^{k}b_{k})$. To compute $Tf$, we need only to fit the boundary condition $u=f$ on $\partial\Omega$ and then evaluate $\partial_{\boldsymbol{\nu}}u=\mathrm{Re}(\nu(z)h^{\prime}(z))$ (since $h^{\prime}(z)=\partial_{x}u-\mathrm{i}\partial_{y}u$ and $\nu(z)=\nu_{1}(z)+\mathrm{i}\nu_{2}(z)$ corresponds to the unit outer normal $\boldsymbol{\nu}=(\nu_{1},\nu_{2})$ on $\partial\Omega$), or more precisely

$$\partial_{\boldsymbol{\nu}}u=\sum_{k=1}^{n}(\nu_{1}(z)\mathrm{Re}z^{k-1}k-\nu_{2}(z)\mathrm{Im}z^{k-1}k)a_{k}-\sum_{k=1}^{n}(\nu_{2}(z)\mathrm{Re}z^{k-1}k+\nu_{1}(z)\mathrm{Im}z^{k-1}k)b_{k}.$$

A concrete example is demonstrated in Figure 4, where $m=10n$ points with the argument of $z\in\partial\Omega$ equally spaced are used for the fitting. We see that the Hermite orthogonalized basis leads to a bit more accurate results than the function values orthogonalized basis does.

This is a continuation of the previous section. We look for the eigenvalue $\lambda$ of $T$ and the assoicated eigenfunction $f:\partial\Omega\to\mathbb{R}$ such that $Tf=\lambda f$. In other words, we want the real part $u$ of an analytic function in $\Omega$ to satisfy $\partial_{\boldsymbol{\nu}}u=\lambda u$ on $\partial\Omega$. Let $Q$ be the matrix from the Arnoldi process on the confluent Vandermonde matrix; see the Introduction section. We partition the $2m$ rows of $Q$ in half: $Q=(Q_{0};Q_{1})$ with $Q_{0}$ containing the first $m$ rows that represents the sampling of basis function values and $Q_{1}$ representing the sampling of derivatives of basis functions. Let the matrices be separated into real parts and imaginary parts: $Q_{0}=Q_{0,r}+\mathrm{i}Q_{0,i}$, $Q_{1}=Q_{1,r}+\mathrm{i}Q_{1,i}$, and $\nu_{1}$, $\nu_{2}$ be diagonal matrices containing the sampled values of the unit outer normal of $\partial\Omega$. Then, the discrete Steklov eigenvalue problem is to find $\mathbf{a}\in\mathbb{R}^{n+1}$, $\mathbf{b}=(0,b_{1},\ldots,b_{n})^{T}$ and $\lambda\in\mathbb{R}$ such that

$$(\nu_{1}Q_{1,r}-\nu_{2}Q_{1,i})\mathbf{a}-(\nu_{1}Q_{1,i}+\nu_{2}Q_{1,r})\mathbf{b}=\lambda(Q_{0,r}\mathbf{a}-Q_{0,i}\mathbf{b}),$$

denoted by $Q_{\nu}\boldsymbol{\beta}=\lambda Q_{00}\boldsymbol{\beta}$ with $\boldsymbol{\beta}=(\mathbf{a};-\hat{\mathbf{b}})$ ($\hat{\mathbf{b}}=(b_{1},\ldots,b_{n})^{T}$), which is a rectangular eigenvalue problem when $m>2n+1$. Following [6], we first do the economic QR decomposition $Q_{00}=\tilde{Q}_{00}\tilde{R}_{00}$, then multiply the eigenvalue equation with $\tilde{Q}_{00}^{*}$ to get $\tilde{Q}_{00}^{*}Q_{\nu}\boldsymbol{\beta}=\lambda\tilde{R}_{00}\boldsymbol{\beta}$, which is a generalized eigenvalue problem of square matrices and can be solved by the standard method (using the MATLAB command eig). For the solved example, see Figure 5. The $m=10n+1$ points are sampled at equally spaced elliptic angle coordinates. The benchmark solution uses $n=400$.

This example is taken from [8]. The domain is $\Omega=(0,1)^{2}$. On the left, right and bottom sides, the homogeneous Neumann condition $\partial_{\boldsymbol{\nu}}u=0$ is imposed, while on the top side is the eigenvalue equation $\partial_{\boldsymbol{\nu}}u=\lambda u$. The analytic solution is $\lambda_{n}=n\pi\tanh(n\pi)$, $u_{n}=\cos(n\pi x)\sinh(n\pi y)$; see [8]. Let $I:u\mapsto u$ be the identity operator on the top side and $O:u\mapsto 0u$ be the zero operator on the other sides. Define $B$ on $\partial\Omega$ piecewise equal to $I$ and $O$. Then, we can write the boundary equation uniformly as $\partial_{\boldsymbol{\nu}}u=\lambda Bu$ on $\partial\Omega$. Using the notation from the previous section, the discrete rectangular eigenvalue problem becomes $Q_{\nu}\boldsymbol{\beta}=\lambda BQ_{00}\boldsymbol{\beta}$. Now, testing the equation with the column space basis of $BQ_{00}$ as before is not a good idea because $BQ_{00}$ is rank deficient. Note that the two sides of the equation live in the sum of the column spaces of $Q_{\nu}$ and $BQ_{00}$. So, we follow [5] to do the economic SVD: $(Q_{\nu},BQ_{00})=USV^{*}$, but we keep the first $2n+1$ columns of $U$ as $\hat{U}$ to match the length of $\boldsymbol{\beta}$, then multiply the eigenvalue equation with $\hat{U}^{*}$ to get $\hat{U}^{*}Q_{\nu}\boldsymbol{\beta}=\hat{U}^{*}BQ_{00}\boldsymbol{\beta}$. The results are shown in Figure 6. On each side of the square, we used $20(n+1)$ Chebyshev points of the first kind for discretization of the boundary equation. From the bottom right, we can see that the function values orthogonalized basis leads to instability. Convergence of the eigenfunctions from the method with Hermite Arnoldi basis slowers down after some $n$, which could be fixed with more sampling points. For example, for $n=60$ and $100(n+1)$ points per side, we got the errors $2.4\times 10^{-14}$, $1.2\times 10^{-13}$ for $p_{5}$, $p_{10}$.

We want to find a polynomial approximation of $\int f(x)\mathrm{d}x\approx p(x)+C$, $C$ is a general constant. For example, if $g(x)=\int_{a}^{x}u(t)~{}\mathrm{d}t$, we have $g^{\prime}(x)=f(x)$, $g(a)=0$, the simplest differential equation with a point value condition. The goal is to find $p(x)\approx g(x)$. The implementation is mostly the same as polyfitAh but the function header and the last line: function [d, H] = polyfitAih(x,f,n) and d = Q(m+1:2*m,2:n+1)\f; d = [0; d]; alternatively, one can orthogonalize only the function values and from the resulting Hessenberg matrix generate the derivative matrix. The solved example is shown in Figure 7. From the bottom right, we can see that, for both the Hermite Arnoldi basis and the function values Arnoldi basis, the error of $p$ can grow with $n$ in some range but still decreases to $10^{-15}$ afterwards. Such transient instability does not appear for $f(x)=1/(1+25x^{2})$ on the Chebyshev grids of $[-1,1]$.

This work is greatly inspired by the fascinating paper Vandermonde with Arnoldi [1] by Brubeck, Nakatsukasa and Trefethen. We saw a door open to the well-conditioned bases of polynomial spaces on arbitrary grids, and to many interesting applications, e.g., the lightning solvers [2, 3, 9, 10]. In the spirit of Chebfun [4], it is natural then to consider differential, integral and other operations of polynomials under the Arnoldi basis. Under the standard basis, the differential operation or Hermite interpolation leads to the confluent Vandermonde matrix, which is ill-conditioned [7]. The Arnoldi process can come to the rescue. However, it is not so straightforward to apply the method from [1]. The key idea behind the new method is the discovery of a connection of the confluent Vandermonde matrix to a Krylov subspace. This enables us to represent and compute the derivatives in a stable way. Two approaches are proposed to accomplish this task. One is constructing the basis for the column space of the confluent Vandermonde matrix directly by the Confluent Vandermonde with Arnoldi method based on the confluent Krylov subspace that we discovered. Another approach is using the basis generated by Vandermonde with Arnoldi and constructing the derivative matrix of the basis by exploiting the confluent Krylov subspace. With the new functionality of computing derivatives, we are able to do the Hermite interpolation/evaluation to high precision, to accomodate the Neumann/Robin boundary conditions, to compute the Dirichlet-to-Neumann map and its eigenmodes, and to do the indefinite integral.

We realized just upon finalizing this note that a formula of the derivative matrix has been given in [9, 2]. It is based on the recursive relation of the polynomial basis from Vandermonde with Arnoldi. Recall that the Arnoldi process produces the basis matrix $Q$ and the Hessenberg matrix $H$ that satisfy $XQ_{-}=QH$. If we imagine that the number of rows of $X$ and $Q$ goes to infinity, then we would get the polynomial recursion $xp_{k}(x)=\sum_{j=1}^{k}p_{j}(x)h_{jk}+p_{k+1}(x)h_{k+1,k}$ or $p_{k+1}(x)=(xp_{k}(x)-\sum_{j=1}^{k}p_{j}(x)h_{jk})/h_{k+1,k}$. Taking derivative gives $p_{k+1}^{\prime}(x)=(p_{k}(x)+xp_{k}^{\prime}(x)-\sum_{j=1}^{k}p_{j}^{\prime}(x)h_{jk})/h_{k+1,k}$. It can be seen that our second approach (e.g. polyfitA + polyvalAh) is equivalent to this formula.

Comparison of the Hermite orthogonalized basis and the function values orthogonalized basis has been made on each example. It is found that they perform equally well for the Hermite interpolation/evaluation task. But the function values orthogonalized basis is unstable for the example of sloshing modes, while the Hermite orthogonalized basis works stably. We also observed instability of both the approaches for the Fourier extension example and the indefinite integral example. A theoretical understanding of the comparison results and the instability phenomena could be interesting.