The roots of a monic polynomial expressed in a Chebyshev basis are known to be
the eigenvalues of the so-called colleague matrix, which is a Hessenberg
matrix that is the sum of a symmetric tridiagonal matrix and a rank-1 matrix.
The rootfinding problem is thus reformulated as an eigenproblem, making the
computation of the eigenvalues of such matrices a subject of significant
practical importance. In this manuscript, we describe an $O(n^{2})$ explicit
structured QR algorithm for colleague matrices and prove that it is
componentwise backward stable, in the sense that the backward error in the
colleague matrix can be represented as relative perturbations to its
components. A recent result of Noferini, Robol, and Vandebril shows that
componentwise backward stability implies that the backward error $\delta c$ in
the vector $c$ of Chebyshev expansion coefficients of the polynomial has the
bound ${\lVert\delta c\rVert}\lesssim{\lVert c\rVert}\mathit{u}$, where
$\mathit{u}$ is machine precision. Thus, the algorithm we describe has both
the optimal backward error in the coefficients and the optimal cost
$O(n^{2})$. We illustrate the performance of the algorithm with several
numerical examples.

A Provably Componentwise Backward Stable $\mathbf{O(n^{2})}$ QR Algorithm for
the Diagonalization of Colleague Matrices

K. Serkh$\mbox{}^{\dagger\,\diamond}$, V. Rokhlin$\mbox{}^{\ddagger\,\oplus}$,

$\mbox{}^{\diamond}$This author’s work was supported in part by the NSERC
Discovery Grants RGPIN-2020-06022 and DGECR-2020-00356.

$\mbox{}^{\oplus}$ This author’s work was supported in part under ONR
N00014-18-1-2353 and NSF DMS-1952751.

$\mbox{}^{\dagger}$ Dept. of Math. and Computer Science, University of
Toronto, Toronto, ON M5S 2E4

$\mbox{}^{\ddagger}$ Dept. of Mathematics, Yale University, New Haven, CT
06511

###### Contents

  1. 1 Introduction
    1. 2 Preliminaries
      1. 2.1 Linear Algebra
      2. 2.2 Error Analysis
        1. 2.2.1 Floating point computation of complex plane rotations
        2. 2.2.2 Multiplication by complex plane rotations
      3. 2.3 Colleague Matrices and Polynomial Rootfinding
      4. 2.4 Stability of Rootfinding Using Linearizations
      5. 2.5 Conventions
        1. 3 The Algorithm
          1. 3.1 Eliminating the Superdiagonal
          2. 3.2 Rotating Back to Hessenberg Form
          3. 3.3 The QR Algorithms
            1. 4 Componentwise Backward Stability
              1. 4.1 Forward Error Analysis of a Single Sweep of $QR$
  2. 4.2 Backward Error Analysis of the $QR$ Algorithms
  3. 5 Numerical Results
    1. 5.1 $p_{\text{rand}}(x)$: Polynomials with Random Coefficients
    2. 5.2 $p_{\text{wilk}}(x)$: Wilkinson’s Polynomial
    3. 5.3 $f_{\text{sin}}(x)$: A Smooth Function
    4. 5.4 $p_{\text{mult}}(x)$: A Polynomial with Multiple Roots
    5. 5.5 $p_{\text{yuji}}(x)$: A Pathological Example from [26]
    6. 5.6 $f_{\text{cas}}(x)$: A Pathological Example from [14]
    7. 5.7 CPU Times
      1. 6 Conclusions and Generalizations

##  1 Introduction

The problem of finding the roots of the polynomial

$\displaystyle p(x)=c_{0}+c_{1}x+\cdots+c_{n-1}x^{n-1}+x^{n}$ (1)

is one of the oldest and most classical problems in mathematics. Countless
methods have been proposed for its solution (see, for example, [28] for a
history, and the two volumes [24] and [25] for a detailed account of such
methods). In the 1800’s, it was observed by Frobenius that the roots of the
polynomial are the eigenvalues of a certain matrix called the companion
matrix, formed using the polynomial coefficients. A matrix whose eigenvalues
are the roots of $p(x)$ is called a linearization of $p(x)$. Given a
linearization of $p(x)$, the roots of the polynomial can thus be recovered by
computing the eigenvalues of the matrix.

If the roots of the polynomial are found numerically, than the computed roots
can be viewed as the exact roots of a perturbed polynomial $p(x)+\delta p(x)$
with coefficients $c_{i}+\delta c_{i}$, where the size of the vector $\delta
c$ is called the backward error in the polynomial coefficients. The best
possible bound on the backward error that a linearization method can have for
general polynomials is

$\displaystyle{\lVert\delta c\rVert}\lesssim{\lVert c\rVert}\mathit{u},$ (2)

or, in other words, that the relative normwise backward error in the
polynomial coefficients is bounded by machine precision $\mathit{u}$ (see, for
example, [26]). The backward error in the companion matrix method was revealed
by the influential paper [19], which analyzed the relationship between
perturbations in the companion matrix and perturbations in the polynomial
coefficients. There, the authors proved that if the companion matrix $C$ is
perturbed by the matrix $E$, then the matrix $C+E$ is a linearization of the
polynomial with perturbed coefficients $c_{i}+\delta c_{i}$, and that the
perturbation satisfies the normwise bound

$\displaystyle{\lVert\delta c\rVert}\lesssim{\lVert c\rVert}{\lVert
E\rVert}\mathit{u}.$ (3)

If the eigenvalues are computed by a standard QR algorithm, which is known to
be backward stable (see, for example, [35]), then the computed eigenvalues are
the exact eigenvalues of $C+E$, where ${\lVert E\rVert}\lesssim{\lVert
C\rVert}\mathit{u}$. Since ${\lVert C\rVert}\approx{\lVert c\rVert}$, it
follows that the backward error in the polynomial coefficients is bounded by

$\displaystyle{\lVert\delta c\rVert}\lesssim{\lVert c\rVert}^{2}\mathit{u}.$
(4)

Thus, as ${\lVert c\rVert}$ get larger, the relative backward error in the
coefficients increases. The companion matrix method, at first glance, would
appear then to have two drawbacks: it falls short of the optimal backward
error bound (2), and it costs $O(n^{3})$ operations as a result of using the
QR algorthim.

The situation improved dramatically in 2007, when Bini, Eidelman, Gemignani,
and Gohberg published a paper [10] describing a stable, $O(n^{2})$ explicit QR
method for companion matrices (around the same time, Chandrasekaran, Gu, Xia,
and Zhu also discovered an $O(n^{2})$ method for companion matrices, see
[15]). The algorithm is based on the observation that the companion matrix and
its QR iterates have a certain structure which allows them to be represented
by a collection of $O(n)$ parameters called generators (specifically, the
companion matrix is a Hessenberg matrix that is the sum of a unitary matrix
and a rank-1 perturbation; matrices of this form are called fellow matrices).
In 2010, an implicit version of this algorithm, also stable and $O(n^{2})$,
and also based on generators, was introduced in [9]. Around the same time, Van
Barel, Vandebril, Van Dooren, and Frederix discovered in [8] an alternative
stable, $O(n^{2})$ implicit QR algorithm based on representing the unitary
part by so-called core transformations, which are rotation matrices acting
only on two adjacent rows at a time (see, for example, [2]). The first example
of a proof of backward stability for an implicit $O(n^{2})$ QR algorithm for
companion matrices was given by Aurentz, Mach, Vandebril, and Watkins in [4];
this algorithm is again based on core transformations. The backward stability
result accompanying this QR algorithm guarantees the sub-optimal bound (2),
but has the optimal complexity of $O(n^{2})$. Amazingly, the authors then
discovered that the algorithm they had constructed, with some minor
modifications, actually yields the optimal bound (4) in practice. An
investigation showed that the reason for this remarkable behavior is that
their algorithm is not just backward stable, but is componentwise backward
stable, meaning that the backward error in the companion matrix can be
decomposed into proportional backward errors in each of its components. They
published a proof of the componentwise backward stability of their algorithm,
together with a proof that componentwise backward stability guarantees the
bound (2), in [3], along with numerical experiments.

Thus, if the coefficients of the polynomial $p(x)$ in the monomial basis are
known, then the algorithm of [3] is optimal in both error and time complexity.
However, if the coefficients are not known, then the companion matrix cannot
be used to the find the roots accurately, since the relationship between the
values of the polynomial $p(x)$ and the coefficients of its monomial expansion
is highly unstable (this fact has been known for many decades, at least as
early as Wilkinson [31]). If the polynomial $p(x)$ is instead expanded in a
basis of Chebyshev polynomials

$\displaystyle p(x)=c_{0}+c_{1}T_{1}(x)+\cdots+c_{n-1}T_{n-1}(x)+T_{n}(x),$
(5)

where $T_{i}(x)$ is the Chebyshev polynomial of order $i$, then the
relationship between the coefficients and the polynomial is perfectly stable
(see, for example, [36]). In fact, this observation is the basis for the
Chebfun software package (see [6] and [16]). An analogue of the companion
matrix, constructed from the Chebyshev expansion coefficients, was discovered
in 1961 by Good [21], who called it the colleague matrix, and independently by
Spect in 1957 [33]–[34]. The first $O(n^{2})$ algorithm for colleague matrices
was discovered by Bini, Gemignani, and Pan in 2005 (even before [10] appeared)
and is a stable, explicit QR algorithm based on generators [11]. Like the
companion matrix, the colleague matrix has a special structure that is
preserved over QR iterations (specifically, the colleague matrix is a
Hessenberg matrix that is the sum of a Hermitian matrix and a rank-1
perturbation). In 2008, Eidelman, Gemignani, and Gohberg, in [20], introduced
a stable, $O(n^{2})$ implicit QR algorithm. The relationship between the
backward error in the Chebyshev expansion coefficients and perturbations to
the colleague matrix was first investigated Nakasukasa and Noferini in [26],
where the authors found a lower bound for the backward error in the
coefficients, showing that a backward stable QR algorithm can do no better
than (4) (around the same time, Lawrence, Van Barel, and Van Dooren published
a general analysis in [23], where they also proved a lower bound for colleague
matrices). In [30], Perez and Noferini improved on this result and found an
upper bound as well, proving that if the perturbation to the colleague matrix
is small, then the bound (2) is achieved. The relationship between
componentwise perturbations to the colleague matrix and the backward error in
the coefficients was described completely in 2019 by Noferini, Robol, and
Vandebril in [27].

Recently (see [32] and [14]), it was observed that certain $O(n^{2})$
structured QR algorithms for colleague matrices are surprisingly stable,
attaining the bound (2) in many cases, an observation that mirrors the
discovery in [3] for the case of companion matrices. However, unlike in [3],
all previously proposed $O(n^{2})$ structured QR algorithms for colleague
matrices have polynomials for which the worst-case bound (4) is attained.
Thus, the question of whether or not there exists a structured $O(n^{2})$ QR
algorithm that, when used to find the roots of a colleague matrix, attains the
optimal bound (2), has remained open. In this manuscript, we answer this
question in the affirmative by presenting a new, explicit $O(n^{2})$ QR
algorithm for colleague matrices (in fact, for all Hessenberg matrices that
have a Hermitian plus rank-1 structure), and prove that our algorithm is
componentwise backward stable. Combined with the result in [27], this amounts
to a proof that the optimal bound (2) is attained for all polynomials $p(x)$.
We demonstrate that this is indeed the case with several numerical
experiments.

The structure of this manuscript is as follows. Section 2 describes the
mathematical and numerical preliminaries. Section 3 describes the algorithm,
and explains the significance of each step. In Section 4, we prove rigorously
that the algorithm is componentwise backward stable. Section 5 presents the
results of several numerical experiments. In Section 6, we discuss possible
extensions and generalizations of the algorithm.

##  2 Preliminaries

In this section, we describe the mathematical and numerical preliminaries.

###  2.1 Linear Algebra

The following lemma states that if the sum of a Hermitian matrix and a rank-1
update $pq^{*}$ is lower Hessenberg, then the matrix is determined entirely by
its diagonal and superdiagonal together with the vectors $p$ and $q$.

######  Lemma 2.1 (Eidelman, Gemignani, Gohberg [20]).

Suppose that $A\in\mathbbm{C}^{n\times n}$ is Hermitian, and let $d$ and
$\beta$ denote the diagonal and superdiagonal of $A$, respectively. Suppose
that $p,q\in\mathbbm{C}^{n}$ and that $A+pq^{*}$ is lower Hessenberg. Then

