BUG: solve_discrete_riccati_system returns NaN instead of raising when the iteration overflows

The bug is in the while error > tolerance loop of solve_discrete_riccati_system in quantecon/_matrix_eqn.py, where a NaN error silently exits the loop. Start by reading that loop and the fail_msg it would raise, then run the reproducer in the issue to confirm NaN is returned. Done means a non-finite error raises fail_msg and the new test in quantecon/tests/test_matrix_eqn.py passes while existing tests stay green.

Summary

solve_discrete_riccati_system loops while error > tolerance, where error sums the largest absolute change in each P(s). When the iterates overflow, error becomes nan. Because nan > tolerance is False, the loop exits as if it had converged and returns non-finite Ps without raising. LQMarkov.stationary_values then fails with SciPy's ValueError: array must not contain infs or NaNs instead of the solver's own failure message, which for beta >= 1 points to stabilizability. With the default max_iter=1000 this happens when P grows by more than about 2.03x per iteration. Slower divergence reaches max_iter and raises as intended.

What the code does

quantecon/_matrix_eqn.py#L294-L319:

    iteration = 0
    while error > tolerance:

        if iteration > max_iter:
            raise ValueError(fail_msg.format(max_iter))

        else:
            error = 0
            for i in range(m):
                ...
                Ps1[i][:, :] = Rs[i] + sum1 - sum2
                error += np.max(np.abs(Ps1[i] - Ps[i]))

            Ps[:, :, :] = Ps1[:, :, :]
            iteration += 1

    return Ps

Once the iterates overflow, the update produces NaN (in the reproducer below, from 0 * inf inside the matrix products), so error is nan. The loop then stops and returns Ps without reaching the raise. As a result fail_msg, which carries the stabilizability hint for beta >= 1, is never raised. LQMarkov.stationary_values builds F from the NaN Ps and calls scipy.linalg.solve at _lqcontrol.py#L545, which rejects the non-finite input.

Reproducer

from quantecon import LQMarkov
from quantecon import _matrix_eqn as qme  # not exported at the top level

# One Markov state, x' = a x + 0 u: the control cannot stabilize x.
# beta=1 (the LQMarkov default) selects the solver's beta>=1 failure message.
for a in (1.1, 2.0):
    lq = LQMarkov([[1.0]], [1.0], [1.0], [a], [0.0], beta=1)
    try:
        lq.stationary_values()
    except ValueError as e:
        print(f"a={a}: LQMarkov raised ValueError: {e}")

# Calling the solver directly for a=2.0: no exception, NaN returned
Ps = qme.solve_discrete_riccati_system(lq.Π, lq.As, lq.Bs, lq.Cs, lq.Qs,
                                       lq.Rs, lq.Ns, lq.beta)
print("solver returned Ps =", Ps.ravel())

Output on main at 70fa215 (NumPy 2.4.6, SciPy 1.18.0, Python 3.13; warning paths shortened):

a=1.1: LQMarkov raised ValueError: Convergence failed after 1000 iterations. When beta>=1, a stationary solution may not exist if the system cannot be stabilized without discounting. Check the model specification.
.../quantecon/_matrix_eqn.py:307: RuntimeWarning: overflow encountered in matmul
  sum1 += beta * Π[i, j] * As[i].T @ Ps[j] @ As[i]
.../quantecon/_matrix_eqn.py:309: RuntimeWarning: invalid value encountered in matmul
  (beta * As[i].T @ Ps[j] @ Bs[i] + Ns[i].T) @ \
.../quantecon/_matrix_eqn.py:310: RuntimeWarning: invalid value encountered in matmul
  solve(Qs[i] + beta * Bs[i].T @ Ps[j] @ Bs[i],
.../quantecon/_matrix_eqn.py:311: RuntimeWarning: invalid value encountered in matmul
  beta * Bs[i].T @ Ps[j] @ As[i] + Ns[i])
a=2.0: LQMarkov raised ValueError: array must not contain infs or NaNs
solver returned Ps = [nan]

Here P grows by beta * a**2 per iteration.

• a = 1.1 (1.21x): P is still finite at max_iter, so the intended message is raised.
• a = 2.0 (4x): error is inf at iteration 512 and nan at iteration 513, where the loop exits. NumPy's RuntimeWarnings are the only signal. The direct call prints none because Python shows each warning only once per location.

Bisecting on this example, growth above about 2.03x per iteration overflows before the default max_iter. The solver also returns NaN in these cases:

• the same model with beta = 0.95 (at iteration 532);
• the same model with two Markov states;
• larger problems (n = 2 or 3, k = 1 or 2) that have an uncontrollable unstable mode.

Suggested fix

Raise the existing fail_msg as soon as error is not finite:

--- a/quantecon/_matrix_eqn.py
+++ b/quantecon/_matrix_eqn.py
@@ -313,6 +313,9 @@
                 Ps1[i][:, :] = Rs[i] + sum1 - sum2
                 error += np.max(np.abs(Ps1[i] - Ps[i]))
 
+            if not np.isfinite(error):
+                raise ValueError(fail_msg.format(iteration + 1))
+
             Ps[:, :, :] = Ps1[:, :, :]
             iteration += 1

Reusing fail_msg keeps the beta >= 1 explanation. The count is the iteration at which the values stopped being finite. With the patch, the reproducer's a = 2.0 case raises this from both LQMarkov and the direct call: Convergence failed after 512 iterations. When beta>=1, a stationary solution may not exist if the system cannot be stabilized without discounting. Check the model specification. Only the line-307 overflow warning is still printed.

The check adds no arithmetic, so converging runs are unchanged:

• On a patched scratch copy of main, quantecon/tests/test_matrix_eqn.py and quantecon/tests/test_lqcontrol.py still pass (19 passed, same as unpatched).
• In a seeded batch of 400 random problems, every converging case returns bitwise-identical Ps. The only change is that cases which returned non-finite Ps now raise.

A possible test for test_matrix_eqn.py. It fails on main with AssertionError: ValueError not raised and passes with the patch:

def test_solve_discrete_riccati_system_overflow_raises():
    # x' = 2x, B = 0: P grows 4x per iteration and overflows before max_iter
    one, zero = np.ones((1, 1, 1)), np.zeros((1, 1, 1))
    with assert_raises(ValueError) as excinfo:
        qme.solve_discrete_riccati_system(np.ones((1, 1)), 2 * one, zero,
                                          zero, one, one, zero, beta=1.0)
    assert "stabilized without discounting" in str(excinfo.exception)

The doubling solver solve_discrete_riccati, used by LQ and Kalman, exits the same way (L209-L225):

• On main, solve_discrete_riccati(2, 0, 1, 1) returns [[nan]], and LQ(1, 1, 2, 0).stationary_values() fails with the same SciPy message.
• With the same guard after its error = ... line (raise ValueError(fail_msg.format(i))), both raise Convergence failed after 10 iterations.
• With both guards, test_matrix_eqn, test_lqcontrol, test_ricatti, test_kalman, test_robustlq and test_lqnash still pass (45 passed).

Context

Found while verifying #984, the fix for #982. #984 rewrites the update inside this loop but keeps the while error > tolerance exit and adds no finiteness check.

🤖 Generated with Claude Code