Matrix Norms and Conditioning¶
Demonstrates vector norms, matrix norms, condition number growth for Hilbert matrices, reciprocal condition number estimation, and the Moore-Penrose pseudo-inverse.
Functions covered: linearsystems.vector_norm(),
linearsystems.matrix_norm(),
linearsystems.condition_number(),
linearsystems.rcond_estimate(),
linearsystems.pseudo_inverse().
Example Code¶
"""Matrix Norms and Conditioning example.
Demonstrates:
- vector_norm: compute 1-, 2-, and inf-norms of several vectors
- matrix_norm: compute 1-, 2- (spectral), inf-, and Frobenius norms
- condition_number: show ill-conditioning of Hilbert matrices (sizes 2..8)
- rcond_estimate: compare with 1/condition_number
- pseudo_inverse: Moore-Penrose pseudo-inverse of a rank-deficient matrix,
verify all four Moore-Penrose conditions
Outputs:
- All norms and condition numbers printed to stdout
- SVG plot: condition number vs matrix size (log scale) + norm bar charts
saved to test_output/example_imsl_norms_conditioning.svg
"""
from __future__ import annotations
from pathlib import Path
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
from linearsystems import (
vector_norm,
matrix_norm,
condition_number,
rcond_estimate,
pseudo_inverse,
)
def _hilbert(n: int) -> np.ndarray:
"""Return the n×n Hilbert matrix H[i,j] = 1/(i+j+1).
Args:
n (int): Size of the Hilbert matrix.
Returns:
np.ndarray: Hilbert matrix of shape (n, n).
"""
i, j = np.meshgrid(np.arange(n), np.arange(n), indexing="ij")
return 1.0 / (i + j + 1)
def run_demo_imsl_norms_conditioning() -> Dict[str, object]:
"""Run matrix norms and conditioning demo.
Exercises vector_norm, matrix_norm, condition_number, rcond_estimate,
and pseudo_inverse. Reports all computed values and verifies the four
Moore-Penrose conditions for the pseudo-inverse.
Args:
None
Returns:
Dict[str, object]: Keys: ``cond_numbers`` (list), ``rcond_estimates``
(list), ``mp_errors`` (list of 4 floats), ``plot_path`` (str).
"""
rng = np.random.default_rng(13)
# ------------------------------------------------------------------
# 1. Vector norms
# ------------------------------------------------------------------
vectors = {
"v1 = [1,2,3,4,5]": np.array([1.0, 2, 3, 4, 5]),
"v2 = ones(8)": np.ones(8),
"v3 = randn(6)": rng.standard_normal(6),
}
print("\nMatrix Norms and Conditioning")
print("=" * 70)
print("\n1. Vector norms")
print(f" {'Vector':<22} {'1-norm':>10} {'2-norm':>10} {'inf-norm':>10}")
print(" " + "-" * 55)
for name, v in vectors.items():
n1 = vector_norm(v, ord=1)
n2 = vector_norm(v, ord=2)
ni = vector_norm(v, ord=np.inf)
print(f" {name:<22} {n1:>10.4f} {n2:>10.4f} {ni:>10.4f}")
# ------------------------------------------------------------------
# 2. Matrix norms
# ------------------------------------------------------------------
matrices = {
"A (3×3 random)": rng.standard_normal((3, 3)),
"I (5×5)": np.eye(5),
"Hilbert(4)": _hilbert(4),
}
print("\n2. Matrix norms")
print(
f" {'Matrix':<18} {'1-norm':>10} {'spectral':>10}"
f" {'inf-norm':>10} {'Frobenius':>10}"
)
print(" " + "-" * 62)
for name, M in matrices.items():
mn1 = matrix_norm(M, ord=1)
mn2 = matrix_norm(M, ord=2)
mni = matrix_norm(M, ord=np.inf)
mnf = matrix_norm(M, ord="fro")
print(f" {name:<18} {mn1:>10.4f} {mn2:>10.4f} {mni:>10.4f} {mnf:>10.4f}")
# ------------------------------------------------------------------
# 3. Condition numbers of Hilbert matrices
# ------------------------------------------------------------------
sizes = list(range(2, 9))
cond_numbers = []
rcond_estimates = []
print("\n3. Condition numbers of Hilbert matrices")
print(f" {'n':>4} {'cond(H)':>18} {'rcond_est':>14} {'1/cond':>14} {'match?':>8}")
print(" " + "-" * 65)
for n in sizes:
H = _hilbert(n)
c = condition_number(H)
rc = rcond_estimate(H)
ref_rc = 1.0 / c if c > 0 else 0.0
match = "yes" if abs(rc - ref_rc) / (ref_rc + 1e-300) < 0.01 else "~"
print(f" {n:>4} {c:>18.4e} {rc:>14.4e} {ref_rc:>14.4e} {match:>8}")
cond_numbers.append(c)
rcond_estimates.append(rc)
# ------------------------------------------------------------------
# 4. Pseudo-inverse of a rank-deficient matrix
# ------------------------------------------------------------------
# Build a 4×3 rank-2 matrix: rank-deficient
U2 = rng.standard_normal((4, 2))
V2 = rng.standard_normal((2, 3))
A_rd = U2 @ V2 # rank ≤ 2
A_pinv = pseudo_inverse(A_rd)
# Moore-Penrose conditions:
# 1. A @ pinv(A) @ A = A
mp1 = float(np.linalg.norm(A_rd @ A_pinv @ A_rd - A_rd))
# 2. pinv(A) @ A @ pinv(A) = pinv(A)
mp2 = float(np.linalg.norm(A_pinv @ A_rd @ A_pinv - A_pinv))
# 3. (A @ pinv(A)).T = A @ pinv(A)
P = A_rd @ A_pinv
mp3 = float(np.linalg.norm(P - P.T))
# 4. (pinv(A) @ A).T = pinv(A) @ A
Q = A_pinv @ A_rd
mp4 = float(np.linalg.norm(Q - Q.T))
mp_errors = [mp1, mp2, mp3, mp4]
print("\n4. Pseudo-inverse of rank-deficient 4×3 matrix (rank ≤ 2)")
print(f" A shape: {A_rd.shape}, pinv(A) shape: {A_pinv.shape}")
print(f" MP cond 1: ‖A pinv(A) A − A‖ = {mp1:.4e}")
print(f" MP cond 2: ‖pinv(A) A pinv(A)−pinv(A)‖ = {mp2:.4e}")
print(f" MP cond 3: ‖(A pinv(A))ᵀ − A pinv(A)‖ = {mp3:.4e}")
print(f" MP cond 4: ‖(pinv(A) A)ᵀ − pinv(A) A‖ = {mp4:.4e}")
# ------------------------------------------------------------------
# Plots
# ------------------------------------------------------------------
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_norms_conditioning.svg"
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
# (0) Condition number vs Hilbert matrix size
ax = axes[0]
ax.semilogy(sizes, cond_numbers, "b-o", linewidth=2, label="cond(H_n)")
ax.semilogy(sizes, [1.0 / r if r > 0 else np.nan for r in rcond_estimates],
"r--s", linewidth=1.5, label="1/rcond_est")
ax.set_xlabel("Matrix size n")
ax.set_ylabel("Condition number (log scale)")
ax.set_title("Condition number of Hilbert matrix H_n")
ax.legend()
ax.grid(True, alpha=0.3)
# (1) Vector norms bar chart
ax2 = axes[1]
vnames = list(vectors.keys())
x = np.arange(len(vnames))
n1s = [vector_norm(v, 1) for v in vectors.values()]
n2s = [vector_norm(v, 2) for v in vectors.values()]
nis = [vector_norm(v, np.inf) for v in vectors.values()]
w = 0.25
ax2.bar(x - w, n1s, w, label="1-norm", color="#1f77b4", alpha=0.85)
ax2.bar(x, n2s, w, label="2-norm", color="#ff7f0e", alpha=0.85)
ax2.bar(x + w, nis, w, label="inf-norm", color="#2ca02c", alpha=0.85)
ax2.set_xticks(x)
ax2.set_xticklabels([f"v{i+1}" for i in range(len(vnames))], fontsize=10)
ax2.set_title("Vector norms")
ax2.set_ylabel("Norm value")
ax2.legend()
ax2.grid(True, axis="y", alpha=0.3)
# (2) Moore-Penrose condition errors
ax3 = axes[2]
mp_labels = ["MP1\nA pinv A = A", "MP2\npinv A pinv=pinv",
"MP3\n(A pinv)ᵀ=A pinv", "MP4\n(pinv A)ᵀ=pinv A"]
bar_colors = ["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728"]
bars = ax3.bar(range(4), [max(e, 1e-18) for e in mp_errors],
color=bar_colors, alpha=0.85)
ax3.set_xticks(range(4))
ax3.set_xticklabels(mp_labels, fontsize=8)
ax3.set_yscale("log")
ax3.set_title("Moore-Penrose conditions (pseudo-inverse)")
ax3.set_ylabel("‖error‖₂ (log scale)")
ax3.grid(True, axis="y", alpha=0.3)
for bar, val in zip(bars, mp_errors):
ax3.text(bar.get_x() + bar.get_width() / 2, bar.get_height() * 1.5,
f"{val:.1e}", ha="center", va="bottom", fontsize=7)
fig.suptitle("Matrix Norms and Conditioning", fontsize=14)
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
print(f"\nPlot saved to: {plot_path}")
return {
"cond_numbers": cond_numbers,
"rcond_estimates": rcond_estimates,
"mp_errors": mp_errors,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_norms_conditioning()
Plot Output¶
Left: condition number of Hilbert matrices vs size (log scale). Centre: vector 1-, 2-, and inf-norms for three vectors. Right: Moore-Penrose condition errors for the pseudo-inverse.¶
Console Output¶
Matrix Norms and Conditioning
======================================================================
1. Vector norms
Vector 1-norm 2-norm inf-norm
-------------------------------------------------------
v1 = [1,2,3,4,5] 15.0000 7.4162 5.0000
v2 = ones(8) 8.0000 2.8284 1.0000
v3 = randn(6) 7.6367 3.9525 3.0783
2. Matrix norms
Matrix 1-norm spectral inf-norm Frobenius
--------------------------------------------------------------
A (3×3 random) 2.6550 2.0635 2.3752 2.3092
I (5×5) 1.0000 1.0000 1.0000 2.2361
Hilbert(4) 2.0833 1.5002 2.0833 1.5097
3. Condition numbers of Hilbert matrices
n cond(H) rcond_est 1/cond match?
-----------------------------------------------------------------
2 1.9281e+01 5.1863e-02 5.1863e-02 yes
3 5.2406e+02 1.9082e-03 1.9082e-03 yes
4 1.5514e+04 6.4459e-05 6.4459e-05 yes
5 4.7661e+05 2.0982e-06 2.0982e-06 yes
6 1.4951e+07 6.6885e-08 6.6885e-08 yes
7 4.7537e+08 2.1036e-09 2.1036e-09 yes
8 1.5258e+10 6.5541e-11 6.5541e-11 yes
4. Pseudo-inverse of rank-deficient 4×3 matrix (rank ≤ 2)
A shape: (4, 3), pinv(A) shape: (3, 4)
MP cond 1: ‖A pinv(A) A − A‖ = 2.3917e-15
MP cond 2: ‖pinv(A) A pinv(A)−pinv(A)‖ = 9.1632e-17
MP cond 3: ‖(A pinv(A))ᵀ − A pinv(A)‖ = 4.5419e-16
MP cond 4: ‖(pinv(A) A)ᵀ − pinv(A) A‖ = 2.1138e-16
Plot saved to: test_output\example_imsl_norms_conditioning.svg