The algebraic Riccati equation, by the Schur method

The continuous-time algebraic Riccati equation

A^T P + P A - P B R^-1 B^T P + W = 0

is quadratic in P, and the thing you want is not just any solution but the one that stabilizes: the symmetric P for which A - B R^-1 B^T P has every eigenvalue in the left half plane. This module computes it the way SLICOT's SB02MD does, which is Laub's method of 1979.

The idea is that the quadratic equation is a linear one in disguise. Build the 2n by 2n Hamiltonian

H = [[A, -B R^-1 B^T], [-W, -A^T]]

and the stabilizing P is read off its stable invariant subspace. If the columns of [X1; X2] are a basis for the eigenvectors of H belonging to the eigenvalues with negative real part, then P = X2 X1^-1. A Hamiltonian's spectrum is symmetric about the imaginary axis, so exactly n of its 2n eigenvalues are stable, and the subspace has exactly the right dimension.

Nothing here forms an eigenvector. lib/schur reduces H to real Schur form, reorders the stable blocks to the top left, and the first n columns of the orthogonal factor are an orthonormal basis for that subspace — better conditioned than eigenvectors, and available whether or not H is diagonalizable.

are_sign solves the same equation by a completely different route, Roberts's matrix sign function, and the two are checked against each other. Two unrelated algorithms landing on the same matrix is a much stronger argument than either one looking right on its own.

They do not land on it equally well, and the test below measures the gap. Over the systems it tries, the Schur method's relative residual stays under 4 × 10⁻¹⁵ while the sign function's reaches 7 × 10⁻¹³, and where the two answers differ it is the sign function that has moved. Its last step solves an overdetermined system through its normal equations, which squares the condition number; that is the price of not computing a Schur form, and it is why SB02MD is the Schur one.

Used with import are.

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 46–488

from schur import reorder_stable, schur
from linalg import inverse, test_matrix

ARE_SIGN_STEPS = 60

#| `B R^-1 B^T`, formed by solving `R Y = B^T` rather than inverting `R`. `B` is
#| `n` by `m` and `R` is `m` by `m`.
def are_gram(B: Arr, R: Arr) {
    n = len(B)
    m = len(R)
    Y = zeros(m, n)
    rhs = zeros(m)
    col = zeros(m)
    for j to n {
        for i to m {
            rhs[i] = B[j, i]
        }
        lin_solve(R, rhs, col)
        for i to m {
            Y[i, j] = col[i]
        }
    }
    return B @ Y
}

#| The Hamiltonian whose stable invariant subspace carries the answer.
def are_hamiltonian(A: Arr, B: Arr, R: Arr, W: Arr) {
    n = len(A)
    #: H = [[A, -G], [-W, -A^T]]
    G = are_gram(B, R)
    H = zeros(2 * n, 2 * n)
    for i to n {
        for j to n {
            H[i, j] = A[i, j]
            H[i, n + j] = -G[i, j]
            H[n + i, j] = -W[i, j]
            H[n + i, n + j] = -A[j, i]
        }
    }
    return H
}

#| Solve the Riccati equation into `P`. Returns the dimension of the stable
#| subspace that was found, which had better be `n`: anything else means the
#| Hamiltonian has an eigenvalue on the imaginary axis, and then there is no
#| stabilizing solution to find.
def are_solve(A: Arr, B: Arr, R: Arr, W: Arr, P: Arr) {
    n = len(A)
    #: stable eigenvalues first; a cut other than n means no stabilizing solution
    H = are_hamiltonian(A, B, R, W)
    T = zeros(2 * n, 2 * n)
    Z = zeros(2 * n, 2 * n)
    junk = schur(H, T, Z)
    cut = reorder_stable(T, Z)
    if cut != n {
        return cut
    }
    # the first n columns of Z are the subspace; split them and solve P X1 = X2
    X1t = zeros(n, n)
    rhs = zeros(n)

from schur import reorder_stable, schur
from linalg import inverse, test_matrix

ARE_SIGN_STEPS = 60

#| `B R^-1 B^T`, formed by solving `R Y = B^T` rather than inverting `R`. `B` is
#| `n` by `m` and `R` is `m` by `m`.
def are_gram(B: Arr, R: Arr) {
    n = len(B)
    m = len(R)
    Y = zeros(m, n)
    rhs = zeros(m)
    col = zeros(m)
    for j to n {
        for i to m {
            rhs[i] = B[j, i]
        }
        lin_solve(R, rhs, col)
        for i to m {
            Y[i, j] = col[i]
        }
    }
    return B @ Y
}

#| The Hamiltonian whose stable invariant subspace carries the answer.
def are_hamiltonian(A: Arr, B: Arr, R: Arr, W: Arr) {
    n = len(A)
    #: H = [[A, -G], [-W, -A^T]]
    G = are_gram(B, R)
    H = zeros(2 * n, 2 * n)
    for i to n {
        for j to n {
            H[i, j] = A[i, j]
            H[i, n + j] = -G[i, j]
            H[n + i, j] = -W[i, j]
            H[n + i, n + j] = -A[j, i]
        }
    }
    return H
}

#| Solve the Riccati equation into `P`. Returns the dimension of the stable
#| subspace that was found, which had better be `n`: anything else means the
#| Hamiltonian has an eigenvalue on the imaginary axis, and then there is no
#| stabilizing solution to find.
def are_solve(A: Arr, B: Arr, R: Arr, W: Arr, P: Arr) {
    n = len(A)
    #: stable eigenvalues first; a cut other than n means no stabilizing solution
    H = are_hamiltonian(A, B, R, W)
    T = zeros(2 * n, 2 * n)
    Z = zeros(2 * n, 2 * n)
    junk = schur(H, T, Z)
    cut = reorder_stable(T, Z)
    if cut != n {
        return cut
    }
    # the first n columns of Z are the subspace; split them and solve P X1 = X2
    X1t = zeros(n, n)
    rhs = zeros(n)