$\displaystyle a_{i,j}=\left\\{\begin{array}[]{ll}-p_{i}q_{j}^{*}&\text{if
$j>i+1$}\\\ \beta_{i}&\text{if $j=i+1$}\\\ d_{i}&\text{if $j=i$}\\\
\overline{\beta_{j}}&\text{if $j=i-1$}\\\ -q_{j}p_{i}^{*}&\text{if
$j<i-1$}\end{array}\right.$ (11)

where $a_{i,j}$ denotes the $(i,j)$-th entry of $A$.

The following lemma states that if the sum of a matrix and a rank-1 update
$pq^{*}$ is lower triangular, then the upper Hessenberg part of the matrix is
determined entirely by its diagonal and subdiagonal, together with the vectors
$p$ and $q$.

######  Lemma 2.2.

Suppose that $B\in\mathbbm{C}^{n\times n}$ and let $d$ and $\gamma$ denote the
diagonal and subdiagonal of $B$, respectively. Suppose that
$p,q\in\mathbbm{C}^{n}$ and that $B+pq^{*}$ is lower triangular. Then

$\displaystyle b_{i,j}=\left\\{\begin{array}[]{ll}-p_{i}q_{j}^{*}&\text{if
$j>i$}\\\ d_{i}&\text{if $j=i$}\\\ \gamma_{j}&\text{if
$j=i-1$}\end{array}\right.$ (15)

where $b_{i,j}$ denotes the $(i,j)$-th entry of $B$.

The following definition introduces two matrix seminorms that we will need in
our error analysis.

######  Definition 2.1.

Suppose that $A\in\mathbbm{C}^{n\times n}$ and let $a_{i,j}$ denote the
$(i,j)$-th entry of $A$. We will use the notation ${\lVert\cdot\rVert}_{H}$ to
denote the square root of the sum of squares of the entries in the upper
Hessenberg part of a matrix, so that

$\displaystyle{\lVert A\rVert}_{H}=\sqrt{\sum_{j\geq i-1}{\lvert
a_{i,j}\rvert}^{2}}.$ (16)

Likewise, we will use the notation ${\lVert\cdot\rVert}_{T}$ to denote the
square root of the sum of squares of the entries in the upper triangular part,
so that

$\displaystyle{\lVert A\rVert}_{T}=\sqrt{\sum_{j\geq i}{\lvert
a_{i,j}\rvert}^{2}}.$ (17)

The following is a straightforward lemma stating that if a certain sequence of
transformations is applied to a matrix on the right, then the upper triangular
part of the result is determined by only the upper Hessenberg part of the
original matrix.

######  Lemma 2.3.

Suppose that $B\in\mathbbm{C}^{n\times n}$, and let
$P_{2},P_{3},\ldots,P_{n}\in\mathbbm{C}^{n\times n}$ be matrices such that
$P_{k}$ only affects the $(k-1,k)$-plane of any vector it is applied to.
Define $P\in\mathbbm{C}^{n\times n}$ by the formula $P=P_{2}P_{3}\cdots
P_{n}$. Then the upper triangular part of $BP^{*}$ is determined entirely by
the upper Hessenberg part of $B$. Furthermore, if $P_{2},P_{3},\ldots,P_{n}$
are unitary, then ${\lVert BP^{*}\rVert}_{T}\leq{\lVert B\rVert}_{H}$.

###  2.2 Error Analysis

The following definition introduces the notation used in the error analysis
that appears in this manuscript. We follow the notation used in [22] and [3].

######  Definition 2.2.

Evaluation of an expression in floating point arithmetic is denoted by
$\mathop{\mathit{fl}}(\cdot)$, and we denote the unit roundoff (or machine
epsilon) by $\mathit{u}$. We assume that

$\displaystyle\mathop{\mathit{fl}}(x\mathop{\mathrm{op}}y)=(x\mathop{\mathrm{op}}y)(1+\delta),\qquad{\lvert\delta\rvert}\leq\mathit{u},$
(18)

where $\mathop{\mathrm{op}}$ stands for any of the basic arithmetic operations
$+,-,*,/$. We denote computed quantities by a hat, so that $\widehat{x}$
denotes the computed approximation to $x$. We use the notation $\lesssim$ to
mean “less than or equal to the right hand side times a modest multiplicative
constant depending on $n$ as a low-degree polynomial”, where the meaning of
$n$ is clear from the context. Whenever a matrix or vector norm appears to the
left or right of $\lesssim$, we omit the particular choice of norm, since in
finite dimensions all norms are equivalent. When $\mathit{u}$ appears in an
expression on the right hand side of $\lesssim$, we ignore all higher order
powers of $\mathit{u}$.

####  2.2.1 Floating point computation of complex plane rotations

The following lemma bounds the forward error of the floating point computation
of a complex plane rotation (see, for example, §20 of [37]).

######  Lemma 2.4.

Suppose that $x=(x_{1},x_{2})^{T}\in\mathbbm{C}^{2}$, and let
$Q\in\mathrm{SU}(2)$ be the complex rotation matrix which eliminates the first
entry, so that $(Qx)_{1}=0$. Let $\widehat{Q}\in\mathbbm{C}^{2\times 2}$ be
the floating point matrix defined by

$\displaystyle\widehat{Q}_{1,1}=c$ $\displaystyle\widehat{Q}_{1,2}=-s$ (19)
$\displaystyle\widehat{Q}_{2,1}=\overline{s}$
$\displaystyle\widehat{Q}_{2,2}=\overline{c}$ (20)

where $c=\mathop{\mathit{fl}}\bigl{(}x_{2}/\sqrt{{\lvert
x_{1}\rvert}^{2}+{\lvert x_{2}\rvert}^{2}}\bigr{)}$ and
$s=\mathop{\mathit{fl}}\bigl{(}x_{1}/\sqrt{{\lvert x_{1}\rvert}^{2}+{\lvert
x_{2}\rvert}^{2}}\bigr{)}$, and where $c=1$ and $s=0$ if ${\lVert x\rVert}=0$.
Then

$\displaystyle{\lVert\widehat{Q}-Q\rVert}\lesssim\mathit{u}.$ (22)

####  2.2.2 Multiplication by complex plane rotations

The following lemma estimates the forward error of applying a plane rotation
to a vector (see, for example, §21 of [37]).

######  Lemma 2.5.

Suppose that $Q\in\mathrm{SU}(2)$ is a complex rotation matrix, and let
$\widehat{Q}$ be a floating point approximation to $Q$ satisfying (22).
Suppose further that $x=(x_{1},x_{2})^{T}\in\mathbbm{C}^{2}$. Then

$\displaystyle{\lVert\mathop{\mathit{fl}}(\widehat{Q}x)-Qx\rVert}\lesssim{\lVert
x\rVert}u.$ (23)

###  2.3 Colleague Matrices and Polynomial Rootfinding

Suppose that $p(x)$ is a monic polynomial of order $n$ represented by

$\displaystyle p(x)=\sum_{j=0}^{n}c_{j}T_{j}(x),$ (24)

where $c_{j}\in\mathbbm{R}$, $c_{n}=1$, and $T_{j}(x)$ is the Chebyshev
polynomial of order $j$. It turns out that the roots of $p(x)$ are the
eigenvalues of the (scaled) $n\times n$ colleague matrix

$\displaystyle C=\left(\begin{array}[]{ccccc}0&\frac{1}{\sqrt{2}}&&&\\\
\frac{1}{\sqrt{2}}&0&\frac{1}{2}&\\\ &\frac{1}{2}&\ddots&\ddots&\\\
&&\ddots&0&\frac{1}{2}\\\
&&&\frac{1}{2}&0\end{array}\right)-\frac{1}{2}e_{n}\left(\begin{array}[]{ccccc}c_{0}\sqrt{2}&c_{1}&c_{2}&\cdots&c_{n-1}\end{array}\right),$
(31)

where $e_{n}$ is the $n$-th unit vector (see, for example, [21]). A matrix $C$
whose eigenvalues are the roots of $p(x)$ is called a linearization of $p(x)$.
Letting

$\displaystyle A=\left(\begin{array}[]{ccccc}0&\frac{1}{\sqrt{2}}&&&\\\
\frac{1}{\sqrt{2}}&0&\frac{1}{2}&\\\ &\frac{1}{2}&\ddots&\ddots&\\\
&&\ddots&0&\frac{1}{2}\\\ &&&\frac{1}{2}&0\end{array}\right)$ (37)

and

$\displaystyle
q^{*}=-\frac{1}{2}\left(\begin{array}[]{ccccc}c_{0}\sqrt{2}&c_{1}&c_{2}&\cdots&c_{n-1}\end{array}\right),$
(39)

we see that the colleague matrix $C$ can be written as

$\displaystyle C=A+e_{n}q^{*},$ (40)

where $A$ is Hermitian and $C$ is lower Hessenberg.

The following beautiful theorem by Noferini, Robol, and Vandebril (see
Corollary 5.4 of [27]) bounds the change in the coefficients of the polynomial
being linearized by the componentwise perturbations of the colleague matrix
$C=A+e_{n}q^{*}$.

######  Theorem 2.6.

Let $C=A+e_{n}q^{*}$ be the linearization (31) of the monic polynomial $p(x)$,
expressed in the Chebyshev basis. Consider the perturbations ${\lVert\delta
A\rVert}\leq\epsilon_{A}$, ${\lVert\delta e_{n}\rVert}\leq\epsilon_{n}$, and
${\lVert\delta q\rVert}\leq\epsilon_{q}$. Then, the matrix

$\displaystyle C+\delta C=A+\delta A+(e_{n}+\delta e_{n})(q+\delta q)^{*}$
(41)

is a linearization of the polynomial

$\displaystyle p(x)+\delta p(x)=\sum_{j=0}^{n}(c_{j}+\delta c_{j})T_{j}(x),$
(42)

where ${\lVert\delta c\rVert}\lesssim\epsilon_{n}+\epsilon_{q}+{\lVert
c\rVert}\epsilon_{A}$.

###  2.4 Stability of Rootfinding Using Linearizations

Suppose that $p(x)$ is a monic polynomial of order $n$ represented in the
Chebshev basis

$\displaystyle p(x)=\sum_{j=0}^{n}c_{j}T_{j}(x),$ (43)

where $c_{j}\in\mathbbm{R}$, $c_{n}=1$, and $T_{j}(x)$ is the Chebyshev
polynomial of order $j$, and let the roots of $p(x)$ be denoted by
$x_{1},x_{2},\ldots x_{n}\in\mathbbm{C}$. Suppose that a rootfinding algorithm
returns the computed roots
$\widehat{x}_{1},\widehat{x}_{2},\ldots,\widehat{x}_{n}\in\mathbbm{C}$. If the
computed roots are the exact roots of some perturbed polynomial

$\displaystyle p(x)+\delta p(x)=\sum_{j=0}^{n}(c_{j}+\delta c_{j})T_{j}(x),$
(44)

where

$\displaystyle\frac{{\lVert\delta c\rVert}}{{\lVert
c\rVert}}\lesssim\mathit{u},$ (45)

then we say that the rootfinding algorithm is backward stable. In fact, this
is the best backward stability bound that can be hoped for, for general
polynomials $p(x)$ (see the discussion in Appendix A of [26]).

######  Remark 2.1.

Suppose that $p(x)$ is a polynomial of order $n$, that is not monic, expressed
in the Chebyshev basis

$\displaystyle p(x)=\sum_{j=0}^{n}a_{j}T_{j}(x),$ (46)

where $a_{j}\in\mathbbm{R}$ and $T_{j}(x)$ is the Chebyshev polynomial of
order $j$. Clearly, the roots of $p(x)$ are identical to the roots of
$p(x)/a_{n}$. Let $c_{j}=a_{j}/a_{n}$, for $j=0,1,\ldots,n$. If a backward
stable rootfinding algorithm is applied to the monic polynomial $p(x)/a_{n}$,
then, letting $\delta a=\delta c\cdot a_{n}$, the algorithm is also backward
stable with respect to the original coefficients $a_{j}$, since

$\displaystyle\frac{{\lVert\delta a\rVert}}{{\lVert
a\rVert}}=\frac{\frac{1}{a_{n}}{\lVert\delta a\rVert}}{\frac{1}{a_{n}}{\lVert
a\rVert}}=\frac{{\lVert\delta c\rVert}}{{\lVert c\rVert}}\lesssim\mathit{u}.$
(47)

When linearization is used as a rootfinding algorithm, the stability of the
computed roots comes from the stability of the eigenvalue algorithm applied to
the colleague matrix $C$. If the eigenvalues of $C$ are computed by an
unstructured QR algorithm, then the backward error $\delta C$ on $C$ is
bounded by ${\lVert\delta C\rVert}\lesssim{\lVert C\rVert}\mathit{u}$. Since
the backward error is unstructured, it follows that ${\lVert\delta
A\rVert}\approx{\lVert C\rVert}\mathit{u}$, so, by Theorem 2.6 together with
the fact that ${\lVert C\rVert}\approx{\lVert c\rVert}$, the backward error in
the coefficients is bounded by ${\lVert\delta c\rVert}\lesssim{\lVert
c\rVert}^{2}\mathit{u}$.

######  Remark 2.2.

This backward error can be reduced by partially balancing the matrix $C$.
Suppose that, instead of computing the eigenvalues of $C$, we compute the
eigenvalues of

$\displaystyle\widetilde{C}=\left(\begin{array}[]{ccccc}0&\frac{1}{\sqrt{2}}&&&\\\
\frac{1}{\sqrt{2}}&0&\frac{1}{2}&\\\ &\frac{1}{2}&\ddots&\ddots&\\\
&&\ddots&0&\frac{{\lVert c\rVert}^{\frac{1}{2}}}{2}\\\ &&&\frac{1}{2{\lVert
c\rVert}^{\frac{1}{2}}}&0\end{array}\right)-\frac{1}{2{\lVert
c\rVert}^{\frac{1}{2}}}e_{n}\left(\begin{array}[]{ccccc}c_{0}\sqrt{2}&c_{1}&c_{2}&\cdots&{\lVert
c\rVert}^{\frac{1}{2}}c_{n-1}\end{array}\right).$ (54)

Provided that the entry in the $(n,n)$-position is small, we have that
${\lVert\widetilde{C}\rVert}\approx{\lVert c\rVert}^{\frac{1}{2}}$, so the
backward error $\delta\widetilde{C}$ of unstructured QR is bounded by
${\lVert\delta\widetilde{C}\rVert}\lesssim{\lVert
c\rVert}^{\frac{1}{2}}\mathit{u}$. In practice, it turns out that
${\lVert\delta A\rVert}\approx{\lVert\delta\widetilde{C}\rVert}\lesssim{\lVert
c\rVert}^{\frac{1}{2}}$, so ${\lVert\delta c\rVert}\lesssim{\lVert
c\rVert}^{\frac{3}{2}}\mathit{u}$. The assumption that the $(n,n)$-th element
is small is not always satisfied. However, usually the norm of $c$ is large
because the last coefficient $a_{n}$ in the non-monic expansion (46) is small.
When this is the case, we can simply raise the order of the expansion by one
by taking an additional term. The last two terms will both be small and
roughly the same size, making the $(n,n)$-th element small. Notice also that,
by adjusting the last row, the matrix $\widetilde{C}$ can be represented as a
symmetric tridiagonal matrix of magnitude ${\lVert c\rVert}^{\frac{1}{2}}$
plus a rank-1 matrix of magnitude ${\lVert c\rVert}^{\frac{1}{2}}$.

######  Remark 2.3.

Let $\texttt{bal}(C)$ denote the matrix $C$ after complete balancing.
Remarkably, in some situations, ${\lVert\texttt{bal}(C)\rVert}\approx 1$ even
when ${\lVert c\rVert}$ is large. Thus, complete balancing can completely
eliminate large entries in $C$, at the expense of destroying its symmetric
tridiagonal plus rank-1 structure. See Remark 5.1, as well as the paper [29],
for a more detailed discussion.

######  Remark 2.4.

While unstructured QR, applied to the colleague matrix, is known to achieve
only the backward error bound ${\lVert\delta c\rVert}\lesssim{\lVert
c\rVert}^{2}\mathit{u}$, the QZ algorithm, applied to an appropriately scaled
matrix pencil, does result in a backward stable rootfinder with the bound
${\lVert\delta c\rVert}\lesssim{\lVert c\rVert}\mathit{u}$ (see, for example,
[26]). This is because the eigenvalue problem for the colleague matrix can be
written as a matrix pencil $A-\lambda B$, where both $A$ and $B$ are small,
and a backward stable QZ algorithm applied to the pencil computes the exact
eigenvalues of a perturbed pencil $(A+\delta A)-\lambda(B+\delta B)$, where
${\lVert\delta A\rVert}\lesssim{\lVert A\rVert}\mathit{u}$ and ${\lVert\delta
B\rVert}\lesssim{\lVert B\rVert}\mathit{u}$. Unfortunately, it appears to be
very difficult to construct structured, $O(n^{2})$ QZ algorithms for colleague
matrices that retain the nice stability properties of the unstructured,
$O(n^{3})$ QZ algorithm.

The following theorem, stated in a slightly different form in [27], says that
if the eigenvalues of the colleague matrix $C=A+e_{n}q^{*}$ are computed using
a componentwise backward stable algorithm, then linearization is backward
stable as a rootfinding algorithm. It follows immediately from Theorem 2.6.

######  Theorem 2.7.

Suppose that the eigenvalues of the colleague matrix $C=A+e_{n}q^{*}$ are
computed by a componentwise backward stable algorithm, in the sense that the
computed eigenvalues are the exact eigenvalues of the matrix

$\displaystyle C+\delta C=A+\delta A+(e_{n}+\delta e_{n})(q+\delta q)^{*},$
(55)

where ${\lVert\delta A\rVert}\lesssim{\lVert
A\rVert}\mathit{u}\approx\mathit{u}$, ${\lVert\delta
e_{n}\rVert}\lesssim{\lVert e_{n}\rVert}\mathit{u}\approx\mathit{u}$, and
${\lVert\delta q\rVert}\lesssim{\lVert q\rVert}\mathit{u}$. Then,
linearization is backward stable as a rootfinding algorithm, with
${\lVert\delta c\rVert}\lesssim{\lVert c\rVert}\mathit{u}$.

###  2.5 Conventions

It was pointed out to the authors that, while the Hessenberg matrices in this
manuscript are all lower Hessenberg, the standard convention in numerical
linear algebra is to study the transpose of the problem, and consider only
upper Hessenberg matrices (see [18]). The upper Hessenberg form is much better
notationally since, in upper Hessenberg form, the first elimination step
eliminates the entry in the $(2,1)$-position, while, in lower Hessenberg form,
the entry in the $(n-1,n)$-position is eliminated first. Furthermore, the
upper Hessenberg form is more convenient when representing polynomials in the
Lagrange basis (see, for example, [17]). Unfortunately, at the time that this
was all pointed out, most of the writing and numerical codes were complete,
and had been written in lower Hessenberg form because of a historical fluke
related to the structure of old explicit QR codes that were used as a template
for our algorithm.

##  3 The Algorithm

In this section, we give an overview of our algorithm. We begin by describing
the class of matrices our algorithm can be applied to. Let
$\mathcal{F}_{n}\subset\mathbbm{C}^{n\times n}$ be the set of lower Hessenberg
matrices of the form

$\displaystyle A+pq^{*},$ (56)

where $A\in C^{n\times n}$ is Hermitian and $p,q^{*}\in\mathbbm{C}^{n}$.
Eidelman, Gemignani, and Gohberg observed in [20] that the matrix $A$ is
determined entirely by:

  1. 1.

The diagonal entries $d_{i}=a_{i,i}$, for $i=1,2,\ldots,n$;

  2. 2.

The superdiagonal entries $\beta_{i}=a_{i,i+1}$ for $i=1,2,\ldots,n-1$;

  3. 3.

The vectors $p$ and $q$

(see Lemma 2.1). Following [20], we call these four vectors the basic elements
or generators of $A$. In [20], the authors construct an implicit QR algorithm
that takes advantage of this structure to achieve a cost of $O(n^{2})$. They
also prove that their algorithm is backward stable, in the sense that, when
the algorithm is used to compute the eigenvalues of a matrix
$C\in\mathcal{F}_{n}$, the computed eigenvalues are the exact eigenvalues of
$C+\delta C$, where ${\lVert\delta C\rVert}\lesssim{\lVert
C\rVert}\mathit{u}$. This is the same backward stability bound that is
provided by an unstructured QR algorithm.

_In this manuscript, we describe a new explicit QR algorithm for matrices
$A+pq^{*}\in\mathcal{F}_{n}$, that also has the cost $O(n^{2})$, and prove
that our algorithm is _componentwise_ backward stable, in the sense that the
computed eigenvalues are the exact eigenvalues of $(A+\delta A)+(p+\delta
p)(q+\delta q)^{*}$, where ${\lVert\delta A\rVert}\lesssim{\lVert
A\rVert}\mathit{u}$, ${\lVert\delta p\rVert}\lesssim{\lVert
p\rVert}\mathit{u}$, and ${\lVert\delta q\rVert}\lesssim{\lVert
q\rVert}\mathit{u}$. _

To motivate our algorithm, consider first the naive unshifted QR algorithm in
exact arithmetic, applied to the matrix $C=A+pq^{*}$. Let the matrix $U_{n}\in
C^{n\times n}$ be the unitary matrix that rotates the $(n-1,n)$-plane so that

$\displaystyle(U_{n}C)_{n-1,n}=0,$ (57)

eliminating the superdiagonal in the $(n-1,n)$-th position. Likewise, let
$U_{n-1}\in C^{n\times n}$ denote the unitary matrix rotating the
$(n-2,n-1)$-plane so that

$\displaystyle(U_{n-1}U_{n}C)_{n-2,n-1}=0,$ (58)

eliminating the superdiagonal in the $(n-2,n-1)$-th position. Continuing in
this fashion, let $U_{n-2},U_{n-3},\ldots,U_{2}$ be the unitary matrices
eliminating the superdiagonal entries in the
$(n-3,n-2),\allowbreak(n-4,n-3),\allowbreak\ldots,(1,2)$ positions of the
matrices
$(U_{n-1}U_{n}C),(U_{n-2}U_{n-1}U_{n}C),\allowbreak\ldots,\allowbreak(U_{3}\cdots
U_{n-1}U_{n}C)$, respectively. Letting $U=U_{2}U_{3}\cdots U_{n}$, we have
that $UC$ is lower triangular. This matrix has the form

$\displaystyle UC=B+(Up)q^{*},$ (59)

where the upper Hessenberg part of the matrix $B=UA$ is determined entirely
by:

  1. 1.

The diagonal entries $\underline{d}_{i}=b_{i,i}$, for $i=1,2,\ldots,n$;

  2. 2.

The subdiagonal entries $\underline{\gamma}_{i}=b_{i+1,i}$, for
$i=1,2,\ldots,n-1$;

  3. 3.

The vectors $\underline{p}=Up$ and $q$

(see Lemma 2.2). Like with the matrix $A$, we call these four vectors the
basic elements or generators of (the upper Hessenberg part of) $B$.

Next, the matrix is multiplied by $U^{*}$ on the right; clearly,

$\displaystyle UCU^{*}=BU^{*}+Up(Uq)^{*},$ (60)

so

$\displaystyle UCU^{*}=UAU^{*}+Up(Uq)^{*}.$ (61)

It’s easy to show that, since $UC$ is lower triangular and
$U^{*}=U_{n}^{*}U_{n-1}^{*}\cdots U_{2}^{*}$, where $U_{k}$ rotates the
$(k-1,k)$-plane, the matrix $UCU^{*}$ is lower Hessenberg. Thus,
$UCU^{*}\in\mathcal{F}_{n}$, and the matrix $\underline{A}=UAU^{*}$ is
determined entirely by its diagonal and superdiagonal, together with
$\underline{p}=Up$ and $\underline{q}=Uq$. Furthermore, the upper triangular
part of $UAU^{*}$ is determined entirely by the upper Hessenberg part of $B$
(see Lemma 2.3), and since $UAU^{*}$ is Hermitian, it follows that the whole
of the matrix $UAU^{*}$ is determined entirely by the upper Hessenberg part of
$B$.

In our algorithm, we use only the basic elements of $A$ and $B$ to represent
our matrices. This results in a single iteration of our QR algorithm requiring
$O(n)$ operations. Furthermore, we prove that the matrix
$\widehat{\underline{A}}$ and vectors $\widehat{\underline{p}}$ and
$\widehat{\underline{q}}$, computed by a single iteration of our QR algorithm,
have the componentwise forward error bounds
${\lVert\widehat{\underline{A}}-\underline{A}\rVert}\lesssim{\lVert
A\rVert}\mathit{u}$,
${\lVert\widehat{\underline{p}}-\underline{p}\rVert}\lesssim{\lVert
p\rVert}\mathit{u}$, and
${\lVert\widehat{\underline{q}}-\underline{q}\rVert}\lesssim{\lVert
p\rVert}\mathit{u}$. We then show that these componentwise forward error
bounds result in componentwise backward stability.

###  3.1 Eliminating the Superdiagonal

In this section, we describe how our algorithm performs a single elimination
of a superdiagonal element (see Algorithm 1). Suppose that we have already
eliminated the superdiagonal elements in the positions
$(n-1,n),(n-2,n-1),\ldots,(k,k+1)$. Let $p^{(k+1)}=U_{k+1}U_{k+2}\cdots
U_{n}p$ and $B^{(k+1)}=U_{k+1}U_{k+2}\cdots U_{n}A$. Suppose further that
$\widehat{p}^{\,(k+1)}$ and $\widehat{B}^{(k+1)}$ are the computed
approximations to $p^{(k+1)}$ and $B^{(k+1)}$, and that the upper Hessenberg
part of the computed matrix $\widehat{B}^{(k+1)}$ is represented by its
generators:

  1. 1.

The diagonal elements $\widehat{d}^{\;(k+1)}_{i}=\widehat{b}^{\,(k+1)}_{i,i}$,
for $i=1,2,\ldots,n$;

  2. 2.

The superdiagional elements
$\widehat{\beta}^{\,(k+1)}_{i}=\widehat{b}^{\,(k+1)}_{i,i+1}$, for
$i=1,2,\ldots,k-1$;

  3. 3.

The subdiagional elements
$\widehat{\gamma}^{(k+1)}_{i}=\widehat{b}^{\,(k+1)}_{i+1,i}$, for
$i=1,2,\ldots,n-1$;

  4. 4.

The vectors $\widehat{p}^{\,(k+1)}$ and $q$, from which the remaining elements
in the upper Hessenberg part are inferred.

Suppose that ${\lVert\widehat{B}^{(k+1)}-B^{(k+1)}\rVert}_{H}\lesssim{\lVert
A\rVert}\mathit{u}$ and
${\lVert\widehat{p}^{\,(k+1)}-p^{(k+1)}\rVert}\lesssim{\lVert
p\rVert}\mathit{u}$. Notice that, if we define
$\widehat{B}^{(n+1)}=B^{(n+1)}=A$ and $\widehat{p}^{\,(n+1)}=p^{(n+1)}=p$,
then this is obviously true for $k=n$.

$\displaystyle\left(\begin{array}[]{ccccccc}\ddots&\widehat{\gamma}_{k-3}^{(k+1)}&\widehat{d}_{k-2}^{\;(k+1)}&\widehat{\beta}_{k-2}^{\,(k+1)}&-\widehat{p}_{k-2}^{\,(k+1)}q_{k}^{*}&-\widehat{p}_{k-2}^{\,(k+1)}q_{k+1}^{*}&\cdots\\\
\hline\cr\cdots&\times&\widehat{\gamma}_{k-2}^{(k+1)}&\widehat{d}_{k-1}^{\;(k+1)}&\widehat{\beta}_{k-1}^{\,(k+1)}&-\widehat{p}_{k-1}^{\,(k+1)}q_{k+1}^{*}&\cdots\\\
\cdots&\times&\widehat{b}_{k,k-2}^{\,(k+1)}&\widehat{\gamma}_{k-1}^{(k+1)}&\widehat{d}_{k}^{\;(k+1)}&-\widehat{p}_{k}^{\,(k+1)}q_{k+1}^{*}&\cdots\\\
\hline\cr\cdots&\times&\times&\times&\widehat{\gamma}_{k}^{(k+1)}&\widehat{d}_{k+1}^{\;(k+1)}&\ddots\end{array}\right)$

Figure 1: The $(k-1)$-th and $k$-th rows of $\widehat{B}^{(k+1)}$, represented
by its generators.

To eliminate the superdiagonal element in the $(k-1,k)$ position of
$\widehat{B}^{(k+1)}+\widehat{p}^{\,(k+1)}q^{*}$, we first compute the
rotation matrix $Q_{k}\in\mathrm{SU}(2)$ that eliminates it by a rotation in
the $(k-1,k)$-plane (see Line 4 of Algorithm 1). Next, we apply the rotation
matrix separately to the generators of $\widehat{B}^{(k+1)}$ and to the vector
$\widehat{p}^{\,(k+1)}$. Since we are only interested in computing the upper
Hessenberg part of $\widehat{B}^{(k)}$, we need to update the subdiagonal
element in the $(k-1,k-2)$ position of $\widehat{B}^{(k+1)}$, represented by
$\widehat{\gamma}^{(k+1)}_{k-2}$ (see Figure 1). However, this calculation
requires the sub-subdiagonal entry in the $(k,k-2)$ position of
$\widehat{B}^{(k+1)}$, which is unknown to us since only the upper Hessenberg
part of $\widehat{B}^{(k+1)}$ is available. Fortunately, it can be recovered
by the following trick.

Since $B^{(k+1)}=U_{k+1}U_{k+2}\cdots U_{n}A$ and $A$ is Hermitian, it follows
that $B^{(k+1)}U_{n}^{*}U_{n-1}^{*}\cdots U_{k+1}^{*}$ is also Hermitian.
Thus,

$\displaystyle(B^{(k+1)}U_{n}^{*}U_{n-1}^{*}\cdots
U_{k+1}^{*})_{k,k-2}=\overline{(B^{(k+1)}U_{n}^{*}U_{n-1}^{*}\cdots
U_{k+1}^{*})}_{k-2,k}.$ (62)

Furthermore, since right-multiplication by $U_{j}^{*}$ only affects columns
$j$ and $j-1$ (see Figure 1), we have that right-multiplication by
$U_{n}^{*}U_{n-1}^{*}\cdots U_{k+1}^{*}$ leaves $b_{k,k-2}^{(k+1)}$ unchanged.
Therefore,

$\displaystyle
b^{(k+1)}_{k,k-2}=\overline{(B^{(k+1)}U_{n}^{*}U_{n-1}^{*}\cdots
U_{k+1}^{*})}_{k-2,k}.$ (63)

We know that the entries in the $(k-2,k),(k-2,k+1),\ldots,(k-2,n)$ positions
of $\widehat{B}^{(k+1)}$ are inferred from $\widehat{p}^{\,(k+1)}$ and $q$ by
the formula

$\displaystyle\widehat{b}^{\,(k+1)}_{k-2,\ell}=-\widehat{p}_{k-2}^{\,(k+1)}q_{\ell}^{*},$
(64)

for $\ell=k,k+1,\ldots,n$. Combining (63) and (64), we thus have that the sub-
subdiagonal entry in the $(k,k-2)$ position of $\widehat{B}^{(k+1)}$ can be
recovered by the formula

$\displaystyle\widehat{b}^{\,(k+1)}_{k,k-2}=\overline{(-\widehat{p}_{k-2}^{\,(k+1)}q^{*}U_{n}^{*}U_{n-1}^{*}\cdots
U_{k+1}^{*})}_{k}.$ (65)

Defining $\widetilde{q}^{\,(k+1)}=U_{k+1}U_{k+2}\cdots U_{n}q$, we have

$\displaystyle\widehat{b}^{\,(k+1)}_{k,k-2}=-\widetilde{q}_{k}^{\,(k+1)}\widehat{p}^{\,(k+1)*}_{k-2}.$
(66)

By computing the vector $\widehat{\widetilde{q}}^{\,(k+1)}$ (see Line 14 of
Algorithm 1), we use this formula to recover the sub-subdiagonal element
$\widehat{b}^{\,(k+1)}_{k,k-2}$.

Thus, the element in the $(k-1,k-2)$ position of $\widehat{B}^{(k+1)}$,
represented by $\widehat{\gamma}_{k-2}^{(k+1)}$, is updated in Line 6 of
Algorithm 1. Next, the elements in the $(k-1,k-1)$ and $(k,k-1)$ positions of
$\widehat{B}^{(k+1)}$, represented by $\widehat{d}_{k-1}^{\;(k+1)}$ and
$\widehat{\gamma}_{k-1}^{(k+1)}$, respectively, are updated in a
straightforward way in Line 8. Finally, the elements in the $(k-1,k)$ and
$(k,k)$ positions of $\widehat{B}^{(k+1)}$, represented by
$\widehat{\beta}_{k-1}^{\,(k+1)}$ and $\widehat{d}_{k}^{\;(k+1)}$,
respectively, are updated in Line 9, and the vector $\widehat{p}^{\,(k+1)}$ is
rotated in Line 10.

Since we’ve eliminated the superdiagonal element in the $(k-1,k)$ position of
$\widehat{B}^{(k+1)}+\widehat{p}^{\,(k+1)}q^{*}$, we have that the $(k-1,k)$
element of the matrix $\widehat{B}^{(k)}$ is inferred from
$\widehat{p}^{\,(k)}$ and $q$ by the formula

$\displaystyle\widehat{b}_{k-1,k}^{\,(k)}=-\widehat{p}_{k-1}^{\,(k)}q_{k}^{*}.$
(67)

Now, we would like the upper Hessenberg part of $\widehat{B}^{(k)}$ to have a
small componentwise error, so that
${\lVert\widehat{B}^{(k)}-B^{(k)}\rVert}_{H}\lesssim{\lVert
A\rVert}\mathit{u}$. However, consider the following scenario. Suppose that
the norm of
$(\widehat{p}_{k-1}^{\,(k+1)}q_{k}^{*},\widehat{p}_{k}^{\,(k+1)}q_{k}^{*})^{T}$
is much larger than ${\lVert A\rVert}$. By Lemma 2.5, the error in
$\widehat{p}_{k-1}^{\,(k)}q_{k}^{*}$ will be approximately
$\Bigl{(}\sqrt{{\lvert\widehat{p}_{k-1}^{\,(k+1)}q_{k}^{*}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}q_{k}^{*}\rvert}^{2}}\Bigr{)}\mathit{u}$,
which will be much larger than ${\lVert A\rVert}\mathit{u}$. In this situation
then, even if ${\lVert\widehat{p}^{\,(k+1)}-p^{(k+1)}\rVert}\lesssim{\lVert
p\rVert}\mathit{u}$ and
${\lVert\widehat{B}^{(k+1)}-B^{(k+1)}\rVert}_{H}\lesssim{\lVert
A\rVert}\mathit{u}$, we will not have
${\lVert\widehat{B}^{(k)}-B^{(k)}\rVert}_{H}\lesssim{\lVert
A\rVert}\mathit{u}$. To remedy this, we must apply a correction to
$\widehat{p}_{k-1}^{\,(k)}$. Recall that the rotation matrix $Q_{k}$ was
defined to be the matrix eliminating the $(k-1,k)$-th entry of
$\widehat{B}^{(k+1)}+\widehat{p}^{\,(k+1)}q^{*}$ in exact arithmetic. If we
let $(\mathring{p}_{k-1}^{(k)},\mathring{p}_{k}^{(k)})^{T}$ denote the result
of applying $Q_{k}$ to
$(\widehat{p}_{k-1}^{\,(k+1)},\widehat{p}_{k}^{\,(k+1)})^{T}$ in exact
arithmetic, and likewise let
$(\mathring{\beta}_{k-1}^{(k)},\mathring{d}_{k}^{(k)})^{T}$ denote the result
of applying $Q_{k}$ to
$(\widehat{\beta}_{k-1}^{\,(k+1)},\widehat{d}_{k}^{\;(k+1)})^{T}$ in exact
arithmetic, then, by the definition of $Q_{k}$, we have

$\displaystyle\mathring{\beta}_{k-1}^{(k)}+\mathring{p}_{k-1}^{(k)}q_{k}^{*}=0.$
(68)

By Lemma 2.5, we have that

$\displaystyle{\lvert\mathring{\beta}_{k-1}^{(k)}-\widehat{\beta}_{k-1}^{\,(k)}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\widehat{\beta}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{d}_{k}^{\;(k+1)}\rvert}^{2}}\Bigr{)}\mathit{u}$
(69)

