Function approximate_spectral_radius 

pub fn approximate_spectral_radius(
    matrix: &dyn Matrix,
    num_iters: usize,
) -> Option<Precision>

Estimate the spectral radius of D⁻¹O (the Neumann iteration matrix) via power iteration. Tighter than the 1 - coherence upper bound that optimal_neumann_terms uses by default.

Why

For DD A = D - O with coherence c, ‖D⁻¹O‖_∞ ≤ 1 - c. That’s a safe upper bound on ρ(D⁻¹O) but often loose — the actual spectral radius can be much smaller, especially on matrices with non-uniform row weights. A tight ρ estimate lets optimal_neumann_terms_with_rho pick a much smaller max_terms:

  k ≥ log(b_inf / (min_diag · tolerance)) / log(1 / ρ)

Smaller ρ → larger log(1/ρ) → smaller k. The auto-tuned closure shrinks proportionally, and per-event SubLinear cost drops.

Algorithm

Standard power iteration on M = D⁻¹O = I - D⁻¹A:

1. Start with a uniform unit vector v_0.
2. For each iteration, compute v_{k+1} = M v_k and the Rayleigh-style ratio ‖M v_k‖_∞ / ‖v_k‖_∞.
3. Renormalise v_{k+1} to unit infinity norm.
4. Return the last ratio (the dominant-eigenvalue estimate).

Cost: O(num_iters · nnz(A)). Convergence is geometric in λ_2 / λ_1 — typically 10–20 iterations suffice on well-conditioned DD matrices. Run once at matrix-build time; amortised across all subsequent events.

Returns

• Some(rho) on success. 0.0 ≤ rho < 1.0 for strict-DD matrices.
• None if num_iters == 0, the matrix is empty, or any diagonal row has zero entry.

Examples

let rho = approximate_spectral_radius(a, /*num_iters=*/ 20).unwrap_or(0.99);
// Use the tight rho instead of the loose (1 - coherence) bound.
let k = optimal_neumann_terms_with_rho(rho, /*b_inf=*/ 10.0, /*min_diag=*/ 5.0, /*tol=*/ 1e-8);