Chapter 159 — Determinants: Computation

# --- Determinants: Cofactor expansion and LU method ---
import numpy as np

def det_cofactor(A):
    """
    Compute determinant via cofactor expansion along first row.
    Recursive, O(n!) — illustrative only, not for large matrices.

    Args:
        A: square 2D numpy array

    Returns:
        scalar determinant value
    """
    n = A.shape[0]
    if n == 1: return A[0, 0]
    if n == 2: return A[0,0]*A[1,1] - A[0,1]*A[1,0]
    det = 0.0
    for j in range(n):
        # Minor: delete row 0, column j
        minor = np.delete(np.delete(A, 0, axis=0), j, axis=1)
        cofactor = ((-1)**j) * det_cofactor(minor)
        det += A[0, j] * cofactor
    return det

def det_via_lu(A):
    """
    Compute determinant via LU decomposition.
    O(n³) and numerically stable.

    Args:
        A: square 2D numpy array

    Returns:
        scalar determinant value
    """
    # PA = LU; det(A) = det(P)⁻¹ * det(L) * det(U)
    # det(L) = product of diagonal (all 1s for unit lower triangular)
    # det(U) = product of diagonal entries
    # det(P) = (-1)^(number of swaps)
    import scipy.linalg
    P, L, U = scipy.linalg.lu(A)
    # Count permutation sign from P
    p_vec = np.argmax(P, axis=0)   # where each column sends the basis
    sign = 1
    visited = [False] * len(p_vec)
    for i in range(len(p_vec)):
        if not visited[i]:
            cycle_len = 0
            j = i
            while not visited[j]:
                visited[j] = True
                j = p_vec[j]
                cycle_len += 1
            if cycle_len % 2 == 0:
                sign *= -1
    return sign * np.prod(np.diag(U))

# Test on small matrices
test_cases = [
    np.array([[2.,1.],[1.,3.]]),
    np.array([[1.,2.,3.],[4.,5.,6.],[7.,8.,9.]]),   # singular!
    np.array([[2.,1.,0.],[1.,3.,1.],[0.,1.,2.]]),
]
for A in test_cases:
    d_cofactor = det_cofactor(A)
    d_numpy    = np.linalg.det(A)
    print(f"n={A.shape[0]}: cofactor={d_cofactor:.4f}, numpy={d_numpy:.4f}, match={np.isclose(d_cofactor, d_numpy)}")

# Characteristic polynomial: roots of det(A - λI) = 0 are eigenvalues
print("\n--- Characteristic polynomial ---")
A = np.array([[4.,1.],[2.,3.]])
lambdas = np.linspace(0, 6, 500)
char_poly = [det_cofactor(A - l*np.eye(2)) for l in lambdas]
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')
fig, ax = plt.subplots(figsize=(8,4))
ax.plot(lambdas, char_poly, 'steelblue', lw=2)
ax.axhline(0, color='k', lw=1)
true_eigs = np.linalg.eigvals(A)
for e in true_eigs:
    ax.axvline(e, color='crimson', linestyle='--', label=f'λ={e:.2f}')
ax.set_xlabel('λ'); ax.set_ylabel('det(A - λI)')
ax.set_title('Characteristic Polynomial: Roots = Eigenvalues')
ax.legend(); plt.tight_layout(); plt.show()