Notebooks · lib/schur

The real Schur decomposition, by Householder and Francis

Every real square matrix can be written A = Q T Q^T with Q orthogonal and T quasi-triangular: upper triangular except for 2 by 2 blocks on the diagonal, one for each complex conjugate pair of eigenvalues. This is the decomposition that control theory runs on. SLICOT's Riccati, Lyapunov and Sylvester solvers all begin by computing it, because it is the only way anyone knows to get at an invariant subspace without inverting something that may be nearly singular.

It is computed in two stages, and this module does both.

hessenberg reduces A to upper Hessenberg form — zero below the subdiagonal -- by n - 2 Householder reflections. This is a finite computation; it ends.

schur then runs the implicit double-shift QR iteration, which is Francis's, on the Hessenberg matrix. Each sweep introduces a bulge at the top left and chases it off the bottom right, and the subdiagonal entries go to zero one after another. This does not end in a fixed number of steps — there is no finite algorithm for eigenvalues of order five and above — so it iterates until each subdiagonal entry is negligible against its neighbors, and deflates.

Every transformation here is orthogonal, applied on both sides. That is what makes the whole thing backward stable, and it is why Q is carried along rather than reconstructed at the end.

Each reflector is applied only to the rows and columns that can be nonzero. That is not for speed. Reach outside the window and the reflector writes a little dust back into the subdiagonal entries deflation has already set to zero, and the iteration never finishes: the residual stays perfect, because a similarity is a similarity, while the factor is quietly not quasi-triangular.

Used with import schur.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 37–817

from linalg import test_matrix

SCHUR_EPS = 2.220446049250313e-16

# A double-shift bulge is three rows tall: the first column of (H - s1)(H - s2)
# has three nonzeros below a Hessenberg matrix's first row.
SCHUR_BULGE = 3

#| Apply the Householder reflector built from `v` to the window of `m` rows and
#| columns starting at `p`: `H` from both sides, so its eigenvalues do not move,
#| and `Q` from the right, so that `A = Q H Q^T` still holds afterwards.
#|
#| The two ranges are not an optimization. From the left it touches columns `c0`
#| and rightwards, and from the right it touches rows `r1` and upwards, because
#| everything outside those is already zero and must stay zero -- a reflector
#| applied to the whole row would write dust back into a subdiagonal entry that
#| deflation had just decided was negligible, and the factor would never finish
#| being quasi-triangular.
def reflect(H, Q, v, p, m, c0, r1) {
    n = len(H)
    w = v[0:m]
    vn = dot(w, w)
    if vn == 0.0 {
        return 0.0
    }
    #: only columns c0 on: the rest are zero and must stay zero
    s = 2.0 * (w @ H[p:p + m, c0:n]) / vn
    H[p:p + m, c0:n] = H[p:p + m, c0:n] - outer(w, s)
    s = 2.0 * (H[0:r1 + 1, p:p + m] @ w) / vn
    H[0:r1 + 1, p:p + m] = H[0:r1 + 1, p:p + m] - outer(s, w)
    s = 2.0 * (Q[0:n, p:p + m] @ w) / vn
    Q[0:n, p:p + m] = Q[0:n, p:p + m] - outer(s, w)
    return 0.0
}

#| Apply the plane rotation `[[cs, -sn], [sn, cs]]` on rows and columns `k` and
#| `k + 1` the same way: `T` from both sides, `Q` from the right.
#|
#| This is only ever called on a block that has already deflated, so everything
#| left of column `k` in those two rows, and everything below row `k + 1` in
#| those two columns, is zero. It stays out of both. A sweep leaves dust of
#| order 1e-17 below the subdiagonal where exact arithmetic leaves zeros, and a rotation
#| over the full row would mix that dust straight into the subdiagonal entry
#| deflation had just set to zero, turning a converged 2 by 2 into two blocks
#| that no longer look converged.
def rotate(T, Q, k, cs, sn) {
    n = len(T)
    for jj to n - k {
        j = k + jj
        t0 = T[k, j]
        t1 = T[k + 1, j]
        T[k, j] = cs * t0 + sn * t1
        T[k + 1, j] = cs * t1 - sn * t0
    }
    for i to k + 2 {
        t0 = T[i, k]
        t1 = T[i, k + 1]
        T[i, k] = cs * t0 + sn * t1
        T[i, k + 1] = cs * t1 - sn * t0
    }

from linalg import test_matrix

SCHUR_EPS = 2.220446049250313e-16

# A double-shift bulge is three rows tall: the first column of (H - s1)(H - s2)
# has three nonzeros below a Hessenberg matrix's first row.
SCHUR_BULGE = 3

#| Apply the Householder reflector built from `v` to the window of `m` rows and
#| columns starting at `p`: `H` from both sides, so its eigenvalues do not move,
#| and `Q` from the right, so that `A = Q H Q^T` still holds afterwards.
#|
#| The two ranges are not an optimization. From the left it touches columns `c0`
#| and rightwards, and from the right it touches rows `r1` and upwards, because
#| everything outside those is already zero and must stay zero -- a reflector
#| applied to the whole row would write dust back into a subdiagonal entry that
#| deflation had just decided was negligible, and the factor would never finish
#| being quasi-triangular.
def reflect(H, Q, v, p, m, c0, r1) {
    n = len(H)
    w = v[0:m]
    vn = dot(w, w)
    if vn == 0.0 {
        return 0.0
    }
    #: only columns c0 on: the rest are zero and must stay zero
    s = 2.0 * (w @ H[p:p + m, c0:n]) / vn
    H[p:p + m, c0:n] = H[p:p + m, c0:n] - outer(w, s)
    s = 2.0 * (H[0:r1 + 1, p:p + m] @ w) / vn
    H[0:r1 + 1, p:p + m] = H[0:r1 + 1, p:p + m] - outer(s, w)
    s = 2.0 * (Q[0:n, p:p + m] @ w) / vn
    Q[0:n, p:p + m] = Q[0:n, p:p + m] - outer(s, w)
    return 0.0
}

#| Apply the plane rotation `[[cs, -sn], [sn, cs]]` on rows and columns `k` and
#| `k + 1` the same way: `T` from both sides, `Q` from the right.
#|
#| This is only ever called on a block that has already deflated, so everything
#| left of column `k` in those two rows, and everything below row `k + 1` in
#| those two columns, is zero. It stays out of both. A sweep leaves dust of
#| order 1e-17 below the subdiagonal where exact arithmetic leaves zeros, and a rotation
#| over the full row would mix that dust straight into the subdiagonal entry
#| deflation had just set to zero, turning a converged 2 by 2 into two blocks
#| that no longer look converged.
def rotate(T, Q, k, cs, sn) {
    n = len(T)
    for jj to n - k {
        j = k + jj
        t0 = T[k, j]
        t1 = T[k + 1, j]
        T[k, j] = cs * t0 + sn * t1
        T[k + 1, j] = cs * t1 - sn * t0
    }
    for i to k + 2 {
        t0 = T[i, k]
        t1 = T[i, k + 1]
        T[i, k] = cs * t0 + sn * t1
        T[i, k + 1] = cs * t1 - sn * t0
    }