and

$\displaystyle{\lvert\mathring{p}_{k-1}^{(k)}q_{k}^{*}-\widehat{p}_{k-1}^{\,(k)}q_{k}^{*}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\widehat{p}_{k-1}^{\,(k+1)}q_{k}^{*}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}q_{k}^{*}\rvert}^{2}}\Bigr{)}\mathit{u}.$
(70)

Thus, if
${\lvert\widehat{p}_{k-1}^{\,(k+1)}q_{k}^{*}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}q_{k}^{*}\rvert}^{2}>{\lvert\widehat{\beta}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{d}_{k}^{\;(k+1)}\rvert}^{2}$,
then we set

$\displaystyle\widehat{p}_{k-1}^{\,(k)}q_{k}^{*}=-\widehat{\beta}_{k-1}^{\,(k)},$
(71)

so

$\displaystyle\widehat{p}_{k-1}^{\,(k)}=-\widehat{\beta}_{k-1}^{\,(k)}/q_{k}^{*}$
(72)

(see Line 12 of Algorithm 1). With this correction to
$\widehat{p}_{k-1}^{\,(k)}$, it is easy to see that
${\lVert\widehat{B}^{(k)}-B^{(k)}\rVert}_{H}\lesssim{\lVert
A\rVert}\mathit{u}$ and
${\lVert\widehat{p}^{\,(k)}-p^{(k)}\rVert}\lesssim{\lVert p\rVert}\mathit{u}$.
If, on the other hand,
${\lvert\widehat{p}_{k-1}^{\,(k+1)}q_{k}^{*}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}q_{k}^{*}\rvert}^{2}\leq{\lvert\widehat{\beta}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{d}_{k}^{\;(k+1)}\rvert}^{2}$,
then the correction is not neccessary, since in this case the error in
$\widehat{p}_{k-1}^{\,(k)}q_{k}^{*}$ is smaller than the error in
$\widehat{\beta}_{k-1}^{\,(k)}$.

This process of eliminating the superdiagonal elements can be repeated, until
the upper Hessenberg part of the matrix $\widehat{B}$, approximating
$B=U_{2}U_{3}\cdots U_{n}A$, is obtained, together with
$\widehat{\underline{p}}$, approximating $\underline{p}=U_{2}U_{3}\cdots
U_{n}p$ (see Algorithm 1). In Section 4.1, Lemma 4.1, we prove that the
forward errors in the upper Hessenberg part of $\widehat{B}$ and in the vector
$\widehat{\underline{p}}$ are proportional to ${\lVert A\rVert}\mathit{u}$ and
${\lVert p\rVert}\mathit{u}$, respectively.

Algorithm 1 (A single elimination of the superdiagonal) Inputs: This algorithm
accepts as inputs two vectors $d$ and $\beta$ representing the diagonal and
superdiagonal, respectively, of an $n\times n$ Hermitian matrix $A$, as well
as two vectors $p$ and $q$ of length $n$, where $A+pq^{*}$ is lower
Hessenberg. Outputs: It returns as its outputs the rotation matrices
$Q_{2},Q_{3},\ldots,Q_{n}\in\mathbbm{C}^{2\times 2}$ so that, letting
$U_{k}\in\mathbbm{C}^{n\times n}$, $k=2,3,\ldots,n$, denote the matrices that
rotate the $(k-1,k)$-plane by $Q_{k}$, $U_{2}U_{3}\cdots U_{n}(A+pq^{*})$ is
lower triangular. It also returns the vectors $\underline{d}$,
$\underline{\gamma}$, and $\underline{p}$, where $\underline{d}$ and
$\underline{\gamma}$ represent the diagonal and subdiagonal, respectively, of
the matrix $U_{2}U_{3}\cdots U_{n}A$, and $\underline{p}=U_{2}U_{3}\cdots
U_{n}p$.

1:Set $\gamma\leftarrow\overline{\beta}$, where $\gamma$ represents the
subdiagonal.

2:Make a copy of $q$, setting $\widetilde{q}\leftarrow q$.

3:for $k=n,n-1,\ldots,2$ do

4: Construct the $2\times 2$ rotation matrix $Q_{k}\in\mathrm{SU}(2)$ so that
$\displaystyle\Bigl{(}Q_{k}\left[\begin{array}[]{c}\beta_{k-1}+p_{k-1}q_{k}^{*}\\\
d_{k}+p_{k}q_{k}^{*}\end{array}\right]\Bigr{)}_{1}=0.$

5: if $k\neq 2$ then

6: Rotate the subdiagonal and the sub-subdiagonal:
$\displaystyle\gamma_{k-2}\leftarrow\Bigl{(}Q_{k}\left[\begin{array}[]{c}\gamma_{k-2}\\\
-\widetilde{q}_{k}p_{k-2}^{*}\end{array}\right]\Bigr{)}_{1}$

7: end if

8: Rotate the diagonal and the subdiagonal:
$\left[\begin{array}[]{c}d_{k-1}\\\ \gamma_{k-1}\end{array}\right]\leftarrow
Q_{k}\left[\begin{array}[]{c}d_{k-1}\\\ \gamma_{k-1}\end{array}\right].$

9: Rotate the superdiagonal and the diagonal:
$\left[\begin{array}[]{c}\beta_{k-1}\\\ d_{k}\end{array}\right]\leftarrow
Q_{k}\left[\begin{array}[]{c}\beta_{k-1}\\\ d_{k}\end{array}\right].$

10: Rotate $p$: $\left[\begin{array}[]{c}p_{k-1}\\\
p_{k}\end{array}\right]\leftarrow Q_{k}\left[\begin{array}[]{c}p_{k-1}\\\
p_{k}\end{array}\right]$

11: if ${\lvert p_{k-1}q_{k}^{*}\rvert}^{2}+{\lvert
p_{k}q_{k}^{*}\rvert}^{2}>{\lvert\beta_{k-1}\rvert}^{2}+{\lvert
d_{k}\rvert}^{2}$ then

12: Correct the vector $p$, setting
$p_{k-1}\leftarrow-\frac{\beta_{k-1}}{q_{k}^{*}}$

13: end if

14: Rotate $\widetilde{q}$: $\left[\begin{array}[]{c}\widetilde{q}_{\,k-1}\\\
\widetilde{q}_{\,k}\end{array}\right]\leftarrow
Q_{k}\left[\begin{array}[]{c}\widetilde{q}_{k-1}\\\
\widetilde{q}_{k}\end{array}\right]$

15:end for

16:Set $\underline{d}\leftarrow d$, $\underline{\gamma}\leftarrow\gamma$, and
$\underline{p}\leftarrow p$.

###  3.2 Rotating Back to Hessenberg Form

In this section, we describe how our algorithm rotates the triangular matrix
produced by an elimination of the superdiagonal back to lower Hessenberg form
(see Algorithm 2). Suppose that $B$ is an $n\times n$ matrix and that $p$ and
$q$ are vectors such that $B+pq^{*}$ is lower triangular. Notice that this
condition is satisfied by the matrix $\widehat{B}$ and the vectors
$\widehat{\underline{p}}$ and $q$ from the preceding section, produced by an
elimination of the superdiagonal. Let $\gamma$ denote the subdiagonal of $B$.
Suppose that we have already applied the rotation matrices
$U_{n}^{*},U_{n-1}^{*},\cdots,U_{k+1}^{*}$ to the right of $B$ and $q^{*}$,
and let $q^{(k+1)}=U_{k+1}U_{k+2}\cdots U_{n}q$ and
$A^{(k+1)}=BU_{n}^{*}U_{n-1}^{*}\cdots U_{k+1}^{*}$. Suppose that
$\widehat{q}^{\,(k+1)}$ and $\widehat{A}^{(k+1)}$ are the computed
approximations to $q^{(k+1)}$ and $A^{(k+1)}$, respectively, and that the
upper triangular part of $\widehat{A}^{(k+1)}$ is represented by its
generators:

  1. 1.

The diagonal entries $\widehat{d}^{\;(k+1)}_{i}=\widehat{a}^{\,(k+1)}_{i,i}$,
for $i=1,2,\ldots,n$;

  2. 2.

The superdiagonal entries
$\widehat{\beta}^{\,(k+1)}_{i}=\widehat{a}^{\,(k+1)}_{i,i+1}$ for
$i=k,k+1,\ldots,n-1$;

  3. 3.

The vectors $p$ and $\widehat{q}^{\,(k+1)}$, from which the remaining elements
in the upper triangular part are inferred.

Suppose that ${\lVert\widehat{A}^{(k+1)}-A^{(k+1)}\rVert}_{T}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$ and
${\lVert\widehat{q}^{\,(k+1)}-q^{(k+1)}\rVert}\lesssim{\lVert
q\rVert}\mathit{u}$. Notice that, if we define
$\widehat{A}^{(n+1)}=A^{(n+1)}=B$ and $\widehat{q}^{\,(n+1)}=q^{(n+1)}=q$,
then this is obviously true for $k=n$.

$\displaystyle\left(\begin{array}[]{c|cc|c}\ddots&\vdots&\vdots&\vdots\\\
\widehat{d}_{k-2}^{\;(k+1)}&-p_{k-2}\widehat{q}_{k-1}^{\,(k+1)*}&-p_{k-2}\widehat{q}_{k}^{\,(k+1)*}&-p_{k-2}\widehat{q}_{k+1}^{\,(k+1)*}\\\
\gamma_{k-2}&\widehat{d}_{k-1}^{\;(k+1)}&-p_{k-1}\widehat{q}_{k}^{\,(k+1)*}&-p_{k-1}\widehat{q}_{k+1}^{\,(k+1)*}\\\
\times&\gamma_{k-1}&\widehat{d}_{k}^{\;(k+1)}&\widehat{\beta}_{k}^{\,(k+1)}\\\
\times&\times&\times&\widehat{d}_{k+1}^{\;(k+1)}\\\
\vdots&\vdots&\vdots&\ddots\end{array}\right)$

Figure 2: The $(k-1)$-th and $k$-th columns of $\widehat{A}^{(k+1)}$,
represented by its generators.

To apply the matrix $U^{*}_{k}$ to
$\widehat{A}^{(k+1)}+p\widehat{q}^{\,(k+1)}$ on the right, we apply the matrix
$Q_{k}^{*}\in\mathrm{SU}(2)$ separately to the generators of
$\widehat{A}^{(k+1)}$ and to the vector $\widehat{q}^{\,(k+1)}$. We start by
rotating the diagonal and superdiagonal elements in the $(k-1,k-1)$ and
$(k-1,k)$ positions of $\widehat{A}^{(k+1)}$, represented by
$\widehat{d}_{k-1}^{\;(k+1)}$ and $-p_{k-1}\widehat{q}_{k}^{\,(k+1)*}$,
respectively, in Line 2 of Algorithm 2, saving the superdiagonal element in
$\widehat{\beta}_{k-1}^{\,(k)}$ (see Figure 2). Next, we rotate the elements
in the $(k,k-1)$ and $(k,k)$ positions, represented by $\gamma_{k-1}$ (the
$(k-1)$-st element of the subdiagonal of $B$) and $\widehat{d}_{k}^{\;(k+1)}$,
respectively, in a straightforward way in Line 3; since we are only interested
in computing the upper triangular part of $\widehat{A}^{(k)}$, we only update
the diagonal entry. Finally, we rotate the vector $\widehat{q}^{\,(k+1)}$ in
Line 4.

The process of applying the rotation matrices on the right can be repeated,
until the upper triangular part of the matrix $\widehat{\underline{A}}$,
approximating $\underline{A}=BU_{n}^{*}U_{n-1}^{*}\cdots U_{2}^{*}$, is
obtained, together with $\widehat{\underline{q}}$, approximating
$\underline{q}=U_{2}U_{3}\cdots U_{n}q$ (see Algorithm 2). In Section 4.1,
Lemma 4.2, we prove that the forward errors in the upper triangular part of
$\widehat{\underline{A}}$ and in the vector $\widehat{\underline{q}}$ are
proportional to ${\lVert B\rVert}_{H}\mathit{u}$ and ${\lVert
q\rVert}\mathit{u}$, respectively.

###  3.3 The QR Algorithms

The elimination of the superdiagonal described in Algorithm 1, followed by the
rotation back to Hessenberg form described in Algorithm 2, can be iterated to
find the eigenvalues of $A+pq^{*}$. Our unshifted explicit QR algorithm, based
on this iteration, is described in Algorithm 3. This unshifted QR algorithm
can be accelerated by the introduction of shifts; our explicit shifted QR
algorithm, with Wilkinson shifts, is described in Algorithm 4.

In Section 4.1, we show that the forward error of one iteration of our QR
algorithm (Algorithm 1 followed by Algorithm 2) satisfies componentwise
forward error bounds. In Section 4.2, we use this result to prove that both
our explicit unshifted QR algorithm (Algorithm 3) and our shifted QR algorithm
(Algorithm 4) are componentwise backward stable.

Algorithm 2 (Rotating the matrix back to Hessenberg form) Inputs: This
algorithm accepts as inputs $n-1$ rotation matrices
$Q_{2},Q_{2},\ldots,Q_{n}\in\mathbbm{C}^{n\times n}$, two vectors $d$ and
$\gamma$ representing the diagonal and subdiagonal, respectively, of an
$n\times n$ complex matrix $B$, and two vectors $p$ and $q$ of length $n$,
where $B+pq^{*}$ is lower triangular. Outputs: Letting
$U_{k}\in\mathbbm{C}^{n\times n}$, $k=2,3,\ldots,n$, denote the matrices that
rotate the $(k-1,k)$-plane by $Q_{k}$, this algorithm returns as its outputs
the vectors $\underline{d}$, $\underline{\beta}$, and $\underline{q}$, where
$\underline{d}$ and $\underline{\beta}$ represent the diagonal and
superdiagonal, respectively, of the matrix $BU_{n}^{*}U_{n-1}^{*}\cdots
U_{2}^{*}$, and $\underline{q}=U_{2}U_{3}\cdots U_{n}q$.

1:for $k=n,n-1,\ldots,2$ do

2: Rotate the diagonal and the superdiagonal:
$\displaystyle\left[\begin{array}[]{c}d_{k-1}\\\
\beta_{k-1}\end{array}\right]\leftarrow\overline{Q_{k}}\left[\begin{array}[]{c}d_{k-1}\\\
-p_{k-1}q_{k}^{*}\end{array}\right].$

3: Rotate the subdiagonal and the diagonal:  $\displaystyle
d_{k}\leftarrow\Bigl{(}\overline{Q_{k}}\left[\begin{array}[]{c}\gamma_{k-1}\\\
d_{k}\end{array}\right]\Bigr{)}_{2}$

4: Rotate $q$: $\left[\begin{array}[]{c}q_{k-1}\\\
q_{k}\end{array}\right]\leftarrow Q_{k}\left[\begin{array}[]{c}q_{k-1}\\\
q_{k}\end{array}\right]$

5:end for

6:Set $\underline{d}\leftarrow d$, $\underline{\beta}\leftarrow\beta$, and
$\underline{q}\leftarrow q$.

Algorithm 3 (Unshifted explicit QR) Inputs: This algorithm accepts as inputs
two vectors $d$ and $\beta$ representing the diagonal and superdiagonal,
respectively, of an $n\times n$ Hermitian matrix $A$, as well as two vectors
$p$ and $q$ of length $n$, where $A+pq^{*}$ is lower Hessenberg. It also
accepts a tolerance $\epsilon>0$, which determines the accuracy the
eigenvalues are computed to. Outputs: It returns as its output the vector
$\lambda$ of length $n$ containing the eigenvalues of the matrix $A+pq^{*}$.

1:for $i=1,2,\ldots,n-1$ do

2: while $\beta_{i}+p_{i}q_{i+1}^{*}\geq\epsilon$ do $\triangleright$ Check if
$(A+pq^{*})_{i,i+1}$ is close to zero

3: Perform one iteration of QR (one step of Algorithm 1 followed by one step
of Algorithm 2) on the submatrix $(A+pq^{*})_{i:n,i:n}$ defined by the vectors
$d_{i:n}$, $\beta_{i:n-1}$, $p_{i:n}$, and $q_{i:n}$.

4: end while

5:end for

6:Set $\lambda\leftarrow d$.

Algorithm 4 (Shifted explicit QR) Inputs: This algorithm accepts as inputs two
vectors $d$ and $\beta$ representing the diagonal and superdiagonal,
respectively, of an $n\times n$ Hermitian matrix $A$, as well as two vectors
$p$ and $q$ of length $n$, where $A+pq^{*}$ is lower Hessenberg. It also
accepts a tolerance $\epsilon>0$, which determines the accuracy the
eigenvalues are computed to. Outputs: It returns as its output the vector
$\lambda$ of length $n$ containing the eigenvalues of the matrix $A+pq^{*}$.

