Chapter 157 — Matrix Inverse

# --- Matrix Inverse: Definition, Computation, Pitfalls ---
import numpy as np

def inverse_gauss_jordan(A):
    """
    Compute matrix inverse via Gauss-Jordan elimination.
    Augments [A | I] and row-reduces to [I | A⁻¹].

    Args:
        A: square 2D numpy array, shape (n,n)

    Returns:
        A_inv: 2D numpy array, shape (n,n)

    Raises:
        ValueError if A is singular.
    """
    n = A.shape[0]
    assert A.shape == (n, n), "Must be square"
    # Augmented matrix [A | I]
    aug = np.hstack([A.astype(float).copy(), np.eye(n)])
    for col in range(n):
        # Find pivot
        pivot_row = col + np.argmax(np.abs(aug[col:, col]))
        if abs(aug[pivot_row, col]) < 1e-12:
            raise ValueError("Matrix is singular or nearly singular")
        aug[[col, pivot_row]] = aug[[pivot_row, col]]   # swap rows
        # Normalize pivot row
        aug[col] /= aug[col, col]
        # Eliminate in all other rows
        for row in range(n):
            if row != col:
                aug[row] -= aug[row, col] * aug[col]
    return aug[:, n:]   # right half is A⁻¹

# Test
np.random.seed(42)
A = np.array([[2., 1., 0.],
              [1., 3., 1.],
              [0., 1., 2.]])

A_inv_scratch = inverse_gauss_jordan(A)
A_inv_numpy   = np.linalg.inv(A)

print("Inverse via Gauss-Jordan:")
print(np.round(A_inv_scratch, 6))
print("\nMax diff from numpy: ", np.max(np.abs(A_inv_scratch - A_inv_numpy)))
print("\nVerify AA⁻¹ = I:", np.allclose(A @ A_inv_scratch, np.eye(3)))

# When NOT to invert: solve instead
b = np.array([1., 2., 3.])
# Bad: x = A⁻¹ @ b  (computes full inverse, then multiplies)
x_inv = np.linalg.inv(A) @ b
# Good: x = np.linalg.solve(A, b)  (uses LU decomposition, faster and more stable)
x_solve = np.linalg.solve(A, b)
print(f"\nSolve via inv: {x_inv}")
print(f"Solve via solve: {x_solve}")
print(f"Equal: {np.allclose(x_inv, x_solve)}")
print("\nRule: Never compute A⁻¹ explicitly if you only need A⁻¹b.")

# --- Pseudoinverse for non-square systems ---
import numpy as np

# Overdetermined system: more equations than unknowns
# A is (5,3) — 5 equations, 3 unknowns
A = np.array([[1,0,0],[0,1,0],[0,0,1],[1,1,0],[0,1,1]], dtype=float)
b = np.array([1., 2., 3., 3.5, 4.5])   # inconsistent — no exact solution

# Pseudoinverse gives least-squares solution
A_pinv = np.linalg.pinv(A)
x_lstsq = A_pinv @ b

print(f"A shape: {A.shape}, A⁺ shape: {A_pinv.shape}")
print(f"Least-squares solution: {x_lstsq.round(4)}")
print(f"Residual ||Ax - b||: {np.linalg.norm(A @ x_lstsq - b):.4f}")

# Underdetermined: fewer equations than unknowns (infinitely many solutions)
# pinv gives minimum-norm solution
A_under = np.array([[1,1,0],[0,1,1]], dtype=float)  # 2 eq, 3 unknowns
b_under = np.array([3., 2.])
x_min_norm = np.linalg.pinv(A_under) @ b_under
print(f"\nUnderdetermined min-norm solution: {x_min_norm.round(4)}")
print(f"||x||: {np.linalg.norm(x_min_norm):.4f}")