DiscreteDP: opt-in iterative (Krylov) solver for sparse evaluate_policy

Summary

For the sparse formulation, DiscreteDP.evaluate_policy solves (I - beta Q_sigma) v = R_sigma with scipy.sparse.linalg.spsolve (SuperLU). On sparsity patterns without exploitable structure, the LU factorization suffers from severe fill-in, making policy evaluation — and hence policy iteration — the bottleneck of the sparse formulation. Krylov methods are immune to fill-in and exploit the structure of this particular system extremely well; this issue proposes adding an opt-in iterative solver option, with the direct solve remaining the default.

Measurements

On a random-pattern instance with n = 3000 states and 5 nonzeros per row of Q_sigma (beta = 0.95, rtol=1e-12):

SuperLU tuning yields at most ~1.4x; the iterative solvers yield ~200-650x on this pattern.

Why Krylov methods work so well here

A = I - beta Q_sigma is strictly (row-)diagonally dominant — each off-diagonal row sums to at most beta < 1 in absolute value — with condition number bounded by (1 + beta) / (1 - beta) independently of the sparsity pattern. Krylov iterations therefore converge in a small, essentially pattern-independent number of steps, each costing one O(nnz) matvec, while the direct factorization's cost is governed by fill-in, which for scattered patterns approaches dense complexity.

The practical consequence shows up in method rankings: benchmarking the three Q representations across methods shows that for sparse Q with a random pattern, policy iteration — the fastest method in the dense world — is often the slowest method (behind even value iteration at low fill), purely because of the LU fill-in in evaluate_policy. An iterative option would restore policy iteration's expected ranking on such models. (Where the pattern is fill-friendly, e.g. k = 1, the direct solve is already fast; the proposal changes nothing there by default.)

Why opt-in rather than a default change

An iterative solve makes policy evaluation tolerance-based rather than exact-up-to-roundoff, which touches policy iteration's finite-termination/exactness property: near-ties in compute_greedy could in principle resolve differently. The measured deviations are ~1e-12 at rtol=1e-12, but the semantic change should be the user's explicit choice.

Proposed API (for discussion)

A solver option on the constructor or solve, e.g. DiscreteDP(..., lineq_solver='direct') with 'direct' (current behavior, default) and 'bicgstab' (or a callable for full flexibility), used by evaluate_policy in the sparse case. Alternatives welcome — the main design questions are where the option lives, whether to expose the tolerance, and whether to fall back to the direct solve on non-convergence.

For cross-implementation reference: QuantEcon.jl is deliberately not adding this for now — its sparse evaluate_policy dispatches to UMFPACK, which handles these patterns ~19x better than SuperLU, and an iterative option there would require adding a Krylov package dependency. If this proposal proves its worth here, the Julia side can revisit.

Context: this came out of the performance work in #855/#857 and QuantEcon/QuantEcon.jl#388; the measurement scripts are reproducible from the description in #857.

🤖 Generated with Claude Code (Claude Fable 5)

Start with the sparse path in DiscreteDP.evaluate_policy and review the reproducible measurement scripts described in #857. Resolve where the opt-in solver and tolerance belong, preserve the direct solver as the default, and compare iterative results and convergence behavior against spsolve on the described sparse benchmarks.