1:for $i=1,2,\ldots,n-1$ do

2: Set $\mu_{\text{sum}}\leftarrow 0$.

3: while $\beta_{i}+p_{i}q_{i+1}^{*}\geq\epsilon$ do $\triangleright$ Check if
$(A+pq^{*})_{i,i+1}$ is close to zero

4: Compute the eigenvalues $\mu_{1}$ and $\mu_{2}$ of the $2\times 2$
submatrix
$\left[\begin{array}[]{cc}d_{i}+p_{i}q_{i}^{*}&\beta_{i}+p_{i}q_{i+1}^{*}\\\
\overline{\beta}_{i}+p_{i+1}q_{i}^{*}&d_{i+1}+p_{i+1}q_{i+1}^{*}\end{array}\right].$
$\triangleright$ This is just ${(A+pq^{*})_{i:i+1,i:i+1}}$

5: Set $\mu$ to whichever of $\mu_{1}$ and $\mu_{2}$ is closest to
$d_{i}+p_{i}q_{i}^{*}$.

6: Set $\mu_{\text{sum}}\leftarrow\mu_{\text{sum}}+\mu$.

7: Set $d_{i:n}\leftarrow d_{i:n}-\mu$.

8: Perform one iteration of QR (one step of Algorithm 1 followed by one step
of Algorithm 2) on the submatrix $(A+pq^{*})_{i:n,i:n}$ defined by the vectors
$d_{i:n}$, $\beta_{i:n-1}$, $p_{i:n}$, and $q_{i:n}$.

9: end while

10: Set $d_{i:n}\leftarrow d_{i:n}+\mu_{\text{sum}}$.

11:end for

12:Set $\lambda_{i}\leftarrow d_{i}+p_{i}q_{i}^{*}$, for $i=1,2,\ldots,n$.

##  4 Componentwise Backward Stability

The principal results of this section are Theorems 4.6 and 4.7, which state
that our unshifted and shifted QR algorithms, respectively, are componentwise
backward stable. In Section 4.1, we prove that the forward error of a single
sweep of our unshifted QR algorithm satisfies componentwise bounds. In Section
4.2, we use these bounds show componentwise backward stability of our QR
algorithms.

###  4.1 Forward Error Analysis of a Single Sweep of $QR$

Suppose that $A$ is Hermitian and $A+pq^{*}$ is lower Hessenberg. In this
section, we prove in Theorem 4.3 that the forward errors in $A$, $p$, and $q$
of single sweep of our explicit $QR$ algorithm are proportional to ${\lVert
A\rVert}\mathit{u}$, ${\lVert p\rVert}\mathit{u}$, and ${\lVert
q\rVert}\mathit{u}$, respectively.

The following lemma bounds the forward error of Algorithm 1 (the elimination
of the superdiagonal).

######  Lemma 4.1.

Suppose that $A\in\mathbbm{C}^{n\times n}$ is a Hermitian matrix, and that
$p,q\in\mathbbm{C}^{n}$. Suppose further that $A+pq^{*}$ is lower Hessenberg,
and let $d$ and $\beta$ denote the diagonal and superdiagonal of $A$,
respectively. Suppose that Algorithm 1 is carried out in floating point
arithmetic with $d$, $\beta$, $p$, and $q$ as inputs, and let
$Q_{2},Q_{3},\ldots,Q_{n}\in\mathrm{SU}(2)$ be the unitary matrices generated
by an exact step of Line 4 of Algorithm 1 applied to the computed vectors at
that step. Let $U_{k}\in\mathbbm{C}^{n\times n}$, $k=2,3,\ldots,n$, denote the
matrices that rotate the $(k-1,k)$-plane by $Q_{k}$, and define
$U\in\mathbbm{C}^{n\times n}$ by the formula $U=U_{2}U_{3}\cdots U_{n}$.
Suppose finally that $\underline{\widehat{d}}$,
$\underline{\widehat{\gamma}}$, and $\underline{\widehat{p}}$ are the outputs
generated by Algorithm 1, and define the upper Hessenberg part of the matrix
${\widehat{B}}\in\mathbbm{C}^{n\times n}$ by the formula

$\displaystyle{\widehat{b}}_{i,j}=\left\\{\begin{array}[]{ll}-\underline{\widehat{p}}_{i}{q}_{j}^{*}&\text{if
$j>i$},\\\ \underline{\widehat{d}}_{i}&\text{if $j=i$},\\\
\underline{\widehat{\gamma}}_{j}&\text{if $j=i-1$}.\end{array}\right.$ (76)

where $\widehat{b}_{i,j}$ denotes the $(i,j)$-th entry of $\widehat{B}$. Let
$B=UA$ and $\underline{p}=Up$. Then

$\displaystyle{\lVert{\widehat{B}}-B\rVert}_{H}\lesssim{\lVert
A\rVert}\mathit{u}$ (77)

and

$\displaystyle{\lVert\underline{\widehat{p}}-\underline{p}\rVert}\lesssim{\lVert
p\rVert}\mathit{u},$ (78)

where ${\lVert\cdot\rVert}_{H}$ denotes the square root of the sum of squares
of the entries in the upper Hessenberg part of its argument (see Definition
2.1).

Proof. Suppose that $\widehat{d}^{\,(k)}$, $\widehat{\gamma}^{(k)}$,
$\widehat{\beta}^{\,(k)}$, $\widehat{p}^{\,(k)}$, and
$\widehat{\widetilde{q}}^{\,(k)}$ denote the computed vectors in Algorithm 1
after the elimination of the superdiagonal elements in the positions
$(n-1,n),(n-2,n-1),\ldots,(k-1,k)$. Suppose further that the upper Hessenberg
part of the matrix $\widehat{B}^{(k)}\in\mathbbm{C}^{n\times n}$ is defined by
the formula

$\displaystyle\widehat{b}^{\,(k)}_{i,j}=\left\\{\begin{array}[]{ll}-\widehat{p}_{i}^{\,(k)}q_{j}^{*}&\text{if
$j>i+1$ or if $j=i+1$ and $j\geq k$},\\\ \widehat{\beta}_{i}^{\,(k)}&\text{if
$j=i+1$ and $j<k$},\\\ \widehat{d}_{i}^{\;(k)}&\text{if $j=i$},\\\
\widehat{\gamma}_{i}^{\,(k)}&\text{if $j=i-1$},\end{array}\right.$ (83)

where $\widehat{b}_{i,j}^{\,(k)}$ denotes the $(i,j)$-th entry of
$\widehat{B}^{(k)}$. Clearly, $\widehat{\underline{d}}=\widehat{d}^{\,(2)}$,
$\widehat{\underline{\gamma}}=\widehat{\gamma}^{(2)}$,
$\widehat{\underline{p}}=\widehat{p}^{\,(2)}$, and
$\widehat{B}=\widehat{B}^{(2)}$. Let $B^{(k)}=U_{k}U_{k+1}\cdots U_{n}A$ and
$p^{(k)}=U_{k}U_{k+1}\cdots U_{n}p$. We will prove that
${\lVert\widehat{B}^{(k)}-B^{(k)}\rVert}_{H}\lesssim{\lVert
A\rVert}\mathit{u}$ and
${\lVert\widehat{p}^{\,(k)}-p^{(k)}\rVert}\lesssim{\lVert p\rVert}\mathit{u}$,
for each $k=n,n-1,\ldots,2$.

We begin by proving this statement for $k=n$. From Line 4, we have that the
matrix $Q_{n}\in\mathrm{SU}(2)$ satisfies

$\displaystyle\biggl{(}Q_{n}\left[\begin{array}[]{c}\beta_{n-1}+p_{n-1}q_{n}^{*}\\\
d_{n}+p_{n}q_{n}^{*}\end{array}\right]\biggr{)}_{1}=0,$ (86)

with the computed matrix $\widehat{Q}_{n}$ satisfying
${\lVert\widehat{Q}_{n}-Q_{n}\rVert}\lesssim u$ by Lemma 2.4. In Line 6, we
have

$\displaystyle\widehat{\gamma}_{n-2}^{(n)}=\mathop{\mathit{fl}}\biggl{(}\widehat{Q}_{n}\left[\begin{array}[]{c}\gamma_{n-2}\\\
-\widetilde{q}_{n}p_{n-2}^{*}\end{array}\right]\biggr{)}_{1}.$ (89)

At this stage $\widetilde{q}$ is still equal to $q$ and, according to Lemma
2.1, $a_{n,n-2}=-q_{n}p_{n-2}^{*}$. By definition, $a_{n-1,n-2}=\gamma_{n-2}$.
Therefore, by Lemma 2.5, we have that
${\lvert\widehat{\gamma}_{n-2}^{(n)}-b_{n-1,n-2}^{(n)}\rvert}\lesssim{\lVert
A\rVert}u$, where $b_{i,j}^{(n)}$ denotes the $(i,j)$-th entry of $B^{(n)}$.
In Line 8, we have

$\displaystyle\left[\begin{array}[]{c}\widehat{d}_{n-1}^{\;(n)}\\\
\widehat{\gamma}_{n-1}^{(n)}\end{array}\right]=\mathop{\mathit{fl}}\biggl{(}\widehat{Q}_{n}\left[\begin{array}[]{c}d_{n-1}\\\
\gamma_{n-1}\end{array}\right]\biggr{)}.$ (94)

Since $a_{n-1,n-1}=d_{n-1}$ and $a_{n,n-1}=\gamma_{n-1}$, by Lemma 2.5, we
have that
${\lvert\widehat{d}_{n-1}^{\;(n)}-b_{n-1,n-1}^{(n)}\rvert}\lesssim{\lVert
A\rVert}u$ and
${\lvert\widehat{\gamma}_{n-1}^{(n)}-b_{n,n-1}^{(n)}\rvert}\lesssim{\lVert
A\rVert}u$. In Line 9, we have

$\displaystyle\left[\begin{array}[]{c}\widehat{\beta}_{n-1}^{\,(n)}\\\
\widehat{d}_{n}^{\;(n)}\end{array}\right]=\mathop{\mathit{fl}}\biggl{(}\widehat{Q}_{n}\left[\begin{array}[]{c}\beta_{n-1}\\\
d_{n}\end{array}\right]\biggr{)}.$ (99)

Since, by definition, $a_{n-1,n}=\beta_{n-1}$ and $a_{n,n}=d_{n}$, it follows
from Lemma 2.5 that
${\lvert\widehat{\beta}_{n-1}^{\,(n)}-b_{n-1,n}^{(n)}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\beta_{n-1}\rvert}^{2}+{\lvert
d_{n}\rvert}^{2}}\Bigr{)}\mathit{u}\leq{\lVert A\rVert}\mathit{u}$ and
${\lvert\widehat{d}_{n}^{\;(n)}-b_{n,n}^{(n)}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\beta_{n-1}\rvert}^{2}+{\lvert
d_{n}\rvert}^{2}}\Bigr{)}\mathit{u}\leq{\lVert A\rVert}\mathit{u}$. In Line
10, we have

$\displaystyle\left[\begin{array}[]{c}\widehat{p}_{n-1}^{\,(n)\dagger}\\\
\widehat{p}_{n}^{\,(n)}\end{array}\right]=\mathop{\mathit{fl}}\biggl{(}\widehat{Q}_{n}\left[\begin{array}[]{c}p_{n-1}\\\
p_{n}\end{array}\right]\biggr{)},$ (104)

where $\widehat{p}_{n-1}^{\,(n)\dagger}$ is a temporary value that will be
corrected later. Once again, Lemma 2.5 tells us that
${\lvert\widehat{p}_{n-1}^{\,(n)\dagger}-p_{n-1}^{(n)}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert
p_{n-1}\rvert}^{2}+{\lvert p_{n}\rvert}^{2}}\Bigr{)}\mathit{u}\leq{\lVert
p\rVert}\mathit{u}$ and
${\lvert\widehat{p}_{n}^{\,(n)}-p_{n}^{(n)}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert
p_{n-1}\rvert}^{2}+{\lvert p_{n}\rvert}^{2}}\Bigr{)}\mathit{u}\leq{\lVert
p\rVert}\mathit{u}$. Next, we observe that, by the definition of $Q_{n}$, we
have that

$\displaystyle b_{n-1,n}^{(n)}+p_{n-1}^{(n)}q_{n}^{*}=0.$ (105)

In Line 12, we apply a correction to $\widehat{p}_{n-1}^{\,(n)\dagger}$, so
that

