Matrix Utility Functions — IS_SYMMETRIC, IS_POSITIVE_DEFINITE, SYMMETRIZE, ORTHOGONALIZE, PSEUDO_INVERSE¶
This example exercises the five IMSL-compatible matrix utility functions:
eigensystems.is_symmetric()— checks \(\max|A - A^H| < \varepsilon\)eigensystems.is_positive_definite()— attempts Cholesky decompositioneigensystems.symmetrize()— returns the Hermitian part \((A + A^H)/2\); verified witheigensystems.is_symmetric()eigensystems.orthogonalize()— economy QR delivers orthonormal columns; verified via \(\|Q^T Q - I\|_\infty\)eigensystems.pseudo_inverse()— Moore–Penrose pseudo-inverse of a rank-deficient matrix; verified via \(\|A \cdot A^+ \cdot A - A\|_\infty\)
Example Code¶
"""IMSL matrix utility examples: is_symmetric, is_positive_definite,
symmetrize, orthogonalize, pseudo_inverse.
Demonstrates:
- is_symmetric and is_positive_definite on various matrices
- symmetrize: asymmetric → symmetric, verify result
- orthogonalize: 4 random vectors → orthonormal, verify Q^T Q = I
- pseudo_inverse: rank-deficient matrix, verify A @ pinv(A) @ A = A
Outputs:
- Boolean check results and reconstruction errors
- SVG heatmaps: before/after for symmetrize and orthogonalize
- Save SVG to test_output/example_imsl_matrix_utils.svg
- Console output captured to docs/_generated_text/matrix_utils_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 (
is_positive_definite,
is_symmetric,
orthogonalize,
pseudo_inverse,
symmetrize,
)
def run_demo_matrix_utils() -> dict:
"""Run all matrix utility demonstrations.
Args:
None
Returns:
dict: Keys ``sym_result`` (np.ndarray), ``Q`` (np.ndarray),
``pinv_A`` (np.ndarray), ``orth_err`` (float),
``pinv_err`` (float), ``plot_path`` (str).
"""
# ------------------------------------------------------------------
# is_symmetric / is_positive_definite — multiple test matrices
# ------------------------------------------------------------------
S_sym = np.array([[4.0, 2.0, 1.0],
[2.0, 5.0, 3.0],
[1.0, 3.0, 6.0]]) # symmetric SPD
S_asym = np.array([[1.0, 2.0, 3.0],
[0.0, 4.0, 5.0],
[0.0, 0.0, 6.0]]) # upper triangular, not symmetric
S_indef = np.array([[ 1.0, 2.0],
[ 2.0, -1.0]]) # indefinite
S_herm = np.array([[3.0 + 0j, 1.0 + 2j],
[1.0 - 2j, 5.0 + 0j]]) # complex Hermitian
print("\nIMSL IS_SYMMETRIC / IS_POSITIVE_DEFINITE Checks")
print("-" * 55)
matrices = {
"S_sym (3×3 SPD)": S_sym,
"S_asym (upper tri)": S_asym,
"S_indef (indefinite)": S_indef,
"S_herm (Hermitian)": S_herm,
}
for name, mat in matrices.items():
sym = is_symmetric(mat)
spd = is_positive_definite(mat)
print(f" {name:<28} is_symmetric={sym!s:<5} is_positive_definite={spd!s}")
# ------------------------------------------------------------------
# symmetrize — make an asymmetric matrix symmetric
# ------------------------------------------------------------------
A_raw = np.array([[1.0, 4.0, 7.0],
[2.0, 5.0, 8.0],
[3.0, 6.0, 9.0]])
A_sym = symmetrize(A_raw)
sym_error = np.max(np.abs(A_sym - A_sym.T))
print("\nIMSL SYMMETRIZE Example")
print("Original (asymmetric):")
print(A_raw)
print("Symmetrized (A + A^T)/2:")
print(A_sym)
print(f"Max asymmetry after symmetrize: {sym_error:.2e} (should be ~0)")
print(f"is_symmetric(symmetrize(A_raw)) = {is_symmetric(A_sym)}")
# ------------------------------------------------------------------
# orthogonalize — 4 random 6-D vectors → orthonormal Q
# ------------------------------------------------------------------
rng = np.random.default_rng(7)
V_raw = rng.standard_normal((6, 4)) # 6 rows, 4 columns
Q = orthogonalize(V_raw)
orth_err = np.max(np.abs(Q.T @ Q - np.eye(4)))
print("\nIMSL ORTHOGONALIZE Example")
print(f"Input V shape: {V_raw.shape} → Q shape: {Q.shape}")
print(f"Q^T Q ≈ I, max error = {orth_err:.2e} (should be ~0)")
# ------------------------------------------------------------------
# pseudo_inverse — rank-deficient 3×4 matrix
# ------------------------------------------------------------------
# Rank-2 matrix: columns 3 and 4 are linear combinations of 1 and 2
C = np.array([
[1.0, 0.0, 1.0, 2.0],
[0.0, 1.0, 1.0, -1.0],
[1.0, 1.0, 2.0, 1.0],
])
C_pinv = pseudo_inverse(C)
recon_err = np.max(np.abs(C @ C_pinv @ C - C))
print("\nIMSL PSEUDO_INVERSE Example")
print(f"C shape: {C.shape}, pseudo_inverse shape: {C_pinv.shape}")
print(f"Reconstruction error ‖C @ pinv(C) @ C - C‖_max = {recon_err:.2e}")
# Verify rank
from eigensystems import matrix_rank
print(f"matrix_rank(C) = {matrix_rank(C)} (expected 2)")
# ------------------------------------------------------------------
# Plot: four heatmaps — A_raw, A_sym, V_raw columns, Q columns
# ------------------------------------------------------------------
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_matrix_utils.svg"
fig, axes = plt.subplots(2, 2, figsize=(13, 10))
cmap = "RdBu_r"
def _heatmap(ax: object, data: np.ndarray, title: str) -> None:
"""Draw a labelled heatmap on ax.
Args:
ax: Matplotlib Axes object.
data (np.ndarray): 2-D array to visualise.
title (str): Axes title.
"""
im = ax.imshow(data, cmap=cmap, aspect="auto")
fig.colorbar(im, ax=ax, fraction=0.046)
ax.set_title(title)
ax.set_xlabel("Column")
ax.set_ylabel("Row")
for i in range(data.shape[0]):
for j in range(data.shape[1]):
val = data[i, j]
ax.text(j, i, f"{val:.2f}", ha="center", va="center",
fontsize=7,
color="white" if abs(val) > 0.5 * np.max(np.abs(data)) else "black")
_heatmap(axes[0, 0], A_raw, "Before SYMMETRIZE\n(asymmetric A_raw)")
_heatmap(axes[0, 1], A_sym, "After SYMMETRIZE\n(A + A^T)/2")
_heatmap(axes[1, 0], V_raw, "Before ORTHOGONALIZE\n(random 6×4 V)")
_heatmap(axes[1, 1], Q, "After ORTHOGONALIZE\n(Q — orthonormal columns)")
fig.suptitle("Matrix Utility Functions — IMSL Routines", fontsize=13,
fontweight="bold")
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
print(f"\nPlot saved: {plot_path}")
return {
"sym_result": A_sym,
"Q": Q,
"pinv_A": C_pinv,
"orth_err": float(orth_err),
"pinv_err": float(recon_err),
"plot_path": str(plot_path),
}
if __name__ == "__main__":
run_demo_matrix_utils()
Plot Output¶
Top row: asymmetric matrix A (left) and its symmetrized form (right). Bottom row: random 6×4 matrix V (left) and orthonormal Q (right).¶
Console Output¶
IMSL IS_SYMMETRIC / IS_POSITIVE_DEFINITE Checks
-------------------------------------------------------
S_sym (3×3 SPD) is_symmetric=True is_positive_definite=True
S_asym (upper tri) is_symmetric=False is_positive_definite=True
S_indef (indefinite) is_symmetric=True is_positive_definite=False
S_herm (Hermitian) is_symmetric=True is_positive_definite=True
IMSL SYMMETRIZE Example
Original (asymmetric):
[[1. 4. 7.]
[2. 5. 8.]
[3. 6. 9.]]
Symmetrized (A + A^T)/2:
[[1. 3. 5.]
[3. 5. 7.]
[5. 7. 9.]]
Max asymmetry after symmetrize: 0.00e+00 (should be ~0)
is_symmetric(symmetrize(A_raw)) = True
IMSL ORTHOGONALIZE Example
Input V shape: (6, 4) → Q shape: (6, 4)
Q^T Q Γëê I, max error = 5.55e-16 (should be ~0)
IMSL PSEUDO_INVERSE Example
C shape: (3, 4), pseudo_inverse shape: (4, 3)
Reconstruction error ΓÇûC @ pinv(C) @ C - CΓÇû_max = 1.33e-15
matrix_rank(C) = 2 (expected 2)
Plot saved: test_output\example_imsl_matrix_utils.svg