Condition Number Example — Matrix Conditioning and Stability¶
This example builds a sequence of matrices whose spectral condition numbers range from 1 (perfectly conditioned) to 10¹² (nearly singular) and studies the effect on linear system solve accuracy.
Key functions used:
eigensystems.condition_number()— spectral condition number κ(A)eigensystems.matrix_norm()— 2-norm ‖A‖₂
A right-hand side perturbed by Gaussian noise of magnitude 10⁻⁸ is used to expose how the relative solution error grows proportionally to κ(A).
Example Code¶
"""IMSL condition number example: matrix conditioning and linear system stability.
Demonstrates:
- Create matrices with a range of condition numbers (well- to ill-conditioned).
- Compute condition numbers via :func:`eigensystems.condition_number`.
- Compute matrix norms via :func:`eigensystems.matrix_norm`.
- Solve perturbed linear systems and measure solution error.
- Plot solution error vs condition number on a log-log scale.
- Save to test_output/example_imsl_condition_number.svg
Outputs:
- Printed table of condition numbers and solution errors.
- SVG figure saved to test_output/.
- Console output captured to docs/_generated_text/condition_number_output.txt
"""
from __future__ import annotations
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from eigensystems import condition_number, matrix_norm
def _build_matrix(n: int, cond_target: float) -> np.ndarray:
"""Build an n×n matrix with approximate condition number cond_target.
Uses a diagonal matrix with singular values spread logarithmically
from 1 to cond_target, then applies random orthogonal transforms.
Args:
n (int): Matrix dimension.
cond_target (float): Desired spectral condition number (σ_max/σ_min).
Returns:
np.ndarray: Real matrix of shape (n, n) with condition ≈ cond_target.
"""
rng = np.random.default_rng(seed=int(cond_target))
singular_vals = np.logspace(0, np.log10(cond_target), n)
Q1, _ = np.linalg.qr(rng.standard_normal((n, n)))
Q2, _ = np.linalg.qr(rng.standard_normal((n, n)))
return Q1 @ np.diag(singular_vals) @ Q2.T
def run_demo_condition_number() -> dict:
"""Run IMSL condition number example: effect on linear system accuracy.
Builds matrices with condition numbers from 1 to 1e12, solves Ax=b with
a small right-hand-side perturbation, and measures the resulting solution
error.
Args:
None
Returns:
dict: Result dictionary with keys ``cond_numbers`` (list[float]),
``errors`` (list[float]), and ``plot_path`` (str).
"""
n = 6
target_conds = [1.0, 10.0, 1e2, 1e4, 1e6, 1e8, 1e10, 1e12]
cond_numbers: list[float] = []
norms: list[float] = []
errors: list[float] = []
rng = np.random.default_rng(0)
x_true = np.ones(n)
perturbation_scale = 1e-8
print("\nIMSL Condition Number Example: Linear System Stability")
print("-" * 72)
print(f"{'Target κ':>12} {'Actual κ':>14} {'||A||₂':>12} {'||Δx||/||x||':>14}")
print("-" * 72)
for target in target_conds:
A = _build_matrix(n, target)
b_exact = A @ x_true
# Perturb b to simulate noise
b_perturbed = b_exact + perturbation_scale * rng.standard_normal(n)
kappa = condition_number(A)
norm_a = matrix_norm(A, ord=2)
try:
x_computed = np.linalg.solve(A, b_perturbed)
rel_err = float(np.linalg.norm(x_computed - x_true) / np.linalg.norm(x_true))
except np.linalg.LinAlgError:
rel_err = float("nan")
cond_numbers.append(kappa)
norms.append(norm_a)
errors.append(rel_err)
print(f" {target:>12.1e} {kappa:>14.4e} {norm_a:>12.4e} {rel_err:>14.4e}")
print("-" * 72)
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_condition_number.svg"
finite_mask = [not np.isnan(e) for e in errors]
kappa_plot = [k for k, m in zip(cond_numbers, finite_mask) if m]
err_plot = [e for e, m in zip(errors, finite_mask) if m]
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
fig.suptitle(
"IMSL condition_number / matrix_norm: Effect of Condition Number on Solve Accuracy",
fontsize=11,
)
ax = axes[0]
ax.loglog(kappa_plot, err_plot, "o-", color="#be123c", linewidth=1.8,
markersize=7, markerfacecolor="#fecdd3", markeredgecolor="#be123c")
ax.set_xlabel("Condition number κ(A) [log scale]")
ax.set_ylabel("Relative solution error ||Δx||/||x|| [log scale]")
ax.set_title("Solution Error vs Condition Number")
ax.grid(True, which="both", linestyle="--", alpha=0.5)
ax = axes[1]
ax.semilogx(kappa_plot, norms, "s-", color="#0369a1", linewidth=1.8,
markersize=7, markerfacecolor="#bae6fd", markeredgecolor="#0369a1")
ax.set_xlabel("Condition number κ(A) [log scale]")
ax.set_ylabel("2-norm ||A||₂")
ax.set_title("Matrix 2-Norm vs Condition Number")
ax.grid(True, which="both", linestyle="--", alpha=0.5)
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
print(f"\nPlot saved: {plot_path}")
return {
"cond_numbers": cond_numbers,
"errors": errors,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_condition_number()
Plot Output¶
Left: relative solution error ‖Δx‖/‖x‖ vs condition number on a log-log scale — the near-linear trend confirms the classical bound Δx/x ≤ κ·Δb/b. Right: 2-norm ‖A‖₂ vs condition number.¶
Console Output¶
IMSL Condition Number Example: Linear System Stability
------------------------------------------------------------------------
Target κ Actual κ ||A||₂ ||Δx||/||x||
------------------------------------------------------------------------
1.0e+00 1.0000e+00 1.0000e+00 3.8125e-09
1.0e+01 1.0000e+01 1.0000e+01 2.6560e-09
1.0e+02 1.0000e+02 1.0000e+02 2.5839e-09
1.0e+04 1.0000e+04 1.0000e+04 5.6974e-09
1.0e+06 1.0000e+06 1.0000e+06 1.2996e-09
1.0e+08 1.0000e+08 1.0000e+08 1.1888e-09
1.0e+10 1.0000e+10 1.0000e+10 2.0054e-07
1.0e+12 1.0000e+12 1.0000e+12 1.7846e-05
------------------------------------------------------------------------
Plot saved: test_output\example_imsl_condition_number.svg