$\displaystyle\widehat{p}_{n-1}^{\,(n)}=\left\\{\begin{array}[]{ll}-\widehat{\beta}_{n-1}^{\,(n)}/q_{n}^{*}&\text{if
${\lvert p_{n-1}q_{n}^{*}\rvert}^{2}+{\lvert
p_{n}q_{n}^{*}\rvert}^{2}>{\lvert\beta_{n-1}\rvert}^{2}+{\lvert
d_{n}\rvert}^{2}$},\\\ \widehat{p}_{n-1}^{\,(n)\dagger}&\text{if ${\lvert
p_{n-1}q_{n}^{*}\rvert}^{2}+{\lvert
p_{n}q_{n}^{*}\rvert}^{2}\leq{\lvert\beta_{n-1}\rvert}^{2}+{\lvert
d_{n}\rvert}^{2}$}.\end{array}\right.$ (108)

From this, we see that, if ${\lvert p_{n-1}q_{n}^{*}\rvert}^{2}+{\lvert
p_{n}q_{n}^{*}\rvert}^{2}>{\lvert\beta_{n-1}\rvert}^{2}+{\lvert
d_{n}\rvert}^{2}$, then

$\displaystyle{\lvert\widehat{p}_{n-1}^{\,(n)}-p_{n-1}^{(n)}\rvert}={\lvert-\widehat{\beta}_{n-1}^{\,(n)}/q_{n}^{*}-p_{n-1}^{(n)}\rvert}$
$\displaystyle={\lvert
b_{n-1,n}^{(n)}/q_{n}^{*}-\widehat{\beta}_{n-1}^{\,(n)}/q_{n}^{*}\rvert}$
$\displaystyle=\frac{1}{{\lvert q_{n}^{*}\rvert}}{\lvert
b_{n-1,n}^{(n)}-\widehat{\beta}_{n-1}^{\,(n)}\rvert}$
$\displaystyle\lesssim\frac{\sqrt{{\lvert\beta_{n-1}\rvert}^{2}+{\lvert
d_{n}\rvert}^{2}}}{{\lvert q_{n}^{*}\rvert}}\mathit{u}$
$\displaystyle\leq\Bigl{(}\sqrt{{\lvert p_{n-1}\rvert}^{2}+{\lvert
p_{n}\rvert}^{2}}\Bigr{)}\mathit{u}$ $\displaystyle\leq{\lVert
p\rVert}\mathit{u}.$ (109)

where the second equality is due to (105). Furthermore,

$\displaystyle{\lvert-\widehat{p}_{n-1}^{\,(n)}q_{n}^{*}-b_{n-1,n}^{(n)}\rvert}={\lvert\widehat{\beta}_{n-1}^{\,(n)}-b_{n-1,n}^{(n)}\rvert}$
$\displaystyle\lesssim\Bigl{(}\sqrt{{\lvert\beta_{n-1}\rvert}^{2}+{\lvert
d_{n}\rvert}^{2}}\Bigr{)}\mathit{u}$ $\displaystyle\leq{\lVert
A\rVert}\mathit{u}.$ (110)

If, on the other hand, ${\lvert p_{n-1}q_{n}^{*}\rvert}^{2}+{\lvert
p_{n}q_{n}^{*}\rvert}^{2}\leq{\lvert\beta_{n-1}\rvert}^{2}+{\lvert
d_{n}\rvert}^{2}$, then

$\displaystyle{\lvert\widehat{p}_{n-1}^{\,(n)}-p_{n-1}^{(n)}\rvert}={\lvert\widehat{p}_{n-1}^{\,(n)\dagger}-p_{n-1}^{(n)}\rvert}$
$\displaystyle\lesssim\sqrt{{\lvert p_{n-1}\rvert}^{2}+{\lvert
p_{n}\rvert}^{2}}\mathit{u}$ $\displaystyle\leq{\lVert p\rVert}\mathit{u}.$
(111)

Moreover,

$\displaystyle{\lvert-\widehat{p}_{n-1}^{\,(n)}q_{n}^{*}-b_{n-1,n}^{(n)}\rvert}={\lvert-\widehat{p}_{n-1}^{\,(n)\dagger}q_{n}^{*}-b_{n-1,n}^{(n)}\rvert}$
$\displaystyle={\lvert-\widehat{p}_{n-1}^{\,(n)\dagger}q_{n}^{*}+p_{n-1}^{(n)}q_{n}^{*}\rvert}$
$\displaystyle\lesssim{\lvert q_{n}^{*}\rvert}\Bigl{(}\sqrt{{\lvert
p_{n-1}\rvert}^{2}+{\lvert p_{n}\rvert}^{2}}\Bigr{)}\mathit{u}$
$\displaystyle\leq\Bigl{(}\sqrt{{\lvert\beta_{n-1}\rvert}^{2}+{\lvert
d_{n}\rvert}^{2}}\Bigr{)}\mathit{u}$ $\displaystyle\leq{\lVert
A\rVert}\mathit{u}.$ (112)

where the second equality follows from (105). This completes the proof that
${\lVert\widehat{B}^{(n)}-B^{(n)}\rVert}_{H}\lesssim{\lVert A\rVert}u$ and
${\lVert\widehat{p}^{\,(n)}-p^{(n)}\rVert}\lesssim{\lVert p\rVert}u$.

Now, we will show that, if
${\lVert\widehat{B}^{(k+1)}-B^{(k+1)}\rVert}_{H}\lesssim{\lVert A\rVert}u$ and
${\lVert\widehat{p}^{\,(k+1)}-p^{(k+1)}\rVert}\lesssim{\lVert p\rVert}u$, then
${\lVert\widehat{B}^{(k)}-B^{(k)}\rVert}_{H}\lesssim{\lVert A\rVert}u$ and
${\lVert\widehat{p}^{\,(k)}-p^{(k)}\rVert}\lesssim{\lVert p\rVert}u$. From
Line 4, we have that the matrix $Q_{k}\in\mathrm{SU}(2)$ satisfies

$\displaystyle\biggl{(}Q_{k}\left[\begin{array}[]{c}\widehat{\beta}_{k-1}^{\,(k+1)}+\widehat{p}_{k-1}^{\,(k+1)}q_{k}^{*}\\\
\widehat{d}_{k}^{\;(k+1)}+\widehat{p}_{k}^{\,(k+1)}q_{k}^{*}\end{array}\right]\biggr{)}_{1}=0,$
(115)

with the computed matrix $\widehat{Q}_{k}$ satisfying
${\lVert\widehat{Q}_{k}-Q_{k}\rVert}\lesssim\mathit{u}$ by Lemma 2.4. In Line
6, we have

$\displaystyle\widehat{\gamma}_{k-2}^{(k)}=\mathop{\mathit{fl}}\biggl{(}\widehat{Q}_{k}\left[\begin{array}[]{c}\widehat{\gamma}_{k-2}^{(k+1)}\\\
-\widehat{\widetilde{q}}_{k}^{\,(k+1)}\widehat{p}_{k-2}^{\,(k+1)*}\end{array}\right]\biggr{)}_{1}.$
(118)

We must first show that

$\displaystyle{\bigl{\lvert}-\widehat{\widetilde{q}}_{k}^{\,(k+1)}\widehat{p}_{k-2}^{\,(k+1)*}-b_{k,k-2}^{(k+1)}\bigr{\rvert}}\lesssim{\lVert
A\rVert}\mathit{u}.$ (119)

We begin by observing that

$\displaystyle
b_{k,k-2}^{(k+1)}=(B^{(k+1)})_{k,k-2}=(B^{(k+1)}U_{n}^{*}U_{n-1}^{*}\cdots
U_{k+1}^{*})_{k,k-2},$ (120)

since right-multiplication by $U_{j}^{*}$ only affects columns $j$ and $j-1$.
We now observe that

$\displaystyle B^{(k+1)}U_{n}^{*}U_{n-1}^{*}\cdots
U_{k+1}^{*}=U_{k+1}U_{k+2}\cdots U_{n}AU_{n}^{*}U_{n-1}^{*}\cdots U_{k+1}^{*}$
(121)

is Hermitian, so from (120) we have that

$\displaystyle
b_{k,k-2}^{(k+1)}=\overline{(B^{(k+1)}U_{n}^{*}U_{n-1}^{*}\cdots
U_{k+1}^{*})}_{k-2,k}.$ (122)

By the induction hypothesis,

$\displaystyle\widehat{b}_{k-2,\ell}^{\,(k+1)}=-\widehat{p}_{k-2}^{\,(k+1)}q_{\ell}^{*}$
(123)

and

$\displaystyle{\lvert\widehat{b}_{k-2,\ell}^{\,(k+1)}-b_{k-2,\ell}^{(k+1)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u},$ (124)

for all $\ell=k,k+1,\ldots,n$. Thus,

$\displaystyle{\bigl{\lvert}(-\widehat{p}_{k-2}^{\,(k+1)}q^{*}U_{n}^{*}U_{n-1}^{*}\cdots
U_{k+1}^{*})_{k}-(B^{(k+1)}U_{n}^{*}U_{n-1}^{*}\cdots
U_{k+1}^{*})_{k-2,k}\bigr{\rvert}}\lesssim{\lVert A\rVert}\mathit{u}.$ (125)

Combining (125) with (122),

$\displaystyle{\bigl{\lvert}-\widehat{p}_{k-2}^{\,(k+1)}\widetilde{q}_{k}^{\,(k+1)*}-\overline{b_{k,k-2}^{(k+1)}}\bigr{\rvert}}\lesssim{\lVert
A\rVert}\mathit{u},$ (126)

where $\widetilde{q}^{\,(k+1)}=U_{k+1}U_{k+2}\cdots U_{n}q$. From Line 14 we
have

$\displaystyle\widehat{\widetilde{q}}^{\,(k+1)}=\mathop{\mathit{fl}}(\widehat{U}_{k+1}\widehat{U}_{k+2}\cdots\widehat{U}_{n}q),$
(127)

and, by repeated application of Lemma 2.5,

$\displaystyle{\bigl{\lvert}\widehat{\widetilde{q}}_{k}^{\,(k+1)}-{\widetilde{q}}_{k}^{\,(k+1)}\bigr{\rvert}}\lesssim\Bigl{(}\sqrt{{\lvert
q_{k}\rvert}^{2}+{\lvert q_{k+1}\rvert}^{2}+\cdots+{\lvert
q_{n}\rvert}^{2}}\Bigr{)}\mathit{u}.$ (128)

Thus,

$\displaystyle{\bigl{\lvert}\widehat{p}_{k-2}^{\,(k+1)}\widehat{\widetilde{q}}_{k}^{\,(k+1)*}-\widehat{p}_{k-2}^{\,(k+1)}{\widetilde{q}}_{k}^{\,(k+1)*}\bigr{\rvert}}\lesssim\Bigl{(}\sqrt{{\lvert\widehat{p}_{k-2}^{\,(k+1)}q_{k}^{*}\rvert}^{2}+{\lvert\widehat{p}_{k-2}^{\,(k+1)}q_{k+1}^{*}\rvert}^{2}+\cdots+{\lvert\widehat{p}_{k-2}^{\,(k+1)}q_{n}^{*}\rvert}^{2}}\Bigr{)}\mathit{u}$
$\displaystyle=\Bigl{(}\sqrt{{\lvert\widehat{b}_{k-2,k}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{b}_{k-2,k+1}^{\,(k+1)}\rvert}^{2}+\cdots+{\lvert\widehat{b}_{k-2,n}^{\,(k+1)}\rvert}^{2}}\Bigr{)}\mathit{u}$
$\displaystyle\lesssim{\lVert A\rVert}\mathit{u}.$ (129)

where the first equality follows by (123) and the second inequality by (124).
Combining (126) and (129),

$\displaystyle{\bigl{\lvert}-\widehat{p}_{k-2}^{\,(k+1)}\widehat{\widetilde{q}}_{k}^{\,(k+1)*}-\overline{b_{k,k-2}^{(k+1)}}\bigr{\rvert}}\lesssim{\lVert
A\rVert}\mathit{u},$ (130)

or, equivalently,

$\displaystyle{\bigl{\lvert}-\widehat{\widetilde{q}}_{k}^{\,(k+1)}\widehat{p}_{k-2}^{\,(k+1)*}-b_{k,k-2}^{(k+1)}\bigr{\rvert}}\lesssim{\lVert
A\rVert}\mathit{u}.$ (131)

Finally, since, by the induction hypothesis,

$\displaystyle{\bigl{\lvert}\widehat{\gamma}_{k-2}^{(k+1)}-b_{k-1,k-2}^{(k+1)}\bigr{\rvert}}\lesssim{\lVert
A\rVert}\mathit{u},$ (132)

we use (131) and (132) and apply Lemma 2.5 to (118) to find that
${\lvert\widehat{\gamma}_{k-2}^{(k)}-b_{k-1,k-2}^{(k)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$. In Line 8, we have

$\displaystyle\left[\begin{array}[]{c}\widehat{d}_{k-1}^{\;(k)}\\\
\widehat{\gamma}_{k-1}^{(k)}\end{array}\right]=\mathop{\mathit{fl}}\biggl{(}\widehat{Q}_{k}\left[\begin{array}[]{c}\widehat{d}_{k-1}^{\;(k+1)}\\\
\widehat{\gamma}_{k-1}^{(k+1)}\end{array}\right]\biggr{)}.$ (137)

By the induction hypothesis,
${\lvert\widehat{d}_{k-1}^{\;(k+1)}-b_{k-1,k-1}^{(k+1)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$ and
${\lvert\widehat{\gamma}_{k-1}^{(k+1)}-b_{k,k-1}^{(k+1)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$. Thus, another application of Lemma 2.5 shows that
${\lvert\widehat{d}_{k-1}^{\;(k)}-b_{k-1,k-1}^{(k)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$ and
${\lvert\widehat{\gamma}_{k-1}^{(k)}-b_{k,k-1}^{(k)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$. In Line 9, we have

$\displaystyle\left[\begin{array}[]{c}\widehat{\beta}_{k-1}^{\,(k)}\\\
\widehat{d}_{k}^{\;(k)}\end{array}\right]=\mathop{\mathit{fl}}\biggl{(}\widehat{Q}_{k}\left[\begin{array}[]{c}\widehat{\beta}_{k-1}^{\,(k+1)}\\\
\widehat{d}_{k}^{\;(k+1)}\end{array}\right]\biggr{)}.$ (142)

By the induction hypothesis,
${\lvert\widehat{\beta}_{k-1}^{\,(k+1)}-b_{k-1,k}^{(k+1)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$ and
${\lvert\widehat{d}_{k}^{\;(k+1)}-b_{k,k}^{(k+1)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$, so it follows from Lemma 2.5 that
${\lvert\widehat{\beta}_{k-1}^{\,(k)}-b_{k-1,k}^{(k)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$ and
${\lvert\widehat{d}_{k}^{\;(k)}-b_{k,k}^{(k)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$. In Line 10, we then have

$\displaystyle\left[\begin{array}[]{c}\widehat{p}_{k-1}^{\,(k)\dagger}\\\
\widehat{p}_{k}^{\,(k)}\end{array}\right]=\mathop{\mathit{fl}}\biggl{(}\widehat{Q}_{k}\left[\begin{array}[]{c}\widehat{p}_{k-1}^{\,(k+1)}\\\
\widehat{p}_{k}^{\,(k+1)}\end{array}\right]\biggr{)},$ (147)

where $\widehat{p}_{k-1}^{\,(k)\dagger}$ is a temporary value that will be
corrected later. By the induction hypothesis,
${\lvert\widehat{p}_{k-1}^{\,(k+1)}-p_{k-1}^{(k+1)}\rvert}\lesssim{\lVert
p\rVert}\mathit{u}$ and
${\lvert\widehat{p}_{k}^{\,(k+1)}-p_{k}^{(k+1)}\rvert}\lesssim{\lVert
p\rVert}\mathit{u}$, so it follows from Lemma 2.5 that
${\lvert\widehat{p}_{k-1}^{\,(k)\dagger}-p_{k-1}^{(k)}\rvert}\lesssim{\lVert
p\rVert}\mathit{u}$ and
${\lvert\widehat{p}_{k}^{\,(k)}-p_{k}^{(k)}\rvert}\lesssim{\lVert
p\rVert}\mathit{u}$. Define $\mathring{\beta}^{(k)}_{k-1}$ and
$\mathring{d}_{k}^{(k)}$ by the formula

$\displaystyle\left[\begin{array}[]{c}\mathring{\beta}_{k-1}^{(k)}\\\
\mathring{d}_{k}^{(k)}\end{array}\right]=Q_{k}\left[\begin{array}[]{c}\widehat{\beta}_{k-1}^{\,(k+1)}\\\
\widehat{d}_{k}^{\;(k+1)}\end{array}\right],$ (152)

and define $\mathring{p}_{k-1}^{(k)}$ and $\mathring{p}_{k}^{(k)}$ by

$\displaystyle\left[\begin{array}[]{c}\mathring{p}_{k-1}^{(k)}\\\
\mathring{p}_{k}^{(k)}\end{array}\right]=Q_{k}\left[\begin{array}[]{c}\widehat{p}_{k-1}^{\,(k+1)}\\\
\widehat{p}_{k}^{\,(k+1)}\end{array}\right].$ (157)

Clearly,

$\displaystyle\mathring{\beta}_{k-1}^{(k)}+\mathring{p}_{k-1}^{(k)}q_{k}^{*}=0,$
(158)

by the definition of $Q_{k}$ (see (115)). Also, by Lemma 2.5, we have that
${\lvert\mathring{\beta}_{k-1}^{(k)}-\widehat{\beta}_{k-1}^{\,(k)}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\widehat{\beta}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{d}_{k}^{\;(k+1)}\rvert}^{2}}\Bigr{)}\mathit{u}\leq{\lVert
A\rVert}\mathit{u}$ and
${\lvert\mathring{p}_{k-1}^{(k)}-\widehat{p}_{k-1}^{\,(k)\dagger}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\widehat{p}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}\rvert}^{2}}\Bigr{)}\mathit{u}\leq{\lVert
p\rVert}\mathit{u}$. In Line 12, we apply a correction to
$\widehat{p}_{k-1}^{\,(k)\dagger}$, so that

$\displaystyle\widehat{p}_{k-1}^{\,(k)}=\left\\{\begin{array}[]{ll}-\widehat{\beta}_{k-1}^{\,(k)}/q_{k}^{*}&\text{if
${\lvert\widehat{p}_{k-1}^{\,(k+1)}q_{k}^{*}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}q_{k}^{*}\rvert}^{2}>{\lvert\widehat{\beta}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{d}_{k}^{\;(k+1)}\rvert}^{2}$},\\\
\widehat{p}_{k-1}^{\,(k)\dagger}&\text{if
${\lvert\widehat{p}_{k-1}^{\,(k+1)}q_{k}^{*}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}q_{k}^{*}\rvert}^{2}\leq{\lvert\widehat{\beta}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{d}_{k}^{\;(k+1)}\rvert}^{2}$}.\end{array}\right.$
(161)

Thus, if
${\lvert\widehat{p}_{k-1}^{\,(k+1)}q_{k}^{*}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}q_{k}^{*}\rvert}^{2}>{\lvert\widehat{\beta}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{d}_{k}^{\;(k+1)}\rvert}^{2}$,
then

$\displaystyle{\bigl{\lvert}\widehat{p}_{k-1}^{\,(k)}-p_{k-1}^{(k)}\bigr{\rvert}}={\bigl{\lvert}-\widehat{\beta}_{k-1}^{\,(k)}/q_{k}^{*}-p_{k-1}^{(k)}\bigr{\rvert}}.$
(162)

We then observe that

$\displaystyle{\bigl{\lvert}-\mathring{\beta}_{k-1}^{(k)}/q_{k}^{*}-p_{k-1}^{(k)}\bigr{\rvert}}={\bigl{\lvert}\mathring{p}_{k-1}^{(k)}-p_{k-1}^{(k)}\bigr{\rvert}}\lesssim{\lVert
p\rVert}\mathit{u},$ (163)

and

$\displaystyle{\bigl{\lvert}-\mathring{\beta}_{k-1}^{(k)}/q_{k}^{*}+\widehat{\beta}_{k-1}^{\,(k)}/q_{k}^{*}\bigr{\rvert}}=\frac{1}{{\lvert
q_{k}^{*}\rvert}}{\bigl{\lvert}\widehat{\beta}_{k-1}^{\,(k)}-\mathring{\beta}_{k-1}^{(k)}\bigr{\rvert}}$
$\displaystyle\lesssim\frac{\sqrt{{\lvert\widehat{\beta}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{d}_{k}^{\;(k+1)}\rvert}^{2}}}{{\lvert
q_{k}^{*}\rvert}}\mathit{u}$
$\displaystyle\leq\Bigl{(}\sqrt{{\lvert\widehat{p}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}\rvert}^{2}}\Bigr{)}\mathit{u}$
$\displaystyle\leq{\lVert p\rVert}\mathit{u}.$ (164)

Finally, combining (162), (163), and (164), we find that
${\lvert\widehat{p}_{k-1}^{\,(k)}-p_{k-1}^{(k)}\rvert}\lesssim{\lVert
p\rVert}\mathit{u}$. Furthermore,

$\displaystyle{\bigl{\lvert}-\widehat{p}_{k-1}^{\,(k)}q_{k}^{*}-b_{k-1,k}^{(k)}\bigr{\rvert}}={\bigl{\lvert}\widehat{\beta}_{k-1}^{\,(k)}-b_{k-1,k}^{(k)}\bigr{\rvert}}\lesssim{\lVert
A\rVert}\mathit{u}.$ (165)

If, conversely,
${\lvert\widehat{p}_{k-1}^{\,(k+1)}q_{k}^{*}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}q_{k}^{*}\rvert}^{2}\leq{\lvert\widehat{\beta}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{d}_{k}^{\;(k+1)}\rvert}^{2}$,
then

$\displaystyle{\bigl{\lvert}\widehat{p}_{k-1}^{\,(k)}-p_{k-1}^{(k)}\bigr{\rvert}}={\bigl{\lvert}\widehat{p}_{k-1}^{\,(k)\dagger}-p_{k-1}^{(k)}\bigr{\rvert}}\lesssim{\lVert
p\rVert}\mathit{u}.$ (166)

Next, we observe that

$\displaystyle{\bigl{\lvert}-\widehat{p}_{k-1}^{\,(k)}q_{k}^{*}-b_{k-1,k}^{(k)}\bigr{\rvert}}={\bigl{\lvert}-\widehat{p}_{k-1}^{\,(k)\dagger}q_{k}^{*}-b_{k-1,k}^{(k)}\bigr{\rvert}}.$
(167)

Since

$\displaystyle{\bigl{\lvert}-\mathring{p}_{k-1}^{(k)}q_{k}^{*}-b_{k-1,k}^{(k)}\bigr{\rvert}}={\bigl{\lvert}\mathring{\beta}_{k-1}^{(k)}-b_{k-1,k}^{(k)}\bigr{\rvert}}\lesssim{\lVert
A\rVert}\mathit{u},$ (168)

and

$\displaystyle{\bigl{\lvert}-\mathring{p}_{k-1}^{(k)}q_{k}^{*}+\widehat{p}_{k-1}^{\,(k)\dagger}q_{k}^{*}\bigr{\rvert}}={\lvert
q_{k}^{*}\rvert}{\bigl{\lvert}-\mathring{p}_{k-1}^{(k)}+\widehat{p}_{k-1}^{\,(k)\dagger}\bigr{\rvert}}$
$\displaystyle\lesssim{\lvert
q_{k}^{*}\rvert}\Bigl{(}\sqrt{{\lvert\widehat{p}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}\rvert}^{2}}\Bigr{)}\mathit{u}$
$\displaystyle\leq\Bigl{(}\sqrt{{\lvert\widehat{\beta}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{d}_{k}^{\;(k+1)}\rvert}^{2}}\Bigr{)}\mathit{u}$
$\displaystyle\leq{\lVert A\rVert}\mathit{u},$ (169)

we combine (167), (168), and (169) to see that
${\lvert-\widehat{p}_{k-1}^{\,(k)}q_{k}^{*}-b_{k-1,k}^{(k)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$.

Now all that’s left is to show that
${\lvert-\widehat{p}_{k-1}^{\,(k)}q_{\ell}^{*}-b_{k-1,\ell}^{(k)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$ and
${\lvert-\widehat{p}_{k}^{\,(k)}q_{\ell}^{*}-b_{k,\ell}^{(k)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$, for all $\ell=k+1,k+2,\ldots,n$. By the induction
hypothesis,

$\displaystyle{\bigl{\lvert}-\widehat{p}_{k-1}^{\,(k+1)}q_{\ell}^{*}-b_{k-1,\ell}^{(k+1)}\bigr{\rvert}}\lesssim{\lVert
A\rVert}\mathit{u}$ (170)

and

$\displaystyle{\bigl{\lvert}-\widehat{p}_{k}^{\,(k+1)}q_{\ell}^{*}-b_{k,\ell}^{(k+1)}\bigr{\rvert}}\lesssim{\lVert
A\rVert}\mathit{u},$ (171)

for all $\ell=k+1,k+2,\ldots,n$. Multiplying (157) by $q_{\ell}^{*}$, we have

$\displaystyle\left[\begin{array}[]{c}\mathring{p}_{k-1}^{(k)}q_{\ell}^{*}\\\
\mathring{p}_{k}^{(k)}q_{\ell}^{*}\end{array}\right]=Q_{k}\left[\begin{array}[]{c}\widehat{p}_{k-1}^{\,(k+1)}q_{\ell}^{*}\\\
\widehat{p}_{k}^{\,(k+1)}q_{\ell}^{*}\end{array}\right],$ (176)

which, combined with (170) and (171), means that

$\displaystyle{\bigl{\lvert}-\mathring{p}_{k-1}^{(k)}q_{\ell}^{*}-b_{k-1,\ell}^{(k)}\bigr{\rvert}}\lesssim{\lVert
A\rVert}\mathit{u}$ (177)

and

$\displaystyle{\bigl{\lvert}-\mathring{p}_{k}^{(k)}q_{\ell}^{*}-b_{k,\ell}^{(k)}\bigr{\rvert}}\lesssim{\lVert
A\rVert}\mathit{u},$ (178)

for all $\ell=k+1,k+2,\ldots,n$. It is not difficult to show (see (164)) that

$\displaystyle{\lvert\mathring{p}_{k-1}^{(k)}-\widehat{p}_{k-1}^{\,(k)}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\widehat{p}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}\rvert}^{2}}\Bigr{)}\mathit{u}$
(179)

and

$\displaystyle{\lvert\mathring{p}_{k}^{(k)}-\widehat{p}_{k}^{\,(k)}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\widehat{p}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}\rvert}^{2}}\Bigr{)}\mathit{u},$
(180)

from which it follows that

$\displaystyle{\lvert\mathring{p}_{k-1}^{(k)}q_{\ell}^{*}-\widehat{p}_{k-1}^{\,(k)}q_{\ell}^{*}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\widehat{p}_{k-1}^{\,(k+1)}q_{\ell}^{*}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}q_{\ell}^{*}\rvert}^{2}}\Bigr{)}\mathit{u}\lesssim{\lVert
A\rVert}\mathit{u}$ (181)

and

$\displaystyle{\lvert\mathring{p}_{k}^{(k)}q_{\ell}^{*}-\widehat{p}_{k}^{\,(k)}q_{\ell}^{*}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\widehat{p}_{k-1}^{\,(k+1)}q_{\ell}^{*}\rvert}^{2}+{\lvert\widehat{p}_{k}^{\,(k+1)}q_{\ell}^{*}\rvert}^{2}}\Bigr{)}\mathit{u}\lesssim{\lVert
A\rVert}\mathit{u},$ (182)

for all $\ell=k+1,k+2,\ldots,n$, where the second inequality follows from
(170) and (171). Finally, combining (177), (178), (181), and (182), we find
that
${\lvert-\widehat{p}_{k-1}^{\,(k)}q_{\ell}^{*}-b_{k-1,\ell}^{(k)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$ and
${\lvert-\widehat{p}_{k}^{\,(k)}q_{\ell}^{*}-b_{k,\ell}^{(k)}\rvert}\lesssim{\lVert
A\rVert}\mathit{u}$, for all $\ell=k+1,k+2,\ldots,n$, and we are done.

$\blacksquare$

The following lemma bounds the forward error of Algorithm 2 (the rotation back
to Hessenberg form).

######  Lemma 4.2.

Suppose that $B\in\mathbbm{C}^{n\times n}$ and $p,q\in\mathbbm{C}^{n}$.
Suppose further that $B+pq^{*}$ is lower triangular, and let $d$ and $\gamma$
denote the diagonal and subdiagonal of $B$, respectively. Suppose that
$Q_{2},Q_{3},\ldots,Q_{n}\in\mathrm{SU}(2)$, and suppose that Algorithm 2 is
carried out in floating point arithmetic, using $d$, $\gamma$, $p$, $q$, and
$Q_{2},Q_{3},\ldots,Q_{n}$ as inputs. Suppose finally that
$\widehat{\underline{d}}$, $\widehat{\underline{\beta}}$, and
$\widehat{\underline{q}}$ are the outputs generated by Algorithm 2, and define
the upper triangular part of the matrix
$\widehat{\underline{A}}\in\mathbbm{C}^{n\times n}$ by the formula

$\displaystyle{\widehat{\underline{a}}}_{i,j}=\left\\{\begin{array}[]{ll}-p_{i}\widehat{\underline{q}}_{j}^{\,*}&\text{if
$j>i+1$},\\\ \underline{\widehat{\beta}}_{i}&\text{if $j=i+1$},\\\
\underline{\widehat{d}}_{i}&\text{if $j=i$},\end{array}\right.$ (186)

where $\widehat{\underline{a}}_{i,j}$ denotes the $(i,j)$-th entry of
$\widehat{\underline{A}}$. Let $U_{k}\in\mathbbm{C}^{n\times n}$,
$k=2,3,\ldots,n$, denote the matrices that rotate the $(k-1,k)$-plane by
$Q_{k}$. Define $U\in\mathbbm{C}^{n\times n}$ by the formula
$U=U_{2}U_{3}\cdots U_{n}$, and let $\underline{A}=BU^{*}$ and
$\underline{q}=Uq$. Then

$\displaystyle{\lVert\widehat{\underline{A}}-\underline{A}\rVert}_{T}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$ (187)

and

$\displaystyle{\lVert\widehat{\underline{q}}-\underline{q}\rVert}\lesssim{\lVert
q\rVert}\mathit{u},$ (188)

where ${\lVert\cdot\rVert}_{T}$ denotes the square root of the sum of squares
of the entries in the upper triangular part of its argument and
${\lVert\cdot\rVert}_{H}$ denotes the square root of the sum of squares of the
upper Hessenberg part (see Definition 2.1).

Proof. Suppose that $\widehat{d}^{\,(k)}$, $\widehat{\beta}^{\,(k)}$, and
$\widehat{q}^{\,(k)}$ denote the computed vectors in Algorithm 2 after
rotations in the positions $(n-1,n),(n-2,n-1),\ldots,(k-1,k)$. Suppose further
that the upper triangular part of the matrix
$\widehat{A}^{(k)}\in\mathbbm{C}^{n\times n}$ is defined by the formula
$\displaystyle{\widehat{a}^{\,(k)}}_{i,j}=\left\\{\begin{array}[]{ll}-p_{i}\widehat{q}_{j}^{\,(k)*}&\text{if
$j>i+1$ or if $j=i+1$ and $j<k$},\\\ \widehat{\beta}^{\,(k)}_{i}&\text{if
$j=i+1$ and $j\geq k$},\\\ \widehat{d}^{\;(k)}_{i}&\text{if
$j=i$},\end{array}\right.$ (192) where $\widehat{a}_{i,j}^{\,(k)}$ denotes the
$(i,j)$-th entry of $\widehat{A}^{(k)}$. Clearly,
$\widehat{\underline{d}}=\widehat{d}^{\,(2)}$,
$\widehat{\underline{\beta}}=\widehat{\beta}^{(2)}$,
$\widehat{\underline{q}}=\widehat{q}^{\,(2)}$, and
$\widehat{\underline{A}}=\widehat{A}^{(2)}$. Let
$A^{(k)}=BU_{n}^{*}U_{n-1}^{*}\cdots U_{k}^{*}$ and
$q^{(k)}=U_{k}U_{k+1}\cdots U_{n}q$. We will prove that
${\lVert\widehat{A}^{(k)}-A^{(k)}\rVert}_{T}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$ and
${\lVert\widehat{q}^{\,(k)}-q^{(k)}\rVert}\lesssim{\lVert q\rVert}u$, for each
$k=n,n-1,\ldots,2$. Define $\widehat{d}^{\,(n+1)}=d$,
$\widehat{q}^{\,(n+1)}=q$, $q^{(n+1)}=q$, and $A^{(n+1)}=B$. Obviously,
$\widehat{A}^{(n+1)}=A^{(n+1)}$ and $\widehat{q}^{\,(n+1)}=q^{(n+1)}$, so the
above statement is true for $k=n+1$. We will prove it for the cases
$k=n,n-1,\ldots,2$ by induction. In Line 2, we have
$\displaystyle\left[\begin{array}[]{c}\widehat{d}_{k-1}^{\;(k)}\\\
\widehat{\beta}_{k-1}^{\,(k)}\end{array}\right]=\mathop{\mathit{fl}}\biggl{(}\widehat{\overline{Q_{k}}}\left[\begin{array}[]{c}\widehat{d}_{k-1}^{\;(k+1)}\\\
-p_{k-1}\widehat{q}_{k}^{\,(k+1)*}\end{array}\right]\biggr{)}.$ (197) By the
induction hypothesis,
${\lvert\widehat{d}_{k-1}^{\;(k+1)}-a_{k-1,k-1}^{(k+1)}\rvert}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$ and ${\lvert-
p_{k-1}\widehat{q}_{k}^{\,(k+1)*}-a_{k-1,k}^{(k+1)}\rvert}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$. Since, by Lemma 2.3, ${\lVert
A^{(k+1)}\rVert}_{T}\leq{\lVert B\rVert}_{H}$, an application of Lemma 2.5
gives us
${\lvert\widehat{d}_{k-1}^{\;(k)}-a_{k-1,k-1}^{(k)}\rvert}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$ and ${\lvert-
p_{k-1}\widehat{q}_{k}^{\,(k)*}-a_{k-1,k}^{(k)}\rvert}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$. In Line 3, we have
$\displaystyle\widehat{d}_{k}^{\;(k)}=\mathop{\mathit{fl}}\biggl{(}\widehat{\overline{Q_{k}}}\left[\begin{array}[]{c}\gamma_{k-1}\\\
\widehat{d}_{k}^{\;(k+1)}\end{array}\right]\biggr{)}_{2}.$ (200) We first
observe that
$\displaystyle\gamma_{k-1}=b_{k,k-1}=(BU_{n}^{*}U_{n-1}^{*}\cdots
U_{k+1}^{*})_{k,k-1}=a^{(k+1)}_{k,k-1},$ (201) since right multiplication by
$U_{j}^{*}$ only affects columns $j$ and $j-1$. Since, by the induction
hypothesis,
${\lvert\widehat{d}_{k}^{\;(k+1)}-a_{k,k}^{(k+1)}\rvert}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$, an application of Lemma 2.5 together with the
inequality ${\lVert A^{(k+1)}\rVert}_{T}\leq{\lVert B\rVert}_{H}$ gives us
${\lvert\widehat{d}_{k}^{\;(k)}-a_{k,k}^{(k)}\rvert}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$. In Line 4,
$\displaystyle\left[\begin{array}[]{c}\widehat{q}_{k-1}^{\,(k)}\\\
\widehat{q}_{k}^{\,(k)}\end{array}\right]=\mathop{\mathit{fl}}\biggl{(}\widehat{Q}_{k}\left[\begin{array}[]{c}\widehat{q}_{k-1}^{\,(k+1)}\\\
\widehat{q}_{k}^{\,(k+1)}\end{array}\right]\biggr{)}.$ (206) By the induction
hypothesis,
${\lvert\widehat{q}_{k-1}^{\,(k+1)}-q_{k-1}^{(k+1)}\rvert}\lesssim{\lVert
q\rVert}\mathit{u}$ and
${\lvert\widehat{q}_{k}^{\,(k+1)}-q_{k}^{(k+1)}\rvert}\lesssim{\lVert
q\rVert}\mathit{u}$, so it follows from Lemma 2.5 that
${\lvert\widehat{q}_{k-1}^{\,(k)}-q_{k-1}^{(k)}\rvert}\lesssim{\lVert
q\rVert}\mathit{u}$ and
${\lvert\widehat{q}_{k}^{\,(k)}-q_{k}^{(k)}\rvert}\lesssim{\lVert
q\rVert}\mathit{u}$. All that’s left now is to prove that ${\lvert-
p_{\ell}\widehat{q}_{k-1}^{\,(k)*}-a_{\ell,k-1}^{(k)}\rvert}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$ and ${\lvert-
p_{\ell}\widehat{q}_{k}^{\,(k)*}-a_{\ell,k}^{(k)}\rvert}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$, for all for all $\ell=1,2,\ldots,k-2$. By the
induction hypothesis,  $\displaystyle{\lvert-
p_{\ell}\widehat{q}_{k-1}^{\,(k+1)*}-a_{\ell,k-1}^{(k+1)}\rvert}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$ (207) and  $\displaystyle{\lvert-
p_{\ell}\widehat{q}_{k}^{\,(k+1)*}-a_{\ell,k}^{(k+1)}\rvert}\lesssim{\lVert
B\rVert}_{H}\mathit{u},$ (208) for all $\ell=1,2,\ldots,k-2$. Define
$\mathring{q}_{k-1}^{(k)}$ and $\mathring{q}_{k}^{(k)}$ by
$\displaystyle\left[\begin{array}[]{c}\mathring{q}_{k-1}^{(k)}\\\
\mathring{q}_{k}^{(k)}\end{array}\right]=Q_{k}\left[\begin{array}[]{c}\widehat{q}_{k-1}^{\,(k+1)}\\\
\widehat{q}_{k}^{\,(k+1)}\end{array}\right].$ (213) Multiplying (213) by
$p_{\ell}^{*}$, we have
$\displaystyle\left[\begin{array}[]{c}\mathring{q}_{k-1}^{(k)}p_{\ell}^{*}\\\
\mathring{q}_{k}^{(k)}p_{\ell}^{*}\end{array}\right]=Q_{k}\left[\begin{array}[]{c}\widehat{q}_{k-1}^{\,(k+1)}p_{\ell}^{*}\\\
\widehat{q}_{k}^{\,(k+1)}p_{\ell}^{*}\end{array}\right],$ (218) which,
combined with (207) and (208) and the fact that ${\lVert
A^{(k+1)}\rVert}_{T}\leq{\lVert B\rVert}_{H}$, means that
$\displaystyle{\bigl{\lvert}-p_{\ell}\mathring{q}_{k-1}^{(k)*}-a_{\ell,k-1}^{(k)}\bigr{\rvert}}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$ (219) and
$\displaystyle{\bigl{\lvert}-p_{\ell}\mathring{q}_{k}^{(k)*}-a_{\ell,k}^{(k)}\bigr{\rvert}}\lesssim{\lVert
B\rVert}_{H}\mathit{u},$ (220) for all $\ell=1,2,\ldots,k-2$. By Lemma 2.5,
$\displaystyle{\lvert\mathring{q}_{k-1}^{(k)}-\widehat{q}_{k-1}^{\,(k)}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\widehat{q}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{q}_{k}^{\,(k+1)}\rvert}^{2}}\Bigr{)}\mathit{u}$
(221) and
$\displaystyle{\lvert\mathring{q}_{k}^{(k)}-\widehat{q}_{k}^{\,(k)}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\widehat{q}_{k-1}^{\,(k+1)}\rvert}^{2}+{\lvert\widehat{q}_{k}^{\,(k+1)}\rvert}^{2}}\Bigr{)}\mathit{u},$
(222) from which it follows that
$\displaystyle{\lvert\mathring{q}_{k-1}^{(k)}p_{\ell}^{*}-\widehat{q}_{k-1}^{\,(k)}p_{\ell}^{*}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\widehat{q}_{k-1}^{\,(k+1)}p_{\ell}^{*}\rvert}^{2}+{\lvert\widehat{q}_{k}^{\,(k+1)}p_{\ell}^{*}\rvert}^{2}}\Bigr{)}\mathit{u}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$ (223) and
$\displaystyle{\lvert\mathring{q}_{k}^{(k)}p_{\ell}^{*}-\widehat{q}_{k}^{\,(k)}p_{\ell}^{*}\rvert}\lesssim\Bigl{(}\sqrt{{\lvert\widehat{q}_{k-1}^{\,(k+1)}p_{\ell}^{*}\rvert}^{2}+{\lvert\widehat{q}_{k}^{\,(k+1)}p_{\ell}^{*}\rvert}^{2}}\Bigr{)}\mathit{u}\lesssim{\lVert
B\rVert}_{H}\mathit{u},$ (224) for all $\ell=1,2,\ldots,k-2$, where the second
inequality follows from (207) and (208) and the inequality ${\lVert
A^{(k+1)}\rVert}_{T}\leq{\lVert B\rVert}_{H}$. Finally, combining (219),
(220), (223), and (224), we find that ${\lvert-
p_{\ell}\widehat{q}_{k-1}^{\,(k)*}-a_{\ell,k-1}^{(k)}\rvert}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$ and ${\lvert-
p_{\ell}\widehat{q}_{k}^{\,(k)*}-a_{\ell,k}^{(k)}\rvert}\lesssim{\lVert
B\rVert}_{H}\mathit{u}$, for all for all $\ell=1,2,\ldots,k-2$, and we are
done. $\blacksquare$

The following theorem bounds the forward errors of a full sweep of our $QR$
algorithm, and is the principal result of this subsection.

######  Theorem 4.3.

Suppose that $A\in\mathbbm{C}^{n\times n}$ is a Hermitian matrix, that
$p,q\in\mathbbm{C}^{n}$, and that $A+pq^{*}$ is lower Hessenberg. Let $d$ and
$\beta$ denote the diagonal and superdiagonal of $A$, respectively. Suppose
that Algorithm 1 is carried out in floating point arithmetic, and let
$Q_{2},Q_{3},\ldots,Q_{n}\in\mathrm{SU}(2)$ be the unitary matrices generated
by an exact step of Line 4 of Algorithm 1 applied to the computed vectors at
that step. Let $U_{k}\in\mathbbm{C}^{n\times n}$, $k=2,3,\ldots,n$, denote the
matrices that rotate the $(k-1,k)$-plane by $Q_{k}$, and define
$U\in\mathbbm{C}^{n\times n}$ by the formula $U=U_{2}U_{3}\cdots U_{n}$.
Suppose that Algorithm 2 is then carried out in floating point arithmetic,
using the outputs of Algorithm 1 as inputs. Suppose finally that
$\widehat{\underline{p}}$ is an output of Algorithm 1 and
$\widehat{\underline{q}}$, $\widehat{\underline{d}}$, and
$\widehat{\underline{\beta}}$ are all outputs of Algorithm 2, and define the
matrix $\widehat{\underline{A}}$ by the formula

$\displaystyle\widehat{\underline{a}}_{i,j}=\left\\{\begin{array}[]{ll}-\widehat{\underline{p}}_{i}\widehat{\underline{q}}_{j}^{*}&\text{if
$j>i+1$}\\\ \widehat{\underline{\beta}}_{i}&\text{if $j=i+1$}\\\
\widehat{\underline{d}}_{i}&\text{if $j=i$}\\\
\overline{(\widehat{\underline{\beta}}_{j})}&\text{if $j=i-1$}\\\
-\widehat{\underline{q}}_{j}\widehat{\underline{p}}_{i}^{*}&\text{if
$j<i-1$}\end{array}\right.$ (230)

where $\widehat{\underline{a}}_{i,j}$ denotes the $(i,j)$-th entry of
$\widehat{\underline{A}}$. Let $\underline{A}=UAU^{*}$, $\underline{p}=Up$,
and $\underline{q}=Uq$. Then

$\displaystyle{\lVert\widehat{\underline{A}}-\underline{A}\rVert}\lesssim{\lVert
A\rVert}\mathit{u},$ (231)
$\displaystyle{\lVert\widehat{\underline{p}}-\underline{p}\rVert}\lesssim{\lVert
p\rVert}\mathit{u},$ (232)

and

$\displaystyle{\lVert\widehat{\underline{q}}-\underline{q}\rVert}\lesssim{\lVert
q\rVert}\mathit{u}.$ (233)

Proof. Suppose that $\widehat{B}$ (defined by (76)),
$\widehat{\underline{p}}$, and
$\widehat{Q}_{2},\widehat{Q}_{3},\ldots,\widehat{Q}_{n}\in\mathbbm{C}^{n\times
n}$ are outputs of Algorithm 1. Let $B=UA$ and $\underline{p}=Up$. By Lemma
4.1,  $\displaystyle{\lVert{\widehat{B}}-B\rVert}_{H}\lesssim{\lVert
A\rVert}\mathit{u}$ (234) and
$\displaystyle{\lVert\underline{\widehat{p}}-\underline{p}\rVert}\lesssim{\lVert
p\rVert}\mathit{u},$ (235) where ${\lVert\cdot\rVert}_{H}$ denotes the square
root of the sum of squares of the entries in the upper Hessenberg part of its
argument (see Definition 2.1). Now suppose that $\widehat{B}$,
$\underline{\widehat{p}}$, $q$, and
$\widehat{Q}_{2},\widehat{Q}_{3},\ldots,\widehat{Q}_{n}\in\mathbbm{C}^{n\times
n}$ are used as inputs to Algorithm 2. Let
$\widehat{U}_{k}\in\mathbbm{C}^{n\times n}$, $k=2,3,\ldots,n$, denote the
matrices that rotate the $(k-1,k)$-plane by $\widehat{Q}_{k}$, and define
$\widehat{U}\in\mathbbm{C}^{n\times n}$ by the formula
$\widehat{U}=\widehat{U}_{2}\widehat{U}_{3}\cdots\widehat{U}_{n}$. Let
$\undertilde{A}=\widehat{B}\widehat{U}^{*}$ and $\undertilde{q}=\widehat{U}q$,
and observe that the upper triangular part of $\undertilde{A}$ is well-defined
due to Lemma 2.3. By Lemma 4.2 we have that
$\displaystyle{\lVert{\widehat{\underline{A}}}-\undertilde{A}\rVert}_{T}\lesssim{\lVert\widehat{B}\rVert}_{H}\mathit{u}$
(236) and
$\displaystyle{\lVert\underline{\widehat{q}}-\undertilde{q}\rVert}\lesssim{\lVert
q\rVert}\mathit{u},$ (237) where ${\lVert\cdot\rVert}_{T}$ denotes the square
root of the sum of squares of the entries in the upper triangular part of its
argument and ${\lVert\cdot\rVert}_{H}$ denotes the square root of the sum of
squares of the upper Hessenberg part (see Definition 2.1). Let
$\underline{A}=BU^{*}=UAU^{*}$ and let $\underline{q}=Uq$. We observe that
$\displaystyle{\lVert\undertilde{A}-\underline{A}\rVert}_{T}={\lVert\widehat{B}\widehat{U}^{*}-BU^{*}\rVert}_{T}$
$\displaystyle\leq{\lVert\widehat{B}\widehat{U}^{*}-B\widehat{U}^{*}\rVert}_{T}+{\lVert
B\widehat{U}^{*}-BU^{*}\rVert}_{T}$
$\displaystyle={\lVert(\widehat{B}-B)\widehat{U}^{*}\rVert}_{T}+{\lVert
B(\widehat{U}^{*}-U^{*})\rVert}_{T}$ $\displaystyle\lesssim{\lVert
A\rVert}\mathit{u},$ (238) where the last inequality follows from (234) and
the fact that ${\lVert\widehat{U}-U\rVert}\lesssim\mathit{u}$. Since, clearly,
${\lVert\widehat{B}\rVert}_{H}\mathit{u}\lesssim{\lVert A\rVert}\mathit{u}$,
we combine (236) and (238) to get
$\displaystyle{\lVert\widehat{\underline{A}}-\underline{A}\rVert}_{T}\lesssim{\lVert
A\rVert}\mathit{u}.$ (239) Now we observe that, since both
$\underline{A}=UAU^{*}$ and $\underline{\widehat{A}}$ are Hermitian,
$\displaystyle{\lVert\widehat{\underline{A}}-\underline{A}\rVert}\lesssim{\lVert
A\rVert}\mathit{u}.$ (240) Next, we observe that
$\displaystyle{\lVert\undertilde{q}-\underline{q}\rVert}={\lVert\widehat{U}q-Uq\rVert}$
$\displaystyle={\lVert(\widehat{U}-U)q\rVert}$ $\displaystyle\lesssim{\lVert
q\rVert}\mathit{u},$ (241) so, combining (237) and (241), we have
$\displaystyle{\lVert\widehat{\underline{q}}-\underline{q}\rVert}\lesssim{\lVert
q\rVert}\mathit{u},$ (242) and we are done. $\blacksquare$

###  4.2 Backward Error Analysis of the $QR$ Algorithms

Suppose that $A$ is Hermitian and $A+pq^{*}$ is lower Hessenberg. In this
section, we prove in Theorems 4.6 and 4.7 that the backward errors in $A$,
$p$, and $q$ of both our explicit unshifted $QR$ algorithm (see Algorithm 3)
and explicit shifted $QR$ algorithm (see Algorithm 4) are proportional to
${\lVert A\rVert}\mathit{u}$, ${\lVert p\rVert}\mathit{u}$, and ${\lVert
q\rVert}\mathit{u}$, respectively.

The following lemma states that a single iteration of our QR algorithm is
componentwise backward stable.

######  Lemma 4.4.

Suppose that $A\in\mathbbm{C}^{n\times n}$ is a Hermitian matrix, that
$p,q\in\mathbbm{C}^{n}$, and that $A+pq^{*}$ is lower Hessenberg. Let $d$ and
$\beta$ denote the diagonal and superdiagonal of $A$, respectively. Suppose
that a single iteration of our QR algorithm (Algorithm 1 followed by Algorithm
2) is carried out in floating point arithmetic, and let
$\widehat{\underline{p}}$, $\widehat{\underline{q}}$,
$\widehat{\underline{d}}$, and $\widehat{\underline{\beta}}$ denote the
outputs of the algorithm. Define the matrix $\widehat{\underline{A}}$ by the
formula (230). Then there exists a unitary matrix $U\in\mathbbm{C}^{n\times
n}$, a matrix $\delta A\in\mathbbm{C}^{n\times n}$, and vectors $\delta
p,\delta q\in\mathbbm{C}^{n}$, such that

$\displaystyle\widehat{\underline{A}}=U(A+\delta A)U^{*},$ (243)
$\displaystyle\widehat{\underline{p}}=U(p+\delta p),$ (244)

and

$\displaystyle\widehat{\underline{q}}=U(q+\delta q),$ (245)

where ${\lVert\delta A\rVert}\lesssim{\lVert A\rVert}\mathit{u}$,
${\lVert\delta p\rVert}\lesssim{\lVert p\rVert}\mathit{u}$, and ${\lVert\delta
q\rVert}\lesssim{\lVert q\rVert}\mathit{u}$.

Proof. Let $U\in\mathbbm{C}^{n\times n}$ be the unitary matrix defined in the
statement of Theorem 4.3, and let $\underline{A}=UAU^{*}$, $\underline{p}=Up$,
and $\underline{q}=Uq$. By Theorem 4.3,

$\displaystyle\widehat{\underline{A}}=\underline{A}+\underline{\delta A},$
(246)

where ${\lVert\underline{\delta A}\rVert}\lesssim{\lVert A\rVert}\mathit{u}$.
Thus,

$\displaystyle\widehat{\underline{A}}=UAU^{*}+\underline{\delta A}=U(A+\delta
A)U^{*},$ (247)

where $\delta A=U^{*}\underline{\delta A}U$. Since $U$ is unitary, clearly
${\lVert\delta A\rVert}\lesssim{\lVert A\rVert}\mathit{u}$. Likewise, by
Theorem 4.3,

$\displaystyle\widehat{\underline{p}}=\underline{p}+\underline{\delta p},$
(248)

where ${\lVert\underline{\delta p}\rVert}\lesssim{\lVert p\rVert}\mathit{u}$,
so

$\displaystyle\widehat{\underline{p}}=U(p+\delta p),$ (249)

where $\delta p=U^{*}\underline{\delta p}$ and ${\lVert\delta
p\rVert}\lesssim{\lVert p\rVert}\mathit{u}$. Similarly,

$\displaystyle\widehat{\underline{q}}=U(q+\delta q),$ (250)

where ${\lVert\delta q\rVert}\lesssim{\lVert q\rVert}\mathit{u}$.

$\blacksquare$

The following lemma states that repeated iterations of our QR algorithm are
componentwise backward stable.

######  Lemma 4.5.

Suppose that $A\in\mathbbm{C}^{n\times n}$ is a Hermitian matrix, that
$p,q\in\mathbbm{C}^{n}$, and that $A+pq^{*}$ is lower Hessenberg. Let $d$ and
$\beta$ denote the diagonal and superdiagonal of $A$, respectively. Suppose
that $k$ iterations of our QR algorithm (Algorithm 1 followed by Algorithm 2)
are carried out in floating point arithmetic, and let $\widehat{p}^{\,(k)}$,
$\widehat{q}^{\,(k)}$, $\widehat{d}^{\,(k)}$, and $\widehat{\beta}^{(k)}$
denote the outputs of the algorithm. Define the matrix $\widehat{A}^{(k)}$ by
the formula (230), making the obvious substitutions. Then there exists a
unitary matrix $U\in\mathbbm{C}^{n\times n}$, a matrix $\delta
A\in\mathbbm{C}^{n\times n}$, and vectors $\delta p,\delta
q\in\mathbbm{C}^{n}$, such that

$\displaystyle\widehat{A}^{(k)}=U(A+\delta A)U^{*},$ (251)
$\displaystyle\widehat{p}^{\,(k)}=U(p+\delta p),$ (252)

and

$\displaystyle\widehat{q}^{\,(k)}=U(q+\delta q),$ (253)

where ${\lVert\delta A\rVert}\lesssim{\lVert A\rVert}\mathit{u}$,
${\lVert\delta p\rVert}\lesssim{\lVert p\rVert}\mathit{u}$, and ${\lVert\delta
q\rVert}\lesssim{\lVert q\rVert}\mathit{u}$.

Proof. We will prove this statement only for the matrix $\widehat{A}^{(k)}$,
since the proofs for $\widehat{p}^{\,(k)}$ and $\widehat{q}^{\,(k)}$ are
essentially identical. By repeated application of Lemma 4.4, we know that
there exist unitary matrices $U^{(1)},U^{(2)},\ldots,U^{(k)}$ and matrices
$\delta
A^{(0)},\delta\widehat{A}^{(1)},\delta\widehat{A}^{(2)},\ldots,\delta\widehat{A}^{(k-1)}$
such that
$\displaystyle\widehat{A}^{(k)}=U^{(k)}(\widehat{A}^{(k-1)}+\delta\widehat{A}^{(k-1)})U^{(k)*},$
(254)
$\displaystyle\widehat{A}^{(k-1)}=U^{(k-1)}(\widehat{A}^{(k-2)}+\delta\widehat{A}^{(k-2)})U^{(k-1)*},$
(255) $\displaystyle\hskip 70.0001pt\vdots$
$\displaystyle\widehat{A}^{(2)}=U^{(2)}(\widehat{A}^{(1)}+\delta\widehat{A}^{(1)})U^{(2)*},$
(256) $\displaystyle\widehat{A}^{(1)}=U^{(1)}(A+\delta A^{(0)})U^{(1)*},$
(257) where ${\lVert\delta A^{(0)}\rVert}\lesssim{\lVert A\rVert}\mathit{u}$
and ${\lVert\delta\widehat{A}^{(\ell)}\rVert}\lesssim{\lVert
A\rVert}\mathit{u}$, for $\ell=1,2,\ldots,k-1$. Combining (254)–(257) and
expanding, we find that
$\displaystyle\widehat{A}^{(k)}=U^{(k)}U^{(k-1)}\cdots
U^{(1)}AU^{(1)*}U^{(2)*}\cdots
U^{(k)*}+U^{(k)}\delta\widehat{A}^{(k-1)}U^{(k)*}$
$\displaystyle+U^{(k)}U^{(k-1)}\delta\widehat{A}^{(k-2)}U^{(k-1)*}U^{(k)*}+\cdots$
$\displaystyle+U^{(k)}U^{(k-1)}\cdots U^{(1)}\delta A^{(0)}U^{(1)*}\cdots
U^{(k-1)*}U^{(k)*}.$ (258) Letting $U=U^{(k)}U^{(k-1)}\cdots U^{(1)}$, this
becomes
$\displaystyle\widehat{A}^{(k)}=UAU^{*}+U^{(k)}\delta\widehat{A}^{(k-1)}U^{(k)*}$
$\displaystyle+U^{(k)}U^{(k-1)}\delta\widehat{A}^{(k-2)}U^{(k-1)*}U^{(k)*}+\cdots$
$\displaystyle+U^{(k)}U^{(k-1)}\cdots U^{(1)}\delta A^{(0)}U^{(1)*}\cdots
U^{(k-1)*}U^{(k)*}.$ (259) Suppose now that the matrix $\delta A$ is defined
by  $\displaystyle\delta
A=U^{*}\bigl{(}U^{(k)}\delta\widehat{A}^{(k-1)}U^{(k)*}+U^{(k)}U^{(k-1)}\delta\widehat{A}^{(k-2)}U^{(k-1)*}U^{(k)*}+\cdots$
$\displaystyle+U^{(k)}U^{(k-1)}\cdots U^{(1)}\delta A^{(0)}U^{(1)*}\cdots
U^{(k-1)*}U^{(k)*}\bigr{)}U.$ (260) Clearly, ${\lVert\delta
A\rVert}\lesssim{\lVert A\rVert}\mathit{u}$. Combining (259) and (260), we
have  $\displaystyle\widehat{A}^{(k)}=U(A+\delta A)U^{*},$ (261) and we are
done. $\blacksquare$

The following theorem states that our explicit unshifted QR algorithm is
componentwise backward stable.

######  Theorem 4.6 (Explicit unshifted QR).

Suppose that $A\in\mathbbm{C}^{n\times n}$ is a Hermitian matrix, that
$p,q\in\mathbbm{C}^{n}$, and that $A+pq^{*}$ is lower Hessenberg. Let $d$ and
$\beta$ denote the diagonal and superdiagonal of $A$, respectively. Suppose
that Algorithm 3 is carried out in floating point arithmetic with
$\epsilon\lesssim{\lVert A\rVert}\mathit{u}$, and let
$\widehat{\lambda}_{1},\widehat{\lambda}_{2},\ldots,\widehat{\lambda}_{n}$
denote the outputs. Then there exist a matrix $\delta A\in\mathbbm{C}^{n\times
n}$ and vectors $\delta p,\delta q\in\mathbbm{C}^{n}$ such that
$\widehat{\lambda}_{1},\widehat{\lambda}_{2},\ldots,\widehat{\lambda}_{n}$ are
the exact eigenvalues of the matrix

$\displaystyle(A+\delta A)+(p+\delta p)(q+\delta q)^{*},$ (262)

where ${\lVert\delta A\rVert}\lesssim{\lVert A\rVert}\mathit{u}$,
${\lVert\delta p\rVert}\lesssim{\lVert p\rVert}\mathit{u}$, and ${\lVert\delta
q\rVert}\lesssim{\lVert q\rVert}\mathit{u}$.

Proof. Suppose that we carry out QR iterations until the entry in the $(1,2)$
position is less than $\epsilon$ in absolute value. Let $\widehat{d}^{\,(1)}$,
$\widehat{\beta}^{(1)\dagger}$, $\widehat{p}^{\,(1)}$, and
$\widehat{q}^{\,(1)}$ denote the resulting outputs, and let
$\widehat{A}^{(1)\dagger}$ be the resulting matrix, defined by formula (230)
(making the obvious substitutions). By Lemma 4.5, there exist a unitary matrix
$U^{(1)}\in\mathbbm{C}^{n\times n}$, a matrix $\delta
A^{(0)\dagger}\in\mathbbm{C}^{n\times n}$, and vectors $\delta p^{(0)},\delta
q^{(0)}\in\mathbbm{C}^{n}$ such that
$\displaystyle\widehat{A}^{(1)\dagger}=U^{(1)}(A+\delta
A^{(0)\dagger})U^{(1)*},$ (263)
$\displaystyle\widehat{p}^{\,(1)}=U^{(1)}(p+\delta p^{(0)}),$ (264) and
$\displaystyle\widehat{q}^{\,(1)}=U^{(1)}(q+\delta q^{(0)}),$ (265) where
${\lVert\delta A^{(0)\dagger}\rVert}\lesssim{\lVert A\rVert}\mathit{u}$,
${\lVert\delta p^{(0)}\rVert}\lesssim{\lVert p\rVert}\mathit{u}$, and
${\lVert\delta q^{(0)}\rVert}\lesssim{\lVert q\rVert}\mathit{u}$. Let
$\widehat{A}^{(1)}$ be equal to $\widehat{A}^{(1)\dagger}$, except that the
entry in the $(1,2)$ position of $\widehat{A}^{(1)}$ is equal to
$-\widehat{p}_{1}^{\,(1)}\widehat{q}_{2}^{\,(1)*}$, so that
$\widehat{A}_{1,2}^{(1)}+\widehat{p}_{1}^{\,(1)}\widehat{q}_{2}^{\,(1)*}=0$.
Since
${\lvert\widehat{A}_{1,2}^{(1)\dagger}+\widehat{p}_{1}^{\,(1)}\widehat{q}_{2}^{\,(1)*}\rvert}<\epsilon$,
we have
$\displaystyle{\lVert\widehat{A}^{(1)\dagger}-\widehat{A}^{(1)}\rVert}_{F}<\epsilon,$
(266) where ${\lVert\cdot\rVert}_{F}$ denotes the Frobenius norm, and since
$\epsilon\lesssim{\lVert A\rVert}\mathit{u}$,
$\displaystyle{\lVert\widehat{A}^{(1)\dagger}-\widehat{A}^{(1)}\rVert}\lesssim{\lVert
A\rVert}\mathit{u}.$ (267) Letting  $\displaystyle\delta A^{(0)}=\delta
A^{(0)\dagger}+U^{(1)*}(\widehat{A}^{(1)}-\widehat{A}^{(1)\dagger})U^{(1)}$
(268) and combining (263) and (268), we have
$\displaystyle\widehat{A}^{(1)}=U^{(1)}(A+\delta A^{(0)})U^{(1)*},$ (269)
where ${\lVert\delta A^{(0)}\rVert}\lesssim{\lVert A\rVert}\mathit{u}$ by
(267). Clearly, since
$\widehat{A}^{\,(1)}_{1,2}+\widehat{p}_{1}^{\,(1)}\widehat{q}_{2}^{\,(1)*}=0$
and $\widehat{A}^{(1)}+\widehat{p}^{\,(1)}\widehat{q}^{\,(1)*}$ is lower
Hessenberg,
$\widehat{\lambda}_{1}=\widehat{A}_{1,1}^{\,(1)}+\widehat{p}_{1}^{\,(1)}\widehat{q}_{1}^{\,(1)*}$
is an eigenvalue of
$\widehat{A}^{(1)}+\widehat{p}^{\,(1)}\widehat{q}^{\,(1)*}$. Thus, from (264),
(265), and (269), we see that $\widehat{\lambda}_{1}$ is an eigenvalue of
$(A+\delta A^{(0)})+(p+\delta p^{(0)})(q+\delta q^{(0)})^{*}$. Now suppose
that we deflate the matrix, and perform QR iterations on the submatrix
$\widehat{A}^{\,(1)}_{2:n,2:n}+\widehat{p}_{2:n}^{\,(1)}\widehat{q}_{2:n}^{\,(1)*}$,
until the entry in the $(1,2)$ position of the deflated matrix is less than
$\epsilon$. Let $\widehat{\undertilde{d}}^{(2)}\in\mathbbm{C}^{n-1}$,
$\widehat{\undertilde{\beta}}^{(2)\dagger}\in\mathbbm{C}^{n-2}$,
$\widehat{\undertilde{p}}^{(2)}\in\mathbbm{C}^{n-1}$, and
$\widehat{\undertilde{q}}^{(2)}\in\mathbbm{C}^{n-1}$ denote the resulting
outputs, and let
$\widehat{\undertilde{A}}^{(2)\dagger}\in\mathbbm{C}^{(n-1)\times(n-1)}$ be
the resulting matrix, defined by formula (230) (again making the obvious
substitutions). By Lemma 4.5, there exist a unitary matrix
$\undertilde{U}^{(2)}\in\mathbbm{C}^{(n-1)\times(n-1)}$, a matrix
$\delta\undertilde{\widehat{A}}^{(1)\dagger}\in\mathbbm{C}^{(n-1)\times(n-1)}$,
and vectors
$\delta\widehat{\undertilde{p}}^{(1)},\delta\widehat{\undertilde{q}}^{(1)}\in\mathbbm{C}^{n-1}$
such that
$\displaystyle\widehat{\undertilde{A}}^{(2)\dagger}=\undertilde{U}^{(2)}(\widehat{A}^{(1)}_{2:n,2:n}+\delta\widehat{\undertilde{A}}^{(1)\dagger})\undertilde{U}^{(2)*},$
(270)
$\displaystyle\widehat{\undertilde{p}}^{\,(2)}=U^{(2)}(\widehat{p}_{2:n}^{\,(1)}+\delta\widehat{\undertilde{p}}^{(1)}),$
(271) and
$\displaystyle\widehat{\undertilde{q}}^{\,(2)}=U^{(2)}(\widehat{q}_{2:n}^{\,(1)}+\delta\widehat{\undertilde{q}}^{(1)}),$
(272) where
${\lVert\delta\widehat{\undertilde{A}}^{(1)\dagger}\rVert}\lesssim{\lVert
A\rVert}\mathit{u}$,
${\lVert\delta\widehat{\undertilde{p}}^{(1)}\rVert}\lesssim{\lVert
p\rVert}\mathit{u}$, and
${\lVert\delta\widehat{\undertilde{q}}^{(1)}\rVert}\lesssim{\lVert
q\rVert}\mathit{u}$. Like before, let
$\widehat{\undertilde{A}}^{(2)}\in\mathbbm{C}^{(n-1)\times(n-1)}$ be equal to
$\widehat{\undertilde{A}}^{(2)\dagger}$, except that the entry in the $(1,2)$
position of $\widehat{\undertilde{A}}^{(2)}$ is equal to
$-\widehat{\undertilde{p}}_{1}^{(2)}\widehat{\undertilde{q}}_{2}^{(2)*}$, so
that
$\widehat{\undertilde{A}}_{1,2}^{(2)}+\widehat{\undertilde{p}}_{1}^{(2)}\widehat{\undertilde{q}}_{2}^{(2)*}=0$.
Since
${\lvert\widehat{\undertilde{A}}_{1,2}^{(2)\dagger}+\widehat{\undertilde{p}}_{1}^{(2)}\widehat{\undertilde{q}}_{2}^{(2)*}\rvert}<\epsilon$,
we have
$\displaystyle{\lVert\widehat{\undertilde{A}}^{(2)\dagger}-\widehat{\undertilde{A}}^{(2)}\rVert}_{F}<\epsilon,$
(273) where ${\lVert\cdot\rVert}_{F}$ denotes the Frobenius norm, and since
$\epsilon\lesssim{\lVert A\rVert}\mathit{u}$,
$\displaystyle{\lVert\widehat{\undertilde{A}}^{(2)\dagger}-\widehat{\undertilde{A}}^{(2)}\rVert}\lesssim{\lVert
A\rVert}\mathit{u}.$ (274) Letting
$\displaystyle\delta\widehat{\undertilde{A}}^{(1)}=\delta\widehat{\undertilde{A}}^{(1)\dagger}+\undertilde{U}^{(2)*}(\widehat{\undertilde{A}}^{(2)}-\widehat{\undertilde{A}}^{(2)\dagger})\undertilde{U}^{(2)},$
(275) we have
$\displaystyle\widehat{\undertilde{A}}^{(2)}=\undertilde{U}^{(2)}(\widehat{A}_{2:n,2:n}^{(1)}+\delta\widehat{\undertilde{A}}^{(1)})\undertilde{U}^{(2)*},$
(276) where ${\lVert\delta\widehat{\undertilde{A}}^{(1)}\rVert}\lesssim{\lVert
A\rVert}\mathit{u}$. Since
$\widehat{\undertilde{A}}^{(2)}_{1,2}+\widehat{\undertilde{p}}_{1}^{(1)}\widehat{\undertilde{q}}_{2}^{(1)*}=0$
and
$\widehat{\undertilde{A}}^{(2)}+\widehat{\undertilde{p}}^{(2)}\widehat{\undertilde{q}}^{(2)*}$
is lower Hessenberg,
$\widehat{\lambda}_{2}=\widehat{\undertilde{A}}_{1,1}^{(2)}+\widehat{\undertilde{p}}_{1}^{(2)}\widehat{\undertilde{q}}_{1}^{(2)*}$
is an eigenvalue of
$\widehat{\undertilde{A}}^{(2)}+\widehat{\undertilde{p}}^{(2)}\widehat{\undertilde{q}}^{(2)*}$.
Now define the unitary matrix $U^{(2)}\in\mathbbm{C}^{n\times n}$ by the
formula  $\displaystyle U^{(2)}=\left[\begin{array}[]{c|ccc}1&0&\cdots&0\\\
\hline\cr
0&\lx@intercol\hfil\hbox{\multirowsetup$\undertilde{U}^{(2)}$}\hfil\lx@intercol\\\
\vdots&&&\\\ 0&&&\end{array}\right],$ (281) the matrix
$\delta\widehat{A}^{(1)}\in\mathbbm{C}^{n\times n}$ by
$\displaystyle\delta\widehat{A}^{(1)}=\left[\begin{array}[]{c|ccc}0&0&\cdots&0\\\
\hline\cr
0&\lx@intercol\hfil\hbox{\multirowsetup$\delta\widehat{\undertilde{A}}^{(1)}$}\hfil\lx@intercol\\\
\vdots&&&\\\ 0&&&\end{array}\right],$ (286) and the vectors
$\delta\widehat{p}^{(1)},\delta\widehat{q}^{(1)}\in\mathbbm{C}^{n}$, by
$\displaystyle\delta\widehat{p}^{\,(1)}=\left[\begin{array}[]{c}0\\\
\hline\cr\delta\widehat{\undertilde{p}}^{(1)}\end{array}\right],$ (289) and
$\displaystyle\delta\widehat{q}^{\,(1)}=\left[\begin{array}[]{c}0\\\
\hline\cr\delta\widehat{\undertilde{q}}^{(1)}\end{array}\right].$ (292)
Clearly, ${\lVert\delta A^{(1)}\rVert}\lesssim{\lVert A\rVert}\mathit{u}$,
${\lVert\delta\widehat{p}^{\,(1)}\rVert}\lesssim{\lVert p\rVert}\mathit{u}$,
and ${\lVert\delta\widehat{q}^{\,(1)}\rVert}\lesssim{\lVert
q\rVert}\mathit{u}$. Let $\widehat{A}^{(2)}\in\mathbbm{C}^{n\times n}$ be
defined by  $\displaystyle\widehat{A}^{(2)}=U^{(2)}(\widehat{A}^{(1)}+\delta
A^{(1)})U^{(2)*},$ (293) and
$\widehat{p}^{\,(2)},\widehat{q}^{\,(2)}\in\mathbbm{C}^{n}$ by
$\displaystyle\widehat{p}^{\,(2)}=U^{(2)}(\widehat{p}^{\,(1)}+\delta\widehat{p}^{\,(1)}),$
(294) and
$\displaystyle\widehat{q}^{\,(2)}=U^{(2)}(\widehat{q}^{\,(1)}+\delta\widehat{q}^{\,(1)}).$
(295) We first notice that
$\widehat{A}_{1,\ell}^{\,(2)}+\widehat{p}_{1}^{\,(2)}\widehat{q}_{\ell}^{\,(2)*}=0$
for $\ell=2,3,\ldots,n$. Next, we observe that
$\widehat{A}_{1,1}^{(2)}+\widehat{p}_{1}^{\,(2)}\widehat{q}_{1}^{\,(2)*}=\widehat{A}_{1,1}^{(1)}+\widehat{p}_{1}^{\,(1)}\widehat{q}_{1}^{\,(1)*}=\widehat{\lambda}_{1}$;
therefore, $\widehat{\lambda}_{1}$ is an eigenvalue of
$\widehat{A}^{(2)}+\widehat{p}^{\,(2)}\widehat{q}^{\,(2)*}$. We then observe
that
$(\widehat{A}^{(2)}+\widehat{p}^{\,(2)}\widehat{q}^{\,(2)*})_{2:n,2:n}=\widehat{\undertilde{A}}^{(2)}+\widehat{\undertilde{p}}^{(2)}\widehat{\undertilde{q}}^{(2)*}$;
therefore, $\widehat{\lambda}_{2}$ is an eigenvalue of
$\widehat{A}^{(2)}+\widehat{p}^{\,(2)}\widehat{q}^{\,(2)*}$. Finally, letting
$U=U^{(2)}U^{(1)}$ and substituting (264), (265), and (269) into (293)–(295)
and expanding, it is straightforward to show that there exist matrices $\delta
A\in\mathbbm{C}^{n\times n}$ and vectors $\delta p,\delta q\in\mathbbm{C}^{n}$
such that
$\displaystyle\widehat{A}^{(2)}+\widehat{p}^{\,(2)}\widehat{q}^{\,(2)*}=U(A+\delta
A)U^{*}+U(p+\delta p)(q+\delta q)^{*}U^{*},$ (296) where ${\lVert\delta
A\rVert}\lesssim{\lVert A\rVert}\mathit{u}$, ${\lVert\delta
p\rVert}\lesssim{\lVert p\rVert}\mathit{u}$, and ${\lVert\delta
q\rVert}\lesssim{\lVert q\rVert}\mathit{u}$. Therefore,
$\widehat{\lambda}_{1}$ and $\widehat{\lambda}_{2}$ are eigenvalues of the
matrix  $\displaystyle(A+\delta A)+(p+\delta p)(q+\delta q),$ (297) where
${\lVert\delta A\rVert}\lesssim{\lVert A\rVert}\mathit{u}$, ${\lVert\delta
p\rVert}\lesssim{\lVert p\rVert}\mathit{u}$, and ${\lVert\delta
q\rVert}\lesssim{\lVert q\rVert}\mathit{u}$. The same proof can be repeated
inductively to show this for all
$\widehat{\lambda}_{1},\widehat{\lambda}_{2},\ldots,\widehat{\lambda}_{n}$.
$\blacksquare$

The following theorem states that our explicit shifted QR algorithm is
componentwise backward stable, for those eigenvalues for which the shifts are
small.

######  Theorem 4.7 (Explicit shifted QR).

Suppose that $A\in\mathbbm{C}^{n\times n}$ is a Hermitian matrix, that
$p,q\in\mathbbm{C}^{n}$, and that $A+pq^{*}$ is lower Hessenberg. Let $d$ and
$\beta$ denote the diagonal and superdiagonal of $A$, respectively. Suppose
that Algorithm 4 is carried out in floating point arithmetic with
$\epsilon\lesssim{\lVert A\rVert}\mathit{u}$, and suppose that $\mu^{(\ell)}$
is the largest total shift encountered at any point during the course of the
algorithm from $i=1,2,\ldots,\ell$ in the outer loop. Let
$\widehat{\lambda}_{1},\widehat{\lambda}_{2},\ldots,\widehat{\lambda}_{n}$
denote the outputs of the algorithm. Then, for each $\ell=1,2,\ldots,n$, there
exist a matrix $\delta A\in\mathbbm{C}^{n\times n}$ and vectors $\delta
p,\delta q\in\mathbbm{C}^{n}$ such that
$\widehat{\lambda}_{1},\widehat{\lambda}_{2},\ldots,\widehat{\lambda}_{\ell}$
are exact eigenvalues of the matrix

$\displaystyle(A+\delta A)+(p+\delta p)(q+\delta q)^{*},$ (298)

where ${\lVert\delta A\rVert}\lesssim({\lVert
A\rVert}+{\lvert\mu^{(\ell)}\rvert})\mathit{u}$, ${\lVert\delta
p\rVert}\lesssim{\lVert p\rVert}\mathit{u}$, and ${\lVert\delta
q\rVert}\lesssim{\lVert q\rVert}\mathit{u}$.

Proof. The proof is essentially identical to the proof of Theorem 4.6, and we
omit it. $\blacksquare$

######  Remark 4.1.

Notice that Theorems 4.6 and 4.7 do not make any mention of convergence. What
they say is that, if the algorithm converges, then it is componentwise
backward stable. We observe that, in practice, Algorithm 4 always converges
rapidly, at least quadratically, for $\epsilon\approx{\lVert
A\rVert}\mathit{u}$.

######  Remark 4.2.

Notice that the bound on $\delta A$ in Theorem 4.7 involves $\mu^{(\ell)}$,
which is the largest total shift encountered at any point during the
calculation of
$\widehat{\lambda}_{1},\widehat{\lambda}_{2},\ldots,\widehat{\lambda}_{\ell}$.
While this bound appears weaker than the corresponding bound in Theorem 4.6,
in practice it turns out to be essentially the same, as follows. We can always
assume that $A$ is much smaller than $p$, or $q$, or both; if this isn’t the
case, then componentwise stability no longer has any special meaning, since it
follows immediately from the usual Bauer-Fike perturbation bounds (see [7]).
Furthermore, we tend to be interested in the componentwise backward stability
of small eigenvalues $\widehat{\lambda}_{i}$, where
${\lvert\widehat{\lambda}_{i}\rvert}\approx{\lVert A\rVert}$. If we perform a
few iterations of unshifted QR on the matrix, then the eigenvalues of the top-
left $2\times 2$ block will approach the two smallest eigenvalues of the
matrix (recalling that our QR algorithm works with lower Hessenberg matrices).
If we now use Algorithm 4, we’ll find that the initial shift is small and, as
a result, all the eigenvalues are computed roughly in order from smallest to
largest. This means that
${\lvert\mu^{(i)}\rvert}\approx{\lvert\widehat{\lambda}_{i}\rvert}$ and
(approximately)
$\widehat{\lambda}_{1}<\widehat{\lambda}_{2}<\cdots<\widehat{\lambda}_{i}$.
For $\widehat{\lambda}_{i}$ such that
${\lvert\widehat{\lambda}_{i}\rvert}\approx{\lVert A\rVert}$, we have then
that the bound ${\lVert\delta A\rVert}\lesssim({\lVert
A\rVert}+{\lvert\mu^{(i)}\rvert})\mathit{u}$ becomes ${\lVert\delta
A\rVert}\lesssim{\lVert A\rVert}\mathit{u}$. Finally, we point out that the
dependence of the bound on $\mu^{(\ell)}$ could likely be removed entirely by
reformulating our QR algorithm as an implicit method.

##  5 Numerical Results

In this section, we demonstrate the componentwise backward stability of our
shifted QR algorithm (see Algorithm 4) by illustrating its stability when it
is used as a rootfinding algorithm (see Sections 2.3 and 2.4). Consider a
polynomial $p(x)$ of order $n$, not necessarily monic, expressed in a
Chebyshev polynomial basis

$\displaystyle p(x)=\sum_{j=0}^{n}a_{j}T_{j}(x),$ (299)

where $a_{j}\in\mathbbm{R}$ and $T_{j}(x)$ is the Chebyshev polynomial of
order $j$. By Theorem 2.7 and Remark 2.1, we have that if the eigenvalues of
the linearization (31), where $c_{j}=a_{j}/a_{n}$, $j=0,1,\ldots,n$, are
computed by a componentwise backward stable algorithm, then the computed roots
$\widehat{x}_{1},\widehat{x}_{2},\ldots,\widehat{x}_{n}$ are the exact roots
of the perturbed polynomial

$\displaystyle p(x)+\delta p(x)=\sum_{j=0}^{n}(a_{j}+\delta a_{j})T_{j}(x),$
(300)

where

$\displaystyle\frac{{\lVert\delta a\rVert}}{{\lVert
a\rVert}}\lesssim\mathit{u}.$ (301)

By applying our QR algorithm to linearizations of various polynomials $p(x)$,
we demonstrate our algorithm’s componentwise backward stability by showing
that the bound (301) always holds.

We estimate the size of the backward error $\delta a$ in the coefficients by
using the following observation (see the discussion accompanying Table 1 in
[26]). By the definition of $p(x)+\delta p(x)$, we have that $(p+\delta
p)(\widehat{x}_{i})=0$ for $i=1,2,\ldots,n$. From (299) and (300), it follows
that

$\displaystyle p(\widehat{x}_{i})=p(\widehat{x}_{i})-(p+\delta
p)(\widehat{x}_{i})=-\sum_{j=0}^{n}\delta a_{j}T_{j}(\widehat{x}_{j}).$ (302)

Since $-1\leq T_{j}(x)\leq 1$ for all $j$ when $x\in[-1,1]$, we have

$\displaystyle p(\widehat{x}_{i})\approx{\lVert\delta a\rVert},$ (303)

whenever $\widehat{x}_{i}\in\mathbbm{C}$ is not too far from the interval
$[-1,1]$.

Even though $\widehat{x}_{i}$ is already a floating point number, the value
$p(\widehat{x}_{i})$ cannot be computed exactly in floating point arithmetic.
Letting $\widehat{p}(\widehat{x}_{i})$ denote the approximation to
$p(\widehat{x}_{i})$ computed in floating point arithmetic, we know that

$\displaystyle\widehat{p}(\widehat{x}_{i})\approx
p(\widehat{x}_{i})+\kappa(p;\widehat{x}_{i})\mathit{u},$ (304)

where

$\displaystyle\kappa(p;\widehat{x}_{i})={\lvert\widehat{x}_{i}\rvert}{\lvert
p^{\prime}(\widehat{x}_{i})\rvert}$ (305)

is the absolute condition number of $p(x)$ at $x=\widehat{x}_{i}$. When
$\kappa(p;\widehat{x}_{i})$ is large, the error in evaluating
$p(\widehat{x}_{i})$ dominates, while when $\kappa(p;\widehat{x}_{i})$ is of
modest size, we have $\widehat{p}(\widehat{x}_{i})\approx p(\widehat{x}_{i})$.
In this section, we investigate the quantity

$\displaystyle\eta(p;\widehat{x}_{i})=\frac{\widehat{p}(\widehat{x}_{i})}{\max(\kappa(p;\widehat{x}_{i}),{\lVert
a\rVert})},$ (306)

for various polynomials $p(x)$. When $\kappa(p;\widehat{x}_{i})\geq{\lVert
a\rVert}$, we have

$\displaystyle{\lvert\eta(p;\widehat{x}_{i})\rvert}=\frac{{\lvert\widehat{p}(\widehat{x}_{i})\rvert}}{\kappa(p;\widehat{x}_{i})}\approx\frac{{\lvert
p(\widehat{x}_{i})\rvert}}{\kappa(p;\widehat{x}_{i})}+u\leq\frac{{\lvert
p(\widehat{x}_{i})\rvert}}{{\lVert a\rVert}}+u.$ (307)

When $\kappa(p;\widehat{x}_{i})\leq{\lVert a\rVert}$,

$\displaystyle{\lvert\eta(p;\widehat{x}_{i})\rvert}=\frac{{\lvert\widehat{p}(\widehat{x}_{i})\rvert}}{{\lVert
a\rVert}}\approx\frac{{\lvert p(\widehat{x}_{i})\rvert}}{{\lVert
a\rVert}}+\frac{\kappa(p;\widehat{x}_{i})}{{\lVert a\rVert}}u\leq\frac{{\lvert
p(\widehat{x}_{i})\rvert}}{{\lVert a\rVert}}+u.$ (308)

Thus, if our QR algorithm is indeed componentwise backward stable and (301) is
satisfied, then, by (303), we expect to find that

$\displaystyle\eta(p;\widehat{x}_{i})\approx\mathit{u},$ (309)

for all polynomials $p(x)$.

For $p(\widehat{x}_{i})$ to be a good approximation to ${\lVert\delta
a\rVert}$ (see (303)), we stated that $\widehat{x}_{i}\in\mathbbm{C}$ should
be “not too far from the interval $[-1,1]$.” We make this notion precise as
follows. Let $z_{i}$, $i=1,2,\ldots,n$ denote the exact roots of the order-$n$
polynomial $p(x)$, and let $\widehat{z}_{i}$ denote the computed roots. We
select roots close to the interval $[-1,1]$ by choosing some $\delta>0$ (for
example, $\delta=10^{-3}$), and letting $\widehat{x}_{i}\in\mathbbm{R}$ denote
the real part of all roots $\widehat{z}_{i}$ that are inside the rectangle

$\displaystyle\\{z\in\mathbbm{C}:1-\delta<\operatorname{Re}(z)<1+\delta,-\delta<\operatorname{Im}(z)<\delta\\}.$
(310)

If the polynomial $p(x)$ has the real root $z_{i}$, then taking the real part
of $\widehat{z}_{i}$ will not result in any additional error. The number of
real roots inside the region (310) will often be less than the order $n$, and
we denote the number of such roots by $n_{\text{roots}}$.

In our numerical experiments, we compute the eigenvalues of the colleague
matrix using three different algorithms: our Algorithm 4; MATLAB’s eig
function; and MATLAB’s eig function with balancing turned off (using the
option ’nobalance’), which we call eig_nb. For our experiments in extended
(quadruple) precision, we use the Advanpix Multiprecision Computing Toolbox
and its implementation of eig (see [1]). Since the Advanpix Multiprecision
Computing Toolbox’s eig function always balances the matrix, and does not
support the ’nobalance’ option, we omit the test of eig_nb in extended
precision.

For each example, we report the degree of the underlying polynomial, the order
$n$ of the Chebyshev expansion used to approximate it, the size of the vector
$c$ in the Euclidean norm, the Frobenius norm of the completely balanced
colleague matrix, which we denote by $\texttt{bal}(C)$, the number
$n_{\text{roots}}$ of computed roots inside the region (310) for the given
value of $\delta>0$, the size $\max_{i}{\lvert z_{i}\rvert}$ of the largest
complex root of the colleague matrix, and the value of
$\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$, where the maximum is taken
over all of the real parts of the computed roots inside the region (310).

We implemented our algorithm in FORTRAN 77, and compiled it using
Lahey/Fujitsu Fortran 95 Express, Release L6.20e. For the timing experiments,
the Fortran codes were compiled using the Intel Fortran Compiler, version
19.0.2.187, with the -fast flag. The MATLAB experiments were performed using
MATLAB R2019b, version 9.7.0.1190202, and the extended precision MATLAB
experiments were performed in quadruple precision (mp.Digits(34)) using the
Advanpix Multiprecision Computing Toolbox, version 4.8.0, Build 14100. All
experiments we conducted on a ThinkPad laptop, with 16GB of RAM and an Intel
Core i7-8550U CPU.

###  5.1 $p_{\text{rand}}(x)$: Polynomials with Random Coefficients

Following [14], we construct polynomials $p_{\text{rand}}(x)$ by sampling
Chebyshev expansion coefficients $a_{i}$ independently from a standard normal
distribution, so that $a_{i}\sim N(0,1)$, for $i=0,1,\ldots,n-1$. Then, we
choose the desired value of ${\lVert c\rVert}$ by setting $a_{n}={\lVert
a\rVert}/{\lVert c\rVert}$, so that the vector of coefficients $c$ appearing
in the colleague matrix (31), where $c_{i}=a_{i}/a_{n}$ for $i=1,2,\ldots,n$,
has the specified norm. For this example, we choose $n=30$ and set
$\delta=10^{-5}$ to extract the real roots (see formula (310)).

We report the results in Figure 3. We see that our algorithm shows the
expected backward stability over the entire range of ${\lVert c\rVert}$, while
MATLAB’s eig, both balanced and unbalanced, shows the expected growth with
${\lVert c\rVert}$ (see the discussion in Section 2.4). Interestingly, for
this example, balancing appears to only improve the error by an order of
magnitude or two, while leaving the growth in the error with respect to
${\lVert c\rVert}$ unchanged. This turns out be completely consistent the
following explanation. Balancing the colleague matrix can reduce the magnitude
of the all elements from ${\lVert c\rVert}$ to ${\lVert
c\rVert}^{\frac{1}{2}}$, except for the element in the $(n,n)$-position, which
balancing cannot change. In this example, all of the elements of the vector
$c$ are around the same size as ${\lVert c\rVert}$, so there are $n$ large
elements in the colleague matrix. Thus, balancing reduces the number of large
elements of size ${\lVert c\rVert}$ from $n$ to $1$, resulting in an $n$-fold
reduction in the norm of the matrix. In this example, $n=30$, which
corresponds well with the approximately $30$-fold reduction in error due to
balancing that we observe in Figure 3.

Figure 3:  The values of $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$ for
various values of ${\lVert c\rVert}$, in double precision (left) and quadruple
precision (right), for the polynomials $p_{\text{rand}}(x)$ of order $n=30$,
computed by our algorithm, eig, and eig_nb, with $\delta=10^{-5}$. The values
are indicated for our algorithm with purple crosses ($+$), for eig with green
x’s ($\times$), and for eig_nb with blue stars ($\star$).

We found that the colleague matrix has a single large eigenvalue of the size
${\lVert c\rVert}$, and the rest of the eigenvalues are small. This is not
surprising, since there is an entry of size ${\lVert c\rVert}$ is the
$(n,n)$-position of the matrix (from which it follows that
$Ce_{n}\approx{\lVert c\rVert}e_{n}$). Thus, for all three algorithms,
$\max_{i}{\lvert\widehat{z}_{i}\rvert}\approx{\lVert c\rVert}$.

###  5.2 $p_{\text{wilk}}(x)$: Wilkinson’s Polynomial

Here we consider the famous Wilkinson polynomial, normalized so that all of
its roots are inside the interval $[-1,1]$, defined by the formula

$\displaystyle
p_{\text{wilk}}(x)=\prod_{i=1}^{m}\Bigl{(}x-\Bigl{(}\frac{2i}{m+1}-1\Bigr{)}\Bigr{)}.$
(311)

We construct an order-$n$ Chebyshev expansion of this degree-$m$ polynomial,
sampling it at $n$ Chebshev points and applying a linear transformation to
obtain the expansion coefficients (see [36]). We then compute the eigenvalues
of the colleague matrix, and set $\delta=10^{-3}$ to extract the real roots
(see formula (310)). The results of our numerical experiment are shown in
Tables 1 and 2. We observe that our algorithm is backward stable, while eig
loses accuracy whenever the roots of the colleague matrix are large. Plots of
the real and complex roots of the order-$100$ Chebyshev expansion are shown
for various degrees of $p_{\text{wilk}}(x)$ in Figure 4. When the order of the
Wilkinson polynomial becomes large, spurious real roots begin appearing in the
middle of the interval $[-1,1]$. It turns out that the roots computed by our
algorithm are still backward stable, even in this situation; the function is
so small near the middle of the interval that a small relative perturbation in
the Chebyshev coefficients causes additional roots to appear.

Figure 4:  The roots of the Chebyshev expansion of order $100$ of the Wilkinson polynomial $p_{\text{wilk}}(x)$, of various degrees, computed by our algorithm. The complex roots $\widehat{z}_{i}$ are plotted with purple crosses ($+$) and the real roots $\widehat{x}_{i}$ are plotted with blue stars ($\star$). The Wilkinson polynomial has, in order from top to bottom, orders $24$, $34$, $44$, and $54$. Observe that the spurious complex roots are well-separated from the interval $[-1,1]$ when the order is low, but eventually meet the interval when the order is large.  |  |  |  |  | eig |  Algorithm 4
---|---|---|---|---|---|---
Degree | $n$ | ${\lVert c\rVert}$ | ${\lVert\texttt{bal}(C)\rVert}$ | $\max_{i}{\lvert z_{i}\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$
14 | 100 | $0.45\cdot 10^{14}$ | $0.32\cdot 10^{2}$ | $0.11\cdot 10^{1}$ | 14 | $0.16\cdot 10^{-13}$ | 14 | $0.71\cdot 10^{-14}$
24 | 24 | $0.95\cdot 10^{4}$ | $0.65\cdot 10^{1}$ | $0.92\cdot 10^{0}$ | 24 | $0.15\cdot 10^{-14}$ | 24 | $0.32\cdot 10^{-14}$
| 25 | $0.66\cdot 10^{15}$ | $0.35\cdot 10^{11}$ | $0.35\cdot 10^{11}$ |  14† | $0.11\cdot 10^{-4}$ | 24 | $0.19\cdot 10^{-14}$
| 26 | $0.22\cdot 10^{15}$ | $0.11\cdot 10^{6}$ | $0.76\cdot 10^{5}$ | 24 | $0.51\cdot 10^{-10}$ | 24 | $0.24\cdot 10^{-14}$
| 27 | $0.40\cdot 10^{16}$ | $0.66\cdot 10^{4}$ | $0.38\cdot 10^{4}$ | 24 | $0.90\cdot 10^{-12}$ | 24 | $0.19\cdot 10^{-14}$
| 28 | $0.12\cdot 10^{15}$ | $0.37\cdot 10^{3}$ | $0.17\cdot 10^{3}$ | 24 | $0.58\cdot 10^{-13}$ | 24 | $0.14\cdot 10^{-14}$
| 100 | $0.38\cdot 10^{14}$ | $0.24\cdot 10^{2}$ | $0.11\cdot 10^{1}$ | 24 | $0.10\cdot 10^{-13}$ | 24 | $0.24\cdot 10^{-14}$
34 | 100 | $0.34\cdot 10^{14}$ | $0.17\cdot 10^{2}$ | $0.11\cdot 10^{1}$ | 34 | $0.14\cdot 10^{-13}$ | 34 | $0.12\cdot 10^{-13}$
44 | 100 | $0.32\cdot 10^{14}$ | $0.16\cdot 10^{2}$ | $0.11\cdot 10^{1}$ | 44 | $0.53\cdot 10^{-13}$ | 44 | $0.41\cdot 10^{-14}$
54 | 100 | $0.30\cdot 10^{14}$ | $0.16\cdot 10^{2}$ | $0.11\cdot 10^{1}$ | 60 | $0.15\cdot 10^{-13}$ | 60 | $0.28\cdot 10^{-13}$
Table 1:  The results of computing the roots of the Wilkinson polynomial $p_{\text{wilk}}(x)$, using our algorithm and eig, with $\delta=10^{-3}$. †The error was so large here that some roots were outside of the region (310).  |  |  |  |  | eig |  Algorithm 4
---|---|---|---|---|---|---
Degree | $n$ | ${\lVert c\rVert}$ | ${\lVert\texttt{bal}(C)\rVert}$ | $\max_{i}{\lvert z_{i}\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$ | $n_{\text{roots}}$ | $\max_{i}{\lvert\eta(p;\widehat{x}_{i})\rvert}$
14 | 100 | $0.12\cdot 10^{34}$ | $0.38\cdot 10^{3}$ | $0.14\cdot 10^{1}$ | 14 | $0.22\cdot 10^{-31}$ | 14 | $0.15\cdot 10^{-32}$
24 | 24 | $0.95\cdot 10^{4}$ | $0.65\cdot 10^{1}$ | $0.92\cdot 10^{0}$ | 24 | $0.44\cdot 10^{-32}$ | 24 | $0.66\cdot 10^{-33}$
| 25 | $0.95\cdot 10^{33}$ | $0.50\cdot 10^{29}$ | $0.50\cdot 10^{29}$ |  12† | $0.29\cdot 10^{-4}$ | 24 | $0.78\cdot 10^{-32}$
| 26 | $0.16\cdot 10^{33}$ | $0.94\cdot 10^{14}$ | $0.66\cdot 10^{14}$ | 24 | $0.11\cdot 10^{-19}$ | 24 | $0.18\cdot 10^{-32}$
| 27 | $0.11\cdot 10^{36}$ | $0.22\cdot 10^{11}$ | $0.11\cdot 10^{11}$ | 24 | $0.55\cdot 10^{-23}$ | 24 | $0.76\cdot 10^{-32}$