Proofs: MRRR algorithm

MRRR assigns every eigenvalue one vector, its twisted vectors have residuals bounded by the eigenvalue error, and the representation tree is finite full

This proof shows: correctness, termination.

Exact arithmetic throughout, except where an approximation $\hat\lambda$ of an eigenvalue is fed in; the floating-point behaviour is cited at the end. Notation of the page: $T$ is real symmetric tridiagonal of order $n$, unreduced on each block; a representation is $L D L^T$ with $L$ unit lower bidiagonal ($l_i$ in position $(i + 1, i)$) and $D = \operatorname{diag}(d_1, \dots, d_n)$; its tridiagonal has diagonal $\alpha_i = d_i + d_{i-1} l_{i-1}^2$ (with $d_0 l_0^2 = 0$) and off-diagonal $\beta_i = d_i l_i$. We assume that no pivot $D^+_i$ or $D^-_i$ below is zero (the code replaces an exact zero by a tiny number, which is a perturbation of size $2^{-1000}$).

Lemma 1 (stationary qd). With $s_1 = -\tau$, $D^+_i = d_i + s_i$, $L^+_i = d_i l_i / D^+_i$, $s_{i+1} = L^+_i l_i s_i - \tau$ and $D^+_n = d_n + s_n$, we have $L^+ D^+ (L^+)^T = L D L^T - \tau I$.

Proof. Both sides are symmetric tridiagonal. Off-diagonal: $D^+_i L^+_i = d_i l_i = \beta_i$. Diagonal $1$: $D^+_1 = d_1 - \tau = \alpha_1 - \tau$. Diagonal $i + 1$: we need $D^+_{i+1} + D^+_i (L^+_i)^2 = d_{i+1} + d_i l_i^2 - \tau$. Now $d_i l_i^2 - D^+_i (L^+_i)^2 = d_i l_i^2 - d_i l_i L^+_i = l_i (D^+_i L^+_i - d_i L^+_i) = l_i L^+_i (D^+_i - d_i) = L^+_i l_i s_i$, so the requirement is $D^+_{i+1} = d_{i+1} + L^+_i l_i s_i - \tau = d_{i+1} + s_{i+1}$, which is the definition. $\square$

Lemma 2 (progressive qd from the bottom). With $p_n = d_n - \mu$ and, for $i = n - 1, \dots, 1$, $D^-_{i+1} = d_i l_i^2 + p_{i+1}$, $U^-_i = d_i l_i / D^-_{i+1}$, $p_i = p_{i+1} d_i / D^-_{i+1} - \mu$, and $D^-_1 = p_1$, we have $U^- D^- (U^-)^T = L D L^T - \mu I$ with $U^-$ unit upper bidiagonal ($U^-_i$ in position $(i, i + 1)$).

Proof. The product has off-diagonal $D^-_{i+1} U^-_i = d_i l_i$ and diagonal $D^-_i + D^-_{i+1} (U^-_i)^2$, with the second term absent for $i = n$. For $i = n$: $D^-_n = d_{n-1} l_{n-1}^2 + d_n - \mu = \alpha_n - \mu$. For $i < n$, using $D^-_i = d_{i-1} l_{i-1}^2 + p_i$ (or $p_1$ for $i = 1$), the requirement $D^-_i + (d_i l_i)^2 / D^-_{i+1} = \alpha_i - \mu$ becomes $p_i = d_i - \mu - d_i^2 l_i^2 / D^-_{i+1} = d_i (D^-_{i+1} - d_i l_i^2) / D^-_{i+1} - \mu = d_i p_{i+1} / D^-_{i+1} - \mu$, the definition. $\square$

Lemma 3 (twisted factorization). Let $M = L D L^T - \mu I$, $\gamma_k = s_k + p_k + \mu$, and fix $r$. Define $z_r = 1$, $z_i = -L^+_i z_{i+1}$ for $i = r - 1, \dots, 1$, and $z_{i+1} = -U^-_i z_i$ for $i = r, \dots, n - 1$. Then $M z = \gamma_r e_r$. If $M$ is nonsingular, $\gamma_r = 1 / (M^{-1})_{rr}$.

