Iterative Solvers Example — CG, GMRES, and BiCG

Demonstrates solving a 50×50 sparse tridiagonal system with three iterative methods — linearsystems.lsitcg() (Conjugate Gradient), linearsystems.lsigmr() (GMRES), and linearsystems.lsibcg() (BiConjugate Gradient) — and compares their convergence histories.

Example Code

"""IMSL Iterative Solvers example: CG, GMRES, and BiCG on a sparse tridiagonal system.

Creates a diagonally dominant 50×50 tridiagonal system and solves it with three
iterative methods: Conjugate Gradient (lsitcg), GMRES (lsigmr), and BiConjugate
Gradient (lsibcg). Convergence histories are collected and plotted.

Outputs:
- Table printed to stdout
- SVG convergence plot saved to test_output/example_imsl_iterative_solvers.svg
"""

from __future__ import annotations

from pathlib import Path
from typing import Dict, List

import matplotlib.pyplot as plt
import numpy as np
import scipy.sparse.linalg as spla

from linearsystems import lsitcg, lsibcg, lsigmr


def _build_tridiagonal(n: int) -> np.ndarray:
    """Build a symmetric diagonally dominant n×n tridiagonal matrix.

    Main diagonal = 4, sub/super diagonals = -1.

    Args:
        n (int): System size.

    Returns:
        np.ndarray: Dense matrix of shape (n, n).
    """
    a = 4.0 * np.eye(n) - np.diag(np.ones(n - 1), 1) - np.diag(np.ones(n - 1), -1)
    return a


def _residuals_cg(
    a: np.ndarray,
    b: np.ndarray,
    tol: float = 1e-10,
) -> List[float]:
    """Run scipy CG and collect per-iteration residual norms.

    Args:
        a (np.ndarray): SPD matrix of shape (n, n).
        b (np.ndarray): Right-hand side vector of shape (n,).
        tol (float): Convergence tolerance. Defaults to ``1e-10``.

    Returns:
        List[float]: Residual norm at each iteration.
    """
    history: List[float] = [float(np.linalg.norm(b))]

    def cb(xk: np.ndarray) -> None:
        history.append(float(np.linalg.norm(a @ xk - b)))

    spla.cg(a, b, rtol=tol, atol=0.0, callback=cb)
    return history


def _residuals_gmres(
    a: np.ndarray,
    b: np.ndarray,
    tol: float = 1e-10,
) -> List[float]:
    """Run scipy GMRES and collect per-iteration residual norms.

    Args:
        a (np.ndarray): Square matrix of shape (n, n).
        b (np.ndarray): Right-hand side vector of shape (n,).
        tol (float): Convergence tolerance. Defaults to ``1e-10``.

    Returns:
        List[float]: Residual norm at each iteration.
    """
    history: List[float] = [float(np.linalg.norm(b))]

    def cb(pr_norm: float) -> None:
        history.append(float(pr_norm))

    spla.gmres(a, b, rtol=tol, atol=0.0, callback=cb, callback_type="pr_norm")
    return history


def _residuals_bicg(
    a: np.ndarray,
    b: np.ndarray,
    tol: float = 1e-10,
) -> List[float]:
    """Run scipy BiCG and collect per-iteration residual norms.

    Args:
        a (np.ndarray): Square matrix of shape (n, n).
        b (np.ndarray): Right-hand side vector of shape (n,).
        tol (float): Convergence tolerance. Defaults to ``1e-10``.

    Returns:
        List[float]: Residual norm at each iteration.
    """
    history: List[float] = [float(np.linalg.norm(b))]

    def cb(xk: np.ndarray) -> None:
        history.append(float(np.linalg.norm(a @ xk - b)))

    spla.bicg(a, b, rtol=tol, atol=0.0, callback=cb)
    return history


def run_demo_imsl_iterative_solvers() -> Dict[str, object]:
    """Run IMSL iterative solvers demo on a 50×50 tridiagonal system.

    Uses lsitcg (CG), lsigmr (GMRES), and lsibcg (BiCG) from the linearsystems
    package to solve the same system and compares convergence and accuracy.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``n_iter_cg``, ``n_iter_gmres``,
            ``n_iter_bicg``, ``residual_cg``, ``residual_gmres``,
            ``residual_bicg``, and ``plot_path`` (str).
    """
    n = 50
    a = _build_tridiagonal(n)
    rng = np.random.default_rng(0)
    x_true = rng.standard_normal(n)
    b = a @ x_true

    tol = 1e-10

    # Solve with all three iterative methods
    res_cg = lsitcg(a, b, tol=tol)
    res_gmres = lsigmr(a, b, tol=tol)
    res_bicg = lsibcg(a, b, tol=tol)

    # Collect convergence histories for plotting
    hist_cg = _residuals_cg(a, b, tol=tol)
    hist_gmres = _residuals_gmres(a, b, tol=tol)
    hist_bicg = _residuals_bicg(a, b, tol=tol)

    print("\nIMSL Iterative Solvers: 50×50 tridiagonal system (seed=0)")
    print("-" * 60)
    print(f"{'Method':<20} {'Iters':>8} {'Final ‖r‖':>14} {'‖x-x_true‖':>14}")
    print("-" * 60)
    methods = [
        ("CG  (lsitcg)", res_cg),
        ("GMRES (lsigmr)", res_gmres),
        ("BiCG (lsibcg)", res_bicg),
    ]
    for name, res in methods:
        err = float(np.linalg.norm(res.x - x_true))
        print(f"{name:<20} {res.n_iter:>8d} {res.residual:>14.4e} {err:>14.4e}")
    print("-" * 60)

    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_iterative_solvers.svg"

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Left: convergence histories
    ax = axes[0]
    ax.semilogy(hist_cg, "b-o", markersize=3, label=f"CG ({len(hist_cg)-1} iters)")
    ax.semilogy(hist_gmres, "r-s", markersize=3, label=f"GMRES ({len(hist_gmres)-1} iters)")
    ax.semilogy(hist_bicg, "g-^", markersize=3, label=f"BiCG ({len(hist_bicg)-1} iters)")
    ax.axhline(tol, color="gray", linestyle="--", linewidth=0.8, label=f"tol={tol:.0e}")
    ax.set_xlabel("Iteration")
    ax.set_ylabel("Residual norm ‖r‖")
    ax.set_title("Convergence: 50×50 Tridiagonal System")
    ax.legend()
    ax.grid(True, which="both", alpha=0.3)

    # Right: final residuals comparison
    ax2 = axes[1]
    names = ["CG\n(lsitcg)", "GMRES\n(lsigmr)", "BiCG\n(lsibcg)"]
    residuals = [res_cg.residual, res_gmres.residual, res_bicg.residual]
    colors = ["#1f77b4", "#d62728", "#2ca02c"]
    bars = ax2.bar(names, residuals, color=colors, alpha=0.85)
    ax2.set_yscale("log")
    ax2.set_title("Final Residual ‖Ax−b‖ by Method")
    ax2.set_ylabel("Residual norm (log scale)")
    for bar, val in zip(bars, residuals):
        ax2.text(
            bar.get_x() + bar.get_width() / 2,
            val * 2,
            f"{val:.2e}",
            ha="center",
            fontsize=9,
        )
    ax2.grid(True, axis="y", alpha=0.3)

    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)
    print(f"\nPlot saved to: {plot_path}")

    return {
        "n_iter_cg": res_cg.n_iter,
        "n_iter_gmres": res_gmres.n_iter,
        "n_iter_bicg": res_bicg.n_iter,
        "residual_cg": res_cg.residual,
        "residual_gmres": res_gmres.residual,
        "residual_bicg": res_bicg.residual,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_imsl_iterative_solvers()

Plot Output

Left: residual norm vs iteration for all three methods on the 50×50 tridiagonal system. Right: final residual ‖Ax−b‖ comparison (log scale).

Console Output

E:\ITDS\itds\packages\linearsystems\examples\example_imsl_iterative_solvers.py:194: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  fig.tight_layout()

IMSL Iterative Solvers: 50×50 tridiagonal system (seed=0)
------------------------------------------------------------
Method                  Iters      Final ΓÇûrΓÇû     ΓÇûx-x_trueΓÇû
------------------------------------------------------------
CG  (lsitcg)               18     1.0891e-09     3.2894e-10
GMRES (lsigmr)             18     1.0396e-09     3.7173e-10
BiCG (lsibcg)              18     1.0891e-09     3.2894e-10
------------------------------------------------------------

Plot saved to: test_output\example_imsl_iterative_solvers.svg