Proofs: QZ algorithm

The QZ algorithm keeps an exact unitary equivalence, preserves the Hessenberg-triangular structure, and terminates with a generalized Schur form full

This proof shows: correctness, termination.

Exact arithmetic, except that the entries set to zero by the tests of lines L4 and L6 are tracked as perturbations. Notation of the page: $H$, $T$, $Q$, $Z$ are the working matrices, indices run from $1$ to $n$, $\varepsilon = 2u$, and a rotation is $G = \begin{pmatrix} c & s \\ -\bar s & c \end{pmatrix}$ with $c \ge 0$ real and $c^2 + |s|^2 = 1$, acting on two rows (from the left) or two columns (from the right). The quadratic convergence of the shifted iteration is not proved here; see the remark at the end.

Lemma 1 (rotations). For any $f, g$, givens(f, g) returns a unitary $G$ with $G (f, g)^T = (r, 0)^T$. Consequently: (a) a row rotation of rows $i, j$ built from $(f, g) = (m_{ik}, m_{jk})$ zeroes $m_{jk}$; (b) a column rotation of columns $i, j$ with the matrix $G$ built from givens(y, x), where $(x, y) = (m_{ki}, m_{kj})$, maps that row pair to $(0, r)$.

Proof. If $g = 0$, $G = I$. If $f = 0$, $c = 0$, $s = 1$ and $G (0, g)^T = (g, 0)^T$. Otherwise, with $\nu = \sqrt{|f|^2 + |g|^2}$ and $\phi = f / |f|$, $c = |f| / \nu$ and $s = \phi \bar g / \nu$; then $c^2 + |s|^2 = 1$, so $G G^* = I$, and the second component is $-\bar s f + c g = -\bar\phi g f / \nu + |f| g / \nu = -|f| g / \nu + |f| g / \nu = 0$, using $\bar\phi f = |f|$. For (b), the new entry in column $i$ is $x c - \bar s y$, which is the second component of $G (y, x)^T$, hence zero. $\square$

Lemma 2 (the invariant). Let $H_0 = A$, $T_0 = B$, $Q = Z = I$ at the start. After any sequence of the operations of the algorithm, $Q$ and $Z$ are unitary and $Q^* (A + E) Z = H$, $Q^* (B + F) Z = T$, where $E = \sum_\tau Q_\tau (\delta_\tau e_{k_\tau} e_{k_\tau - 1}^T) Z_\tau^*$ and $F = \sum_\sigma Q_\sigma (\gamma_\sigma e_{k_\sigma} e_{k_\sigma}^T) Z_\sigma^*$ collect the entries $\delta_\tau$ of $H$ and $\gamma_\sigma$ of $T$ that lines L4 and L6 set to zero, with $Q_\tau, Z_\tau$ the values of $Q, Z$ at that moment. Hence $\|E\|_F \le \sum_\tau |\delta_\tau|$ and $\|F\|_F \le \sum_\sigma |\gamma_\sigma|$.

Proof. Induction on the operations. A row rotation $G$ replaces $H, T, Q$ by $G H, G T, Q G^*$; then $(Q G^*)^* (A + E) Z = G H$, and $Q G^*$ is unitary. A column rotation $R$ replaces $H, T, Z$ by $H R, T R, Z R$, and $Q^* (A + E) (Z R) = H R$. The assignments to zero of entries that the rotation has just made zero (lines L1, L7, L9, L10) change nothing in exact arithmetic. Setting $h_{k,k-1} = \delta$ to zero replaces $H$ by $H - \delta e_k e_{k-1}^T = Q^* (A + E + Q \delta e_k e_{k-1}^T Z^*) Z$, which extends $E$ by one term; likewise for $T$. Each term is $Q_\tau$ times a matrix of Frobenius norm $|\delta_\tau|$ times $Z_\tau^*$, and unitary factors preserve the Frobenius norm, so the triangle inequality gives the bounds. $\square$

Lemma 3 (Hessenberg-triangular reduction). After line L1, $T$ is upper triangular and $H$ is upper Hessenberg.

Proof. Phase 1 processes columns $j = 1, \dots, n - 1$ and, for $i = n, \dots, j + 1$, rotates rows $i - 1, i$ to zero $t_{ij}$ (Lemma 1a). Rows $i - 1, i \ge j$ have zeros in every column $j' < j$ already, so these zeros stay; entries $t_{i'j}$ with $i' > i$ were zeroed before and are not touched. So at the end $T$ is upper triangular. Phase 2 processes columns $j = 1, \dots, n - 2$ and, for $i = n, \dots, j + 2$: (i) rotates rows $i - 1, i$ to zero $h_{ij}$; in columns $j' < j$ both rows are zero (rows $i - 1 \ge j + 1 > j' + 1$, and those columns are already Hessenberg), so they stay zero; in $T$, row $i - 1$ is zero before column $i - 1$ and row $i$ before column $i$, so the only new nonzero is $t_{i,i-1}$. (ii) Rotates columns $i - 1, i$ to zero $t_{i,i-1}$ (Lemma 1b); rows below $i$ of $T$ are zero in both columns, so $T$ is upper triangular again; in $H$ the columns $i - 1, i \ge j + 1$ are not yet processed and column $j$ is untouched, so $h_{ij} = 0$ and all earlier zeros of $H$ stay. At the end $H$ is upper Hessenberg. $\square$

Lemma 4 (the window). Suppose $H$ is Hessenberg, $T$ triangular, $h_{lo,lo-1} = 0$ (or $lo = 1$) and $h_{hi+1,hi} = 0$ (or $hi = n$). Then any rotation of two rows or two columns with indices in $[lo, hi]$ leaves $h_{lo,lo-1}$, $h_{hi+1,hi}$ and every entry of $H$ and $T$ in rows $> hi$ and columns $< lo$ unchanged and zero.

Proof. A row rotation in $[lo, hi]$ changes only rows in $[lo, hi]$; in columns $< lo$ these rows of $H$ and $T$ are zero (Hessenberg with $h_{lo,lo-1} = 0$, and triangular), and a rotation of two zeros gives zeros. A column rotation in $[lo, hi]$ changes only columns in $[lo, hi]$; in rows $> hi$ these columns are zero, by $h_{hi+1,hi} = 0$ and the structure. $\square$

Lemma 5 (the QZ step keeps the structure). If before line L9 $H$ is Hessenberg, $T$ is triangular with nonzero diagonal on $[lo, hi]$, and $h_{lo,lo-1} = 0$, $h_{hi+1,hi} = 0$ (where they exist), then after line L10 the same holds, with the same zero subdiagonal entries outside $(lo, hi]$.

Proof. By Lemma 4 the entries outside the window keep their zeros. Line L9 rotates rows $lo, lo + 1$. In $H$ both rows are zero before column $lo$ ($h_{lo,lo-1} = 0$), so $H$ stays Hessenberg; in $T$ the only new nonzero is $t_{lo+1,lo}$. We claim: before iteration $k$ of line L10, $H$ is Hessenberg and $T$ is triangular except for the entry $(k + 1, k)$. This holds for $k = lo$. In iteration $k$, the column rotation of columns $k, k + 1$ zeroes $t_{k+1,k}$ (Lemma 1b); the rows below $k + 1$ of $T$ are zero in both columns, so $T$ is triangular. In $H$, the rows below $k + 2$ are zero in both columns, so the only possible new nonzero is $h_{k+2,k}$, and only if $k + 2 \le hi$ (otherwise $h_{k+2,k+1} = h_{hi+1,hi} = 0$ and $h_{k+2,k}$ stays zero). If $k + 2 \le hi$, the row rotation of rows $k + 1, k + 2$ zeroes $h_{k+2,k}$ (Lemma 1a); in columns $< k$ both rows of $H$ are zero, so $H$ is Hessenberg again; in $T$, both rows are zero in columns $\le k$ (because $t_{k+1,k}$ was just zeroed), so the only new nonzero is $t_{k+2,k+1}$. This is the claim for $k + 1$. After the last iteration $k = hi - 1$, $T$ is triangular and $H$ Hessenberg. $\square$

Lemma 6 (deflation of an infinite eigenvalue). If at line L7 $t_{kk} = 0$ for some $k \in [lo, hi]$, then after the chase and the final rotation of rows $lo, lo + 1$, $H$ is Hessenberg, $T$ is triangular, $t_{lo,lo} = 0$ and $h_{lo+1,lo} = 0$.

Proof. We show by downward induction on $j = k, k - 1, \dots, lo$ that $H$ is Hessenberg, $T$ triangular and $t_{jj} = 0$. For $j = k$ this is the hypothesis. Going from $j$ to $j - 1$: the column rotation of columns $j - 1, j$ maps $(t_{j-1,j-1}, t_{j-1,j})$ to $(0, r)$ (Lemma 1b). Row $j$ of $T$ is zero in both columns ($t_{j,j-1} = 0$ and $t_{jj} = 0$), and so are all rows below; so $T$ stays triangular, $t_{jj}$ stays $0$, and now $t_{j-1,j-1} = 0$. In $H$, the only possible new nonzero is $h_{j+1,j-1}$, when $j + 1 \le hi$. Then the row rotation of rows $j, j + 1$ zeroes it (Lemma 1a); in columns $< j - 1$ both rows of $H$ are zero, so $H$ is Hessenberg; in $T$, both rows are zero in columns $\le j$ (as $t_{jj} = t_{j+1,j} = 0$), so $T$ stays triangular. At $j = lo$, the rotation of rows $lo, lo + 1$ zeroes $h_{lo+1,lo}$ (Lemma 1a). In $H$, both rows are zero before column $lo$, so $H$ stays Hessenberg; in $T$, both rows are zero in column $lo$ ($t_{lo,lo} = 0$, $t_{lo+1,lo} = 0$) and before it, so $T$ stays triangular and $t_{lo,lo}$ stays $0$. $\square$

Theorem 7 (correctness and termination). The algorithm terminates. When it returns with converged true, $S = H$ and $T$ are upper triangular, $Q$ and $Z$ are unitary, and $Q^* (A + E) Z = S$, $Q^* (B + F) Z = T$ with $\|E\|_F \le 2 (n - 1) \varepsilon \max_\tau \|H_\tau\|_F$ and $\|F\|_F \le n (n - 1) \varepsilon \|B\|_F$, where $H_\tau$ are the intermediate matrices, so $\|E\|_F \le 2 (n - 1) \varepsilon \|A\|_F (1 + O(n \varepsilon))$. The generalized eigenvalues of the pencil $(A + E, B + F)$ are the pairs $(s_{ii}, t_{ii})$.

Proof. Structure. Lemma 3 gives the structure after line L1. Lines L4 and L6 only set entries to zero, which keeps it. Lemma 6 keeps it at line L7, and Lemma 5 at lines L9 and L10, whose hypotheses hold because line L7 is reached only after all $t_{kk}$, $lo \le k \le hi$, have been tested, and line L9 only when they are all nonzero.

Zero subdiagonal entries stay zero. All rotations act on indices inside a window $[lo, hi]$ bounded by zero subdiagonal entries or the ends, so by Lemma 4 a zero $h_{k,k-1}$ is never changed again. Inside the window all subdiagonal entries are nonzero by the choice of $lo$.

Termination. Let $N$ be the number of nonzero entries $h_{k,k-1}$ with $2 \le k \le hi$, and consider $\Phi = N + hi$. Line L4 does not increase $\Phi$. Line L5 decreases $hi$ by one, and the entry $h_{hi,hi-1}$ it drops is zero, so $\Phi$ drops by one. Line L7 makes the nonzero $h_{lo+1,lo}$ zero (it is nonzero by the choice of $lo$, and $lo + 1 \le hi$ because $h_{hi,hi-1} \ne 0$), so $N$ and $\Phi$ drop by one. A QZ step leaves $\Phi$ unchanged and increases the step counter, which is capped. So the number of passes through the loop is at most $2n$ plus the cap, and each pass is finite.

Triangular result. If the loop ends with converged true, then $hi = 1$, and $hi$ was decreased from each value $m$ only when $h_{m,m-1} = 0$; these zeros persist, so $H$ is upper triangular. $T$ is upper triangular by the structure argument.

Backward error. By Lemma 2, $\|E\|_F \le \sum_\tau |\delta_\tau|$. Each subdiagonal position is set to zero at most once (it stays zero afterwards and the test only acts on nonzero entries), so there are at most $n - 1$ terms, and each satisfies $|\delta_\tau| \le \varepsilon (|h_{k-1,k-1}| + |h_{kk}|) \le 2 \varepsilon \|H_\tau\|_F$, or $|\delta_\tau| \le \varepsilon \|A\|_F$ when both diagonal entries are zero. Each diagonal entry of $T$ set to zero satisfies $|\gamma_\sigma| \le \varepsilon \|B\|_F$, and the test of line L6 is applied to at most $n$ entries in each pass that leads to line L7; by the termination argument there are at most $n - 1$ such passes, so at most $n (n - 1)$ terms. Finally $\|H_\tau\|_F = \|A + E_\tau\|_F$ by unitary invariance, which gives the first-order bound. The last statement holds because $\det((A + E) - \lambda (B + F)) = \det(Q) \det(S - \lambda T) \det(Z^*) = \det(Q Z^*) \prod_i (s_{ii} - \lambda t_{ii})$. $\square$

Remark (what is cited, not proved). That one QZ step equals one step of the shifted QR algorithm on $H T^{-1}$ and on $T^{-1} H$ (when $T$ is nonsingular) follows from the implicit Q theorem; see GVL Section 7.7 and Watkins and Elsner (1994). It gives the quadratic local convergence of $h_{hi,hi-1}$ for a simple eigenvalue with the shift of line L8; the tests measure it. There is no proof that every pencil converges within the cap, which is why the implementation reports converged. Backward stability in floating point, $\|E\|_2 = O(u \|A\|_2)$ and $\|F\|_2 = O(u \|B\|_2)$ including rounding, is due to Moler and Stewart (1973); the tests observe it.