Proof. By Lemma 1, $M z = L^+ D^+ w$ with $w = (L^+)^T z$, and $w_i = z_i + L^+_i z_{i+1} = 0$ for $i < r$. Row $i$ of $L^+ D^+ w$ is $D^+_i w_i + L^+_{i-1} D^+_{i-1} w_{i-1}$, which vanishes for $i < r$. By Lemma 2, $M z = U^- D^- v$ with $v = (U^-)^T z$, and $v_i = z_i + U^-_{i-1} z_{i-1} = 0$ for $i > r$; row $i$ of $U^- D^- v$ is $D^-_i v_i + U^-_i D^-_{i+1} v_{i+1}$, which vanishes for $i > r$. So $M z = c \, e_r$ with $c = (Mz)_r = \beta_{r-1} z_{r-1} + (\alpha_r - \mu) + \beta_r z_{r+1}$. Now $\beta_{r-1} z_{r-1} = -D^+_{r-1} (L^+_{r-1})^2 = D^+_r - (\alpha_r - \mu)$ by the diagonal of Lemma 1, and $\beta_r z_{r+1} = -D^-_{r+1} (U^-_r)^2 = D^-_r - (\alpha_r - \mu)$ by the diagonal of Lemma 2 (terms with index $0$ or $n + 1$ are absent, and the identities still hold with $D^+_1 = \alpha_1 - \mu$, $D^-_n = \alpha_n - \mu$). Hence $c = D^+_r + D^-_r - (\alpha_r - \mu) = (d_r + s_r) + (d_{r-1} l_{r-1}^2 + p_r) - (d_r + d_{r-1} l_{r-1}^2 - \mu) = \gamma_r$. Finally $z = \gamma_r M^{-1} e_r$, and $z_r = 1$ gives $\gamma_r (M^{-1})_{rr} = 1$. $\square$

Lemma 4 (counts). The number of negative $D^+_i$ in Lemma 1 equals the number of eigenvalues of $L D L^T$ less than $\tau$.

Proof. $L D L^T - \tau I = L^+ D^+ (L^+)^T$ is congruent to $D^+$, and by Sylvester's law of inertia it has as many negative eigenvalues as $D^+$ has negative entries. Its eigenvalues are $\lambda - \tau$. $\square$

Lemma 5 (bisection). Let $\nu(x)$ be the count of Lemma 4. The routine bisect(k) ends with $lo < hi$, $\nu(lo) < k \le \nu(hi)$, and returns the midpoint $\hat\lambda$; hence the $k$-th smallest eigenvalue $\lambda_k$ of $L D L^T$ satisfies $|\hat\lambda - \lambda_k| \le (hi - lo)/2 \le u \max(|lo|, |hi|)$ when the width test stops it, and the loop runs at most $400$ steps after finitely many bracket expansions.

Proof. $\nu$ is nondecreasing, $0$ below the spectrum and $n$ above it, so moving $lo$ down by doubling steps reaches $\nu(lo) < k$ and moving $hi$ up reaches $\nu(hi) \ge k$ after finitely many steps. Each bisection step keeps the half on which the invariant holds. The invariant gives $\lambda_k \in [lo, hi)$: at most $k - 1$ eigenvalues lie below $lo$, so $\lambda_k \ge lo$, and at least $k$ lie below $hi$, so $\lambda_k < hi$. The midpoint is then within $(hi - lo)/2$ of $\lambda_k$. $\square$

Lemma 6 (residual and angle of a twisted vector). Let $\lambda_j$ be the eigenvalue of $L D L^T$ nearest to $\mu$, $\delta = |\lambda_j - \mu|$, $g = \min_{i \ne j} |\lambda_i - \mu|$, and suppose $\delta < g / (2n)$. For $r$ minimizing $|\gamma_k|$ and $\hat z = z / \|z\|_2$: $\|M \hat z\|_2 \le |\gamma_r| \le 2 n \delta$, and $\sin\angle(\hat z, v_j) \le 2 n \delta / g$.

Proof. With orthonormal eigenvectors $v_i$, $(M^{-1})_{kk} = \sum_i v_i(k)^2 / (\lambda_i - \mu)$. Pick $k$ with $v_j(k)^2 \ge 1/n$. Then $|(M^{-1})_{kk}| \ge 1/(n\delta) - 1/g \ge 1/(2 n \delta)$, so by Lemma 3 $|\gamma_k| \le 2n\delta$ and the minimum is no larger. By Lemma 3, $\|M \hat z\| = |\gamma_r| / \|z\| \le |\gamma_r|$, since $\|z\| \ge |z_r| = 1$. Write $\hat z = \cos\theta \, v_j + \sin\theta \, w$ with $w \perp v_j$ a unit vector; then $\|M \hat z\| \ge \|M w\| \sin\theta \ge g \sin\theta$, because $M$ maps the orthogonal complement of $v_j$ into itself with smallest singular value $g$. $\square$

Theorem 7 (correctness and termination). On each unreduced block: (a) the loop terminates; (b) every eigenvalue index is assigned exactly one value and one unit vector; (c) every representation on the stack equals $T - \sigma I$ for its accumulated shift $\sigma$; (d) for an index finished as a singleton in a representation where the relative gap condition holds, the returned pair $(\sigma + \hat\lambda, \hat z)$ satisfies $\|(T - (\sigma + \hat\lambda) I) \hat z\|_2 \le 2 n \delta$ and $\sin\angle(\hat z, v) \le 2 n \delta / g$ with $\delta$, $g$ as in Lemma 6 for that representation; (e) vectors from a fallback cluster are orthonormal.

Proof. (c) The root satisfies $L D L^T = T - \sigma I$ by Gaussian elimination without pivoting, whose pivots are all positive by the test in root_ldl (so $T - \sigma I$ is positive definite). A child built from $(L D L^T, \tau)$ satisfies $L^+ D^+ (L^+)^T = L D L^T - \tau I = T - (\sigma + \tau) I$ by Lemma 1, and is pushed with shift $\sigma + \tau$. Induction on the depth.

(b) The root is pushed with all indices $1, \dots, m$. When a node is popped, _groups partitions its index list into consecutive runs; each run of length one is assigned at once, each longer run is either assigned at once (fallback) or pushed as a child with exactly that index list. So the index lists of the nodes ever pushed form a tree in which the children of a node partition a subset of its indices and the remaining indices are assigned in the node; every index is assigned exactly once. The eigenvalues of a child are the same indices because a shift preserves the order: the $k$-th smallest eigenvalue of $T - \sigma I$ minus $\tau$ is the $k$-th smallest of $T - (\sigma + \tau) I$.

(a) A child is pushed only when the depth of its parent is below maxdepth, so the tree has height at most maxdepth, and each node has at most $m/2$ children because each child has at least two indices. So the tree is finite. Each node does finitely many operations: by Lemma 5 each bisection ends, the shift search tries at most $2 \cdot \mathrm{TRIES}$ shifts, and the fallback tries at most $m$ twist indices per eigenvalue.

(d) By (c), $T - (\sigma + \hat\lambda) I = L D L^T - \hat\lambda I = M$, and Lemma 6 applies.

(e) Each kept vector is normalized after being made orthogonal to the earlier ones (twice, which in exact arithmetic is once); the Rayleigh-Ritz step multiplies this orthonormal set by the orthogonal eigenvector matrix of a small symmetric matrix. $\square$

Remarks (what is cited, not proved). (1) The size of $\delta$. By Lemma 5, $\delta \le u |\lambda_j^{(\mathrm{rep})}|$ for an eigenvalue $\lambda_j^{(\mathrm{rep})}$ of the representation, so Lemma 6 gives $\sin\angle \le 2 n u / \mathrm{relgap}$ with $\mathrm{relgap} = g / |\lambda_j^{(\mathrm{rep})}|$; the clustering rule makes $\mathrm{relgap}$ at least about $\mathrm{gaptol}$ for singletons. This is why MRRR shifts a cluster until its members have large relative gaps. (2) Floating point. That the computed stationary and progressive transforms are exact for small relative perturbations of $L$, $D$ and the outputs (mixed relative stability), that a relatively robust representation determines its eigenvalues and eigenvectors to high relative accuracy, and that the computed vectors from different representations are numerically orthogonal, $|\hat z_i^T \hat z_j| = O(n u / \mathrm{gaptol})$, and have residuals $O(n u \|T\|)$, are the theorems of Dhillon and Parlett (2004a, 2004b) and Dhillon, Parlett and Vömel (2006), under the hypothesis that every representation used is relatively robust for its cluster. The element growth test of the code is the heuristic check of that hypothesis used in LAPACK; it is not a proof. The tests observe residuals below $1.3 \, n u \|T\|$, orthogonality below $9 \, n u$ on random matrices and on glued Wilkinson matrices, and $68 \, n u$ on $W_{21}^+$, whose singletons have relative gaps just above $\mathrm{gaptol}$.