Cholesky Factorization Example — cholfact and cholsolve¶
Demonstrates Cholesky factorization on a 5×5 symmetric positive definite
matrix (scaled Hilbert matrix) using linearsystems.cholfact() and
linearsystems.cholsolve(), with accuracy compared against the direct
SPD solver linearsystems.lslsf().
Example Code¶
"""IMSL Cholesky Factorization example: factor and solve with cholfact/cholsolve.
Creates a 5×5 symmetric positive definite matrix (scaled Hilbert matrix),
factors it with cholfact(), solves with cholsolve(), and compares accuracy
against the direct SPD solver lslsf().
Outputs:
- Factorization timing and residual norms printed to stdout
- SVG heatmap of A, L factor, and residual comparison saved to
test_output/example_imsl_cholesky_solve.svg
"""
from __future__ import annotations
import time
from pathlib import Path
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
from linearsystems import cholfact, cholsolve, lslsf
def _hilbert(n: int) -> np.ndarray:
"""Return the n×n Hilbert matrix.
H[i,j] = 1 / (i + j + 1).
Args:
n (int): Matrix dimension.
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.0)
def run_demo_imsl_cholesky_solve() -> Dict[str, object]:
"""Run Cholesky factorization demo on a 5×5 scaled Hilbert matrix.
Factors with cholfact(), solves with cholsolve(), and compares against
lslsf() (direct SPD solver). Both timing and residuals are reported.
Args:
None
Returns:
Dict[str, object]: Result dict with keys ``x_chol``, ``x_direct``,
``residual_chol``, ``residual_direct``, ``time_chol``,
``time_direct``, and ``plot_path`` (str).
"""
n = 5
# Scale to improve conditioning while keeping SPD
a = _hilbert(n)
a = a / np.max(a) * 10.0 + np.eye(n) * n # ensure strict positive definiteness
rng = np.random.default_rng(7)
x_true = rng.standard_normal(n)
b = a @ x_true
# Cholesky: factor then solve
t0 = time.perf_counter()
c_low = cholfact(a)
x_chol = cholsolve(c_low, b)
t_chol = time.perf_counter() - t0
res_chol = float(np.linalg.norm(a @ x_chol - b))
# Direct SPD solve (lslsf)
t0 = time.perf_counter()
result_direct = lslsf(a, b)
t_direct = time.perf_counter() - t0
x_direct = result_direct.x
res_direct = result_direct.residual
# Extract lower-triangular Cholesky factor L
# c_low is (packed_c, lower) from scipy.linalg.cho_factor
packed_c = c_low[0]
L = np.tril(packed_c)
print("\nIMSL Cholesky Solve: 5×5 scaled Hilbert matrix")
print("-" * 60)
print(f"{'Method':<25} {'‖Ax-b‖':>12} {'‖x-x_true‖':>14} {'Time (s)':>10}")
print("-" * 60)
err_chol = float(np.linalg.norm(x_chol - x_true))
err_direct = float(np.linalg.norm(x_direct - x_true))
print(f"{'cholfact+cholsolve':<25} {res_chol:>12.4e} {err_chol:>14.4e} {t_chol:>10.6f}")
print(f"{'lslsf (direct SPD)':<25} {res_direct:>12.4e} {err_direct:>14.4e} {t_direct:>10.6f}")
print("-" * 60)
print("\nCholesky factor L (lower triangle):")
for row in L:
print(" " + " ".join(f"{v:8.4f}" for v in row))
print("-" * 60)
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_cholesky_solve.svg"
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
def _heatmap(ax: plt.Axes, mat: np.ndarray, title: str) -> None:
"""Draw an annotated heatmap on *ax*.
Args:
ax (plt.Axes): Target axes.
mat (np.ndarray): 2D matrix to visualize.
title (str): Axes title.
"""
im = ax.imshow(mat, cmap="viridis", aspect="auto")
fig.colorbar(im, ax=ax, shrink=0.8)
ax.set_title(title)
for (r, c), val in np.ndenumerate(mat):
ax.text(c, r, f"{val:.2f}", ha="center", va="center", fontsize=7, color="white")
_heatmap(axes[0], a, "Matrix A (Scaled Hilbert 5×5)")
_heatmap(axes[1], L, "Cholesky Factor L (lower)")
# Right: component-wise error comparison
idx = np.arange(n)
axes[2].plot(idx, x_true, "ko-", label="x_true", linewidth=1.5)
axes[2].plot(idx, x_chol, "bs--", label="cholfact+cholsolve", linewidth=1.5)
axes[2].plot(idx, x_direct, "r^:", label="lslsf direct", linewidth=1.5)
axes[2].set_xticks(idx)
axes[2].set_xticklabels([f"x[{i}]" for i in idx])
axes[2].set_title("Solution Components Comparison")
axes[2].set_ylabel("Value")
axes[2].legend(fontsize=8)
axes[2].grid(True, alpha=0.3)
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
print(f"\nPlot saved to: {plot_path}")
return {
"x_chol": x_chol,
"x_direct": x_direct,
"residual_chol": res_chol,
"residual_direct": res_direct,
"time_chol": t_chol,
"time_direct": t_direct,
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_imsl_cholesky_solve()
Plot Output¶
Left: heatmap of the 5×5 scaled Hilbert matrix A. Centre: lower Cholesky factor L. Right: per-component solution comparison between cholfact/cholsolve and lslsf.¶
Console Output¶
IMSL Cholesky Solve: 5×5 scaled Hilbert matrix
------------------------------------------------------------
Method ΓÇûAx-bΓÇû ΓÇûx-x_trueΓÇû Time (s)
------------------------------------------------------------
cholfact+cholsolve 2.0138e-15 2.3240e-16 0.000220
lslsf (direct SPD) 2.0138e-15 2.3240e-16 0.000198
------------------------------------------------------------
Cholesky factor L (lower triangle):
3.8730 0.0000 0.0000 0.0000 0.0000
5.0000 2.5820 0.0000 0.0000 0.0000
3.3333 2.5000 2.4433 0.0000 0.0000
2.5000 2.0000 1.6667 2.3836 0.0000
2.0000 1.6667 1.4286 1.2500 2.3503
------------------------------------------------------------
Plot saved to: test_output\example_imsl_cholesky_solve.svg