Matrix Analysis Functions — DET, TRACE, RANK, FROBENIUS, SVD_VALUES, SCHUR

This example applies every scalar matrix analysis and decomposition routine to a 4×4 test matrix and a rank-deficient 3×3 matrix:

Example Code

"""IMSL matrix analysis examples: determinant, trace, matrix_rank,
frobenius_norm, svd_values, schur_decomp.

Demonstrates:
- Scalar properties of a 4×4 matrix (det, trace, rank, Frobenius norm)
- schur_decomp: A = Z T Z^H verification
- svd_values on a rank-deficient matrix vs full svd
- Heatmaps of A and T (Schur form) + singular value bar chart + scalar
  property bar chart

Outputs:
- Printed scalar properties and verification errors
- SVG figure saved to test_output/example_imsl_matrix_analysis.svg
- Console output captured to docs/_generated_text/matrix_analysis_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 (
    determinant,
    frobenius_norm,
    matrix_rank,
    schur_decomp,
    svd,
    svd_values,
    trace,
)


def run_demo_matrix_analysis() -> dict:
    """Run all matrix analysis demonstrations.

    Creates a 4×4 matrix and a rank-deficient matrix, then applies every
    scalar analysis and decomposition function.

    Args:
        None

    Returns:
        dict: Keys ``det``, ``tr``, ``rank``, ``fnorm``, ``singular_values``,
            ``schur_T``, ``schur_Z``, ``plot_path`` (str).
    """
    # ------------------------------------------------------------------
    # 4×4 test matrix
    # ------------------------------------------------------------------
    A = np.array([
        [ 4.0,  3.0, -2.0,  1.0],
        [ 2.0,  5.0,  1.0,  0.0],
        [-1.0,  2.0,  6.0,  3.0],
        [ 0.0,  1.0, -1.0,  4.0],
    ])

    det_A   = determinant(A)
    trace_A = trace(A)
    rank_A  = matrix_rank(A)
    fnorm_A = frobenius_norm(A)

    print("\nIMSL Matrix Analysis: 4×4 Test Matrix")
    print("A =")
    print(A)
    print("-" * 50)
    print(f"determinant(A)     = {det_A:.6f}")
    print(f"trace(A)           = {trace_A:.6f}")
    print(f"matrix_rank(A)     = {rank_A}")
    print(f"frobenius_norm(A)  = {fnorm_A:.6f}")

    # Cross-check: trace equals sum of eigenvalues; det equals product
    eigs_A = np.linalg.eigvals(A)
    trace_check = np.sum(eigs_A).real
    det_check   = np.prod(eigs_A).real
    print(f"Trace check (sum of eigenvalues)  = {trace_check:.6f}")
    print(f"Det check   (prod of eigenvalues) = {det_check:.6f}")

    # ------------------------------------------------------------------
    # Schur decomposition: A = Z T Z^H
    # ------------------------------------------------------------------
    T, Z = schur_decomp(A)
    recon_err = np.max(np.abs(A - Z @ T @ Z.T.conj()))

    print("\nIMSL SCHUR_DECOMP Example")
    print("-" * 50)
    print("Schur form T (quasi-upper-triangular):")
    print(np.round(T, 4))
    print(f"Reconstruction error ‖A - Z T Z^H‖_max = {recon_err:.2e}")

    # ------------------------------------------------------------------
    # SVD — full vs values-only, including rank-deficient matrix
    # ------------------------------------------------------------------
    # Rank-deficient: third row is sum of first two
    B = np.array([
        [1.0, 2.0, 3.0],
        [4.0, 5.0, 6.0],
        [5.0, 7.0, 9.0],  # row3 = row1 + row2
    ])

    s_only  = svd_values(B)
    U, s_full, Vh = svd(B)
    rank_B  = matrix_rank(B)

    print("\nIMSL SVD_VALUES Example: Rank-Deficient 3×3 Matrix")
    print("B (row3 = row1 + row2, expected rank 2):")
    print(B)
    print("-" * 50)
    print(f"Singular values (svd_values): {s_only}")
    print(f"Singular values (svd full):   {s_full}")
    print(f"matrix_rank(B) = {rank_B}  (expected 2)")
    recon_err_b = np.max(np.abs(B - U @ np.diag(s_full) @ Vh))
    print(f"SVD reconstruction error ‖B - U Σ V^H‖_max = {recon_err_b:.2e}")

    # ------------------------------------------------------------------
    # Plot
    # ------------------------------------------------------------------
    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_matrix_analysis.svg"

    fig, axes = plt.subplots(2, 2, figsize=(13, 10))

    cmap = "RdBu_r"

    # Heatmap of A
    ax = axes[0, 0]
    im = ax.imshow(A, cmap=cmap, aspect="auto")
    fig.colorbar(im, ax=ax, fraction=0.046)
    ax.set_title("Matrix A (4×4)")
    ax.set_xlabel("Column")
    ax.set_ylabel("Row")
    for i in range(4):
        for j in range(4):
            ax.text(j, i, f"{A[i, j]:.1f}", ha="center", va="center",
                    fontsize=9, color="white" if abs(A[i, j]) > 3 else "black")

    # Heatmap of Schur form T
    ax = axes[0, 1]
    im = ax.imshow(T.real, cmap=cmap, aspect="auto")
    fig.colorbar(im, ax=ax, fraction=0.046)
    ax.set_title("Schur Form T  (A = Z T Z^H)")
    ax.set_xlabel("Column")
    ax.set_ylabel("Row")
    for i in range(4):
        for j in range(4):
            val = T.real[i, j]
            if abs(val) > 0.001:
                ax.text(j, i, f"{val:.2f}", ha="center", va="center",
                        fontsize=8, color="white" if abs(val) > 3 else "black")

    # Singular value bar chart of B
    ax = axes[1, 0]
    ax.bar(np.arange(1, len(s_only) + 1), s_only, color="#0369a1",
           alpha=0.85, edgecolor="#1e3a5f")
    for k, sv in enumerate(s_only):
        ax.text(k + 1, sv + 0.1, f"{sv:.3f}", ha="center", fontsize=9)
    ax.set_xticks(np.arange(1, len(s_only) + 1))
    ax.set_xticklabels([f{k+1}" for k in range(len(s_only))])
    ax.set_ylabel("Singular value")
    ax.set_title("SVD_VALUES: Rank-Deficient Matrix B\n(σ₃ ≈ 0 → rank 2)")

    # Scalar properties bar chart
    ax = axes[1, 1]
    labels = ["det(A)", "trace(A)", "rank(A)", "‖A‖_F"]
    values = [float(det_A), float(trace_A), float(rank_A), float(fnorm_A)]
    colors_bar = ["#be123c", "#059669", "#d97706", "#7c3aed"]
    bars = ax.bar(labels, values, color=colors_bar, alpha=0.85, edgecolor="white")
    for bar, v in zip(bars, values):
        ax.text(bar.get_x() + bar.get_width() / 2,
                bar.get_height() + 0.3 * np.sign(v),
                f"{v:.2f}", ha="center", va="bottom" if v >= 0 else "top",
                fontsize=9)
    ax.axhline(0, color="gray", linewidth=0.6)
    ax.set_title("Scalar Properties of A")
    ax.set_ylabel("Value")

    fig.suptitle("Matrix Analysis 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 {
        "det": det_A,
        "tr": trace_A,
        "rank": rank_A,
        "fnorm": fnorm_A,
        "singular_values": s_only,
        "schur_T": T,
        "schur_Z": Z,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_matrix_analysis()

Plot Output

Heatmaps of A and Schur form T, singular value bar chart, scalar properties bar chart

Top-left: heatmap of the test matrix A. Top-right: Schur form T showing quasi-upper-triangular structure. Bottom-left: singular values of the rank-deficient matrix B (σ₃ ≈ 0). Bottom-right: scalar properties of A.

Console Output

IMSL Matrix Analysis: 4×4 Test Matrix
A =
[[ 4.  3. -2.  1.]
 [ 2.  5.  1.  0.]
 [-1.  2.  6.  3.]
 [ 0.  1. -1.  4.]]
--------------------------------------------------
determinant(A)     = 308.000000
trace(A)           = 19.000000
matrix_rank(A)     = 4
frobenius_norm(A)  = 11.313708
Trace check (sum of eigenvalues)  = 19.000000
Det check   (prod of eigenvalues) = 308.000000

IMSL SCHUR_DECOMP Example
--------------------------------------------------
Schur form T (quasi-upper-triangular):
[[ 1.7994  1.0329 -2.0627 -1.7783]
 [ 0.      6.718   1.4487  2.9363]
 [ 0.     -0.2323  6.718  -1.0337]
 [ 0.      0.      0.      3.7647]]
Reconstruction error ΓÇûA - Z T Z^HΓÇû_max = 7.99e-15

IMSL SVD_VALUES Example: Rank-Deficient 3×3 Matrix
B (row3 = row1 + row2, expected rank 2):
[[1. 2. 3.]
 [4. 5. 6.]
 [5. 7. 9.]]
--------------------------------------------------
Singular values (svd_values): [1.56633231e+01 8.12593979e-01 1.17807579e-15]
Singular values (svd full):   [1.56633231e+01 8.12593979e-01 1.17807579e-15]
matrix_rank(B) = 2  (expected 2)
SVD reconstruction error ‖B - U Σ V^H‖_max = 4.44e-15

Plot saved: test_output\example_imsl_matrix_analysis.svg