Solving the Riccati equation without integrating it

Adam Zettel · September 23, 2026 · written with an agent

A^{⊤} P + P A - P B R^{- 1} B^{⊤} P + W = 0

Every LQR controller comes from a quadratic matrix equation, and the tempting way to solve it is to integrate the differential version backwards until it stops moving. On a one-state system that method returns 73.783 where the answer is 414.214, with no warning. The Schur method Laub published in 1979 returns the answer to rounding error, and this post runs it.

Open in the workbench Download .vivOutputs below are from the site build’s run; Run all runs it again in your browser.

Look inside this program

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

The gains of a linear-quadratic regulator come out of

A^{⊤} P + P A - P B R^{- 1} B^{⊤} P + W = 0,

quadratic in the matrix $P$ . There is an obvious way to solve it: run $P ˙ = - (A^{⊤} P + P A - P B R^{- 1} B^{⊤} P + W)$ backwards from zero until it stops changing. Each step is one matrix multiply and one matrix add, it converges on most problems, and the earlier LQR post on this site uses it.

Take $A = - 0.001$ , $B = 0.001$ , $W = 1$ , $R = 1$ , one state. The quadratic formula gives $P = 414.214$ exactly. The closed loop it produces decays like $e^{- 0.00141 t}$ , so the Riccati flow settles over roughly seven hundred seconds, and eighty seconds of backward integration gets nowhere near. It returns 73.783, and nothing in its output signals the error.

The fix is to solve a linear eigenvalue problem instead. Stack $A$ and the two weights into the $2 n \times 2 n$ Hamiltonian

H = (A - W - B R^{- 1} B^{⊤} - A^{⊤}),

whose spectrum is symmetric about the imaginary axis: every eigenvalue $λ$ is matched by $- λ$ , so exactly $n$ of the $2 n$ are stable. Take an orthonormal basis $[X_{1}; X_{2}]$ for the invariant subspace of those $n$ , and $P = X_{2} X_{1}^{- 1}$ . No quadratic equation is solved.

Computing that basis is the main work. It is done by reducing $H$ to real Schur form (block upper triangular) and moving the stable blocks to the top left. The animation below shows that reduction: the eleven subdiagonal entries of the Hessenberg matrix, driven to zero by Francis's implicit double-shift QR iteration, one turn of its outer loop at a time. Entry 9 (counting from 0) goes first, and fast. Its size over four successive turns is $3.5 \times 1 0^{- 3}$ , $7.7 \times 1 0^{- 7}$ , $1.9 \times 1 0^{- 13}$ , $3.8 \times 1 0^{- 27}$ : the exponent doubles each time, which is the quadratic convergence a good shift gives. When an entry reaches rounding level, the matrix splits there into two smaller problems and the iteration moves on to the block above. After thirty turns nothing is left to do, and the four entries still nonzero, at $k = 0, 2, 4, 6$ , belong to the 2 by 2 blocks holding the Hamiltonian's four complex-conjugate pairs of eigenvalues.

The last figure compares residuals on ten systems. The Schur method's relative residual stays below $1 0^{- 14}$ on all of them, including those where the largest entry of $P$ runs to five figures and the absolute residual is $5 \times 1 0^{- 7}$ . The matrix sign function, another way to get $P$ from $H$ , is more than thirty times worse on at least one.

Printed

exact_slow: 414.214

schur_slow: [[414.214]]

naive_slow: [[73.783]]

worst_open: 2.0668

dim_stable: 6

0.00496361

8.03144e-16

-0.335031

ham_trace: 0

pmax: [5.7319, 3003.78, 150.669, 50.9556, 8122.15, 63.5384, 56.0792, 2228.53, 4011.09, 22253.8]

the six stable ones, which P is built fromthe six unstable ones, discarded

Each subdiagonal entry shrinking to zero

entry 9, the first to goentry 7, which goes nextentry 3, higher up

The subdiagonal being driven to zero

0.0 / 9.9 s

Relative residual, ten systems, two methods

the Schur methodthe matrix sign function

Tests 10 of 10 pass

the backward integration is wrong by four fifths on the slow system (passes)
the demo system starts unstable and the gains put every mode in the left half plane (passes)
the Hamiltonian has exactly six stable eigenvalues and zero trace (passes)
every eigenvalue of the Hamiltonian is matched by minus itself (passes)
the largest P, five figures, still has a relative residual below 1e-14 (passes)
P solves the Riccati equation with a relative residual below 1e-14 (passes)
thirty turns leaves four subdiagonal entries, one per complex pair (passes)
entry 9 converges quadratically: the exponent doubles each turn (passes)
the Schur method beats the sign function on every one of the ten systems (passes)
some of these systems are nearly uncontrollable, which is the hard case (passes)