On the Backward Error Incurred by the Compact Rational Krylov Linearization

Abstract

One of the most successful methods for solving a polynomial (PEP) or rational eigenvalue problem (REP) is to recast it, by linearization, as an equivalent but larger generalized eigenvalue problem which can be solved by standard eigensolvers. In this work, we investigate the backward errors of the computed eigenpairs incurred by the application of the well-received compact rational Krylov (CORK) linearization. Our treatment is unified for the PEPs or REPs expressed in various commonly used bases, including Taylor, Newton, Lagrange, orthogonal, and rational basis functions. We construct one-sided factorizations that relate the eigenpairs of the CORK linearization and those of the PEPs or REPs. With these factorizations, we establish upper bounds for the backward error of an approximate eigenpair of the PEPs or REPs relative to the backward error of the corresponding eigenpair of the CORK linearization. These bounds suggest a scaling strategy to improve the accuracy of the computed eigenpairs. We show, by numerical experiments, that the actual backward errors can be successfully reduced by scaling and the errors, before and after scaling, are both well predicted by the bounds.

A Proof of (14)

We now show that $G(\lambda )$ given by (14) satisfies (12) with $\mathbf{L} _m(\lambda )$ corresponding to the orthogonal basis and $g(\lambda ) = 1$. To simplify the notation, we drop the argument $\lambda$ of the basis function $b_j^{(k)}(\lambda )$ in the rest of this proof.

We shall concentrate on showing that the product of $G(\lambda )$ and the kth column of $\mathbf{L} _m(\lambda )$ is kth column of $e_1^T\otimes R_m(\lambda )$, that is

$$\begin{aligned} D_{k-1} - \sum _{j=k-1}^m b_{j-(k-1)}^{(k-1)} D_j + \frac{\lambda + \alpha _{k-1}}{\beta _k}\sum _{j=k}^m b_{j-k}^{(k)} D_j + \frac{\gamma _k}{\beta _{k+1}}\sum _{j = k+1}^m b_{j-(k+1)}^{(k+1)} D_j = 0, \end{aligned}$$

(28)

which is true if

$$\begin{aligned}&b_{0}^{(k-1)} = 1, \end{aligned}$$

(29a)

$$\begin{aligned}&\frac{\lambda +\alpha _{k-1}}{\beta _k} b_{j-k}^{(k)} + \frac{\gamma _k}{\beta _{k+1}} b_{j-(k+1)}^{(k+1)} = b_{j-(k-1)}^{(k-1)}, \quad k \leqslant j \leqslant m, \end{aligned}$$

(29b)

both hold.

Equation (29a) is given by the definition of the orthogonal basis. We now show (29b) by induction.

For $j = k$, (29b) reads

$$\begin{aligned} \frac{\lambda +\alpha _{k-1}}{\beta _k} b_0^{(k)} + \frac{\gamma _k}{\beta _{k+1}} b_{-1}^{(k+1)} = b_1^{(k-1)}, \end{aligned}$$

which reduces to

$$\begin{aligned} \frac{\lambda +\alpha _{k-1}}{\beta _k} b_0^{(k-1)} = b_1^{(k-1)}, \end{aligned}$$

(30)

since $b_{-1}^{(k+1)} = 0$ and $b_0^{(k)} = b_0^{(k-1)} = 1$. Equation (30) is the recurrence relation (13) with 0 and $k-1$ in place of j and k, respectively.

For $j = k+1$, (29b) becomes

$$\begin{aligned} \frac{\lambda +\alpha _{k-1}}{\beta _k} b_1^{(k)} + \frac{\gamma _k}{\beta _{k+1}} b_0^{(k+1)} = b_2^{(k-1)}. \end{aligned}$$

(31)

To show (31), we substitute $b_1^{(k)} = \frac{\lambda + \alpha _k}{\beta _{k+1}}b_0^{(k)}$ into to have

$$\begin{aligned} \frac{\lambda + \alpha _k}{\beta _{k+1}}\left( \frac{\lambda + \alpha _{k-1}}{\beta _k} b_0^{(k-1)}\right) + \frac{\gamma _k}{\beta _{k+1}}b_0^{(k-1)} = b_2^{(k-1)}, \end{aligned}$$

where $b_0^{(k-1)} = b_0^{(k)} = b_0^{(k+1)} = 1$ is used. Since the terms in the parentheses, by the recurrence relation, are just $b_1^{(k-1)}$, the task is boiled down to verify

$$\begin{aligned} \frac{\lambda + \alpha _k}{\beta _{k+1}}b_1^{(k-1)} + \frac{\gamma _k}{\beta _{k+1}}b_0^{(k-1)} = b_2^{(k-1)}. \end{aligned}$$

This is, again, a variant of the recurrence relation of (13).

Now suppose that (29b) and

$$\begin{aligned} \frac{\lambda +\alpha _{k-1}}{\beta _k} b_{j+1-k}^{(k)} + \frac{\gamma _k}{\beta _{k+1}} b_{j+1-(k+1)}^{(k+1)} = b_{j+1-(k-1)}^{(k-1)} \end{aligned}$$

(32)

both hold for any $k-1 \leqslant j \leqslant m-2$. By the recurrence relation (13), we have

$$\begin{aligned} (\lambda +\alpha _{j+1}) b^{(k-1)}_{j-k+2} + \gamma _{j+1} b^{(k-1)}_{j-k+1} = \beta _{j+2}b^{(k-1)}_{j-k+3}, \end{aligned}$$

where we now substitute (29b) and (32) for $b_{j-(k-1)}^{(k-1)}$ and $b_{j+1-(k-1)}^{(k-1)}$ respectively to have

$$\begin{aligned} \left( \frac{\lambda +\alpha _{k-1}}{\beta _k} b_{j-k+2}^{(k)} + \frac{\gamma _k}{\beta _{k+1}} b_{j-k+1}^{(k+1)}\right) \beta _{j+2} = b_{j-k+3}^{(k-1)}\beta _{j+2}, \end{aligned}$$

(33)

where we have used two variants of the recurrence relation

$$\begin{aligned} (\lambda +\alpha _{j+1}) b^{(k)}_{j-k+1} + \gamma _{j+1} b^{(k)}_{j-k}&= \beta _{j+2}b^{(k)}_{j-k+2}, \\ (\lambda +\alpha _{j+1}) b^{(k+1)}_{j-k} + \gamma _{j+1} b^{(k+1)}_{j-k-1}&= \beta _{j+2}b^{(k+1)}_{j-k+1}. \end{aligned}$$

Dropping $\beta _{j+2}$’s in (33) yields

$$\begin{aligned} \frac{\lambda +\alpha _{k-1}}{\beta _k} b_{j+2-k}^{(k)} + \frac{\gamma _k}{\beta _{k+1}} b_{j+2-(k+1)}^{(k+1)} = b_{j+2-(k-1)}^{(k-1)}. \end{aligned}$$

This, by induction, shows (29b), which, along with (29a), verifies (28).

The proofs for the first and the last two columns of $G(\lambda )\mathbf{L} _m(\lambda )$ are essentially the same, though the structures of $\mathbf{L} _m(\lambda )$ for these columns are slightly different from that of a general kth column.