Triangular Systems and QR Solve¶
Demonstrates forward and backward substitution, IMSL-style triangular solvers, QR factorisation and least-squares solve, matrix inverse, and matrix power.
Functions covered: linearsystems.lsltr(),
linearsystems.lsltc(),
linearsystems.solve_lower_triangular(),
linearsystems.solve_upper_triangular(),
linearsystems.qrfact(),
linearsystems.qrsolve(),
linearsystems.matrix_inverse(),
linearsystems.matrix_power().
Example Code¶
"""Triangular Systems and QR Solve example.
Demonstrates:
- Forward/back substitution via solve_lower_triangular, solve_upper_triangular
- IMSL-style triangular solvers lsltr (real) and lsltc (complex)
- QR factorization (qrfact) and least-squares solve (qrsolve) for a 5×3 system
- matrix_inverse: verify A @ inv(A) ≈ I
- matrix_power: verify A³ = A @ A @ A
Outputs:
- Residuals and verification errors printed to stdout
- SVG heatmaps of L, U factors and bar chart of solutions saved to
test_output/example_imsl_triangular_qr.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 (
lsltr,
lsltc,
solve_upper_triangular,
solve_lower_triangular,
qrfact,
qrsolve,
matrix_inverse,
matrix_power,
)
def run_demo_imsl_triangular_qr() -> Dict[str, object]:
"""Run triangular-systems and QR example.
Builds explicit lower and upper triangular systems, solves them via
forward/back substitution (solve_lower_triangular, solve_upper_triangular)
and via lsltr. Factorises a 5×3 overdetermined system with qrfact and
solves it with qrsolve. Checks matrix_inverse and matrix_power.
Args:
None
Returns:
Dict[str, object]: Keys: ``res_lower``, ``res_upper``, ``res_lsltr``,
``qr_residual``, ``inv_error``, ``power_error``, ``plot_path``.
"""
rng = np.random.default_rng(7)
# ------------------------------------------------------------------
# 1. Lower triangular solve L x = b
# ------------------------------------------------------------------
n = 5
L = np.tril(rng.uniform(1, 4, (n, n))) # random lower-tri, pos diag
np.fill_diagonal(L, np.abs(np.diag(L)) + 1.0)
x_true = rng.standard_normal(n)
b_lower = L @ x_true
x_lt = solve_lower_triangular(L, b_lower)
res_lower = float(np.linalg.norm(L @ x_lt - b_lower))
# lsltr (lower=True) should give identical result
result_lsltr_lo = lsltr(L, b_lower, lower=True)
res_lsltr_lo = result_lsltr_lo.residual
print("\nTriangular Systems and QR Solve")
print("=" * 60)
print(f"\n1. Lower triangular solve (n={n})")
print(f" solve_lower_triangular residual : {res_lower:.4e}")
print(f" lsltr(lower=True) residual : {res_lsltr_lo:.4e}")
# ------------------------------------------------------------------
# 2. Upper triangular solve U x = b
# ------------------------------------------------------------------
U = np.triu(rng.uniform(1, 4, (n, n)))
np.fill_diagonal(U, np.abs(np.diag(U)) + 1.0)
b_upper = U @ x_true
x_ut = solve_upper_triangular(U, b_upper)
res_upper = float(np.linalg.norm(U @ x_ut - b_upper))
result_lsltr_up = lsltr(U, b_upper, lower=False)
res_lsltr_up = result_lsltr_up.residual
print(f"\n2. Upper triangular solve (n={n})")
print(f" solve_upper_triangular residual : {res_upper:.4e}")
print(f" lsltr(lower=False) residual : {res_lsltr_up:.4e}")
res_lsltr = max(res_lsltr_lo, res_lsltr_up)
# ------------------------------------------------------------------
# 3. Complex triangular solve via lsltc
# ------------------------------------------------------------------
Lc = np.tril(rng.uniform(1, 3, (n, n)) + 1j * rng.uniform(0, 0.5, (n, n)))
np.fill_diagonal(Lc, np.abs(np.diag(Lc)) + 1.0)
xc_true = rng.standard_normal(n) + 1j * rng.standard_normal(n)
bc = Lc @ xc_true
result_ltc = lsltc(Lc, bc, lower=True)
res_ltc = result_ltc.residual
print(f"\n3. Complex triangular solve lsltc (n={n})")
print(f" lsltc residual : {res_ltc:.4e}")
# ------------------------------------------------------------------
# 4. QR factorization and least-squares solve
# ------------------------------------------------------------------
m, k = 5, 3
A = rng.standard_normal((m, k))
x_ls_true = rng.standard_normal(k)
b_ls = A @ x_ls_true + 0.01 * rng.standard_normal(m) # noisy RHS
qr = qrfact(A)
x_qr = qrsolve(qr, b_ls)
qr_residual = float(np.linalg.norm(A @ x_qr - b_ls))
print(f"\n4. QR least-squares ({m}×{k} system)")
print(f" {'Component':<12} {'x_solved':>12}")
for i, v in enumerate(x_qr):
print(f" x[{i}] {v:>12.6f}")
print(f" ‖Ax - b‖₂ = {qr_residual:.4e}")
# ------------------------------------------------------------------
# 5. Matrix inverse: verify A @ inv(A) ≈ I
# ------------------------------------------------------------------
A5 = rng.standard_normal((n, n))
A5 = A5 @ A5.T + n * np.eye(n) # well-conditioned SPD
A5_inv = matrix_inverse(A5)
inv_error = float(np.linalg.norm(A5 @ A5_inv - np.eye(n)))
print(f"\n5. matrix_inverse (n={n})")
print(f" ‖A @ inv(A) − I‖₂ = {inv_error:.4e}")
# ------------------------------------------------------------------
# 6. Matrix power: verify A³ = A @ A @ A
# ------------------------------------------------------------------
A3 = rng.standard_normal((3, 3))
A3 = A3 @ A3.T + 3.0 * np.eye(3)
Ap3 = matrix_power(A3, 3)
Ap3_ref = A3 @ A3 @ A3
power_error = float(np.linalg.norm(Ap3 - Ap3_ref))
print(f"\n6. matrix_power(A, 3) (3×3)")
print(f" ‖A³ − A@A@A‖₂ = {power_error:.4e}")
# ------------------------------------------------------------------
# Plots
# ------------------------------------------------------------------
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_triangular_qr.svg"
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
# (0,0) L heatmap
ax = axes[0, 0]
im = ax.imshow(L, cmap="Blues")
ax.set_title("Lower triangular L (5×5)")
ax.set_xlabel("Column")
ax.set_ylabel("Row")
fig.colorbar(im, ax=ax, shrink=0.8)
# (0,1) U heatmap
ax = axes[0, 1]
im2 = ax.imshow(U, cmap="Oranges")
ax.set_title("Upper triangular U (5×5)")
ax.set_xlabel("Column")
ax.set_ylabel("Row")
fig.colorbar(im2, ax=ax, shrink=0.8)
# (0,2) bar chart: lower-tri solutions
ax = axes[0, 2]
idx = np.arange(n)
ax.bar(idx - 0.2, x_true, width=0.35, label="x_true", color="#1f77b4", alpha=0.85)
ax.bar(idx + 0.2, x_lt, width=0.35, label="x_solved (lower)", color="#ff7f0e", alpha=0.85)
ax.set_xticks(idx)
ax.set_xticklabels([f"x[{i}]" for i in idx])
ax.set_title("Lower-tri solution components")
ax.set_ylabel("Value")
ax.legend(fontsize=8)
ax.grid(True, axis="y", alpha=0.3)
# (1,0) QR factor Q heatmap
Q, R = qr
ax = axes[1, 0]
im3 = ax.imshow(Q, cmap="RdBu", aspect="auto")
ax.set_title("QR factor Q (5×3)")
ax.set_xlabel("Column")
ax.set_ylabel("Row")
fig.colorbar(im3, ax=ax, shrink=0.8)
# (1,1) QR solve: bar chart
ax = axes[1, 1]
idx3 = np.arange(k)
ax.bar(idx3 - 0.2, x_ls_true, width=0.35, label="x_ls_true", color="#2ca02c", alpha=0.85)
ax.bar(idx3 + 0.2, x_qr, width=0.35, label="x_qr", color="#d62728", alpha=0.85)
ax.set_xticks(idx3)
ax.set_xticklabels([f"x[{i}]" for i in idx3])
ax.set_title("QR least-squares solutions (5×3)")
ax.set_ylabel("Value")
ax.legend(fontsize=8)
ax.grid(True, axis="y", alpha=0.3)
# (1,2) residuals summary bar chart
ax = axes[1, 2]
labels = [
"solve_lower\ntriangular",
"solve_upper\ntriangular",
"lsltr\n(lower)",
"lsltr\n(upper)",
"lsltc\n(complex)",
"qrsolve\n(LS)",
"matrix_inv\nerror",
"matrix_pow\nerror",
]
vals = [res_lower, res_upper, res_lsltr_lo, res_lsltr_up, res_ltc,
qr_residual, inv_error, power_error]
colors = ["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
"#8c564b", "#e377c2", "#7f7f7f"]
ax.bar(range(len(labels)), vals, color=colors, alpha=0.85)
ax.set_xticks(range(len(labels)))
ax.set_xticklabels(labels, fontsize=7)
ax.set_yscale("log")
ax.set_title("Residuals / errors summary")
ax.set_ylabel("‖error‖₂ (log scale)")
ax.grid(True, axis="y", alpha=0.3)
fig.suptitle("Triangular Systems and QR Solve", fontsize=14)
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
print(f"\nPlot saved to: {plot_path}")
return {
"res_lower": res_lower,
"res_upper": res_upper,
"res_lsltr": res_lsltr,
"qr_residual": qr_residual,
"inv_error": inv_error,
"power_error": power_error,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_triangular_qr()
Plot Output¶
Top row: heatmaps of L and U factors and solution bar chart. Bottom row: Q factor heatmap, QR least-squares solution, and residuals summary.¶
Console Output¶
Triangular Systems and QR Solve
============================================================
1. Lower triangular solve (n=5)
solve_lower_triangular residual : 0.0000e+00
lsltr(lower=True) residual : 0.0000e+00
2. Upper triangular solve (n=5)
solve_upper_triangular residual : 0.0000e+00
lsltr(lower=False) residual : 0.0000e+00
3. Complex triangular solve lsltc (n=5)
lsltc residual : 5.8287e-16
4. QR least-squares (5×3 system)
Component x_solved
x[0] 0.756560
x[1] -0.857320
x[2] 0.773006
ΓÇûAx - bΓÇûΓéé = 1.1568e-02
5. matrix_inverse (n=5)
‖A @ inv(A) − I‖₂ = 4.0645e-16
6. matrix_power(A, 3) (3×3)
‖A³ − A@A@A‖₂ = 0.0000e+00
Plot saved to: test_output\example_imsl_triangular_qr.svg