NumPy - Linear Algebra (np.linalg) Overview

Key Insights

NumPy’s linalg module provides production-ready implementations of fundamental linear algebra operations including matrix decomposition, eigenvalue computation, and solving linear systems with performance optimized through BLAS/LAPACK libraries
Understanding the distinction between @ operator for matrix multiplication, np.dot() for various products, and specialized functions like np.inner() and np.outer() prevents common calculation errors in numerical computing
Matrix decomposition methods (SVD, QR, Cholesky) enable efficient solutions to overdetermined systems, least squares problems, and numerical stability improvements over direct matrix inversion

Matrix Operations and Products

NumPy distinguishes between element-wise and matrix operations. The @ operator and np.matmul() perform matrix multiplication, while * performs element-wise multiplication.

import numpy as np

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

# Matrix multiplication
result_matmul = A @ B
print("Matrix multiplication:\n", result_matmul)
# [[19 22]
#  [43 50]]

# Element-wise multiplication
result_elementwise = A * B
print("Element-wise:\n", result_elementwise)
# [[ 5 12]
#  [21 32]]

# Dot product for 1D arrays
v1 = np.array([1, 2, 3])
v2 = np.array([4, 5, 6])
dot_product = np.dot(v1, v2)  # Returns scalar: 32

# Inner and outer products
inner = np.inner(v1, v2)  # Same as dot for 1D: 32
outer = np.outer(v1, v2)
print("Outer product:\n", outer)
# [[ 4  5  6]
#  [ 8 10 12]
#  [12 15 18]]

Solving Linear Systems

The np.linalg.solve() function solves systems of linear equations Ax = b more efficiently and accurately than computing the inverse matrix.

# Solve Ax = b
A = np.array([[3, 1], [1, 2]])
b = np.array([9, 8])

x = np.linalg.solve(A, b)
print("Solution:", x)  # [2. 3.]

# Verify solution
print("Verification:", np.allclose(A @ x, b))  # True

# For multiple right-hand sides
B = np.array([[9, 1], [8, 2]])
X = np.linalg.solve(A, B)
print("Multiple solutions:\n", X)

# Least squares for overdetermined systems
# More equations than unknowns
A_overdetermined = np.array([[1, 1], [1, 2], [1, 3]])
b_overdetermined = np.array([2, 3, 4.5])

# Use lstsq instead of solve
result = np.linalg.lstsq(A_overdetermined, b_overdetermined, rcond=None)
x_lstsq = result[0]
residuals = result[1]
print("Least squares solution:", x_lstsq)  # [0.83333333 1.        ]
print("Residuals:", residuals)  # [0.08333333]

Matrix Decompositions

Decomposition methods factorize matrices into simpler components, enabling efficient computation and improved numerical stability.

# Singular Value Decomposition (SVD)
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]])
U, s, Vt = np.linalg.svd(A, full_matrices=False)

print("Shape of U:", U.shape)    # (4, 3)
print("Singular values:", s)      # [25.46... 1.29... 0.0...]
print("Shape of Vt:", Vt.shape)   # (3, 3)

# Reconstruct original matrix
S = np.diag(s)
A_reconstructed = U @ S @ Vt
print("Reconstruction error:", np.linalg.norm(A - A_reconstructed))

# QR decomposition for orthogonalization
A = np.array([[1, 2], [3, 4], [5, 6]], dtype=float)
Q, R = np.linalg.qr(A)

print("Q (orthogonal):\n", Q)
print("R (upper triangular):\n", R)
print("Q is orthogonal:", np.allclose(Q.T @ Q, np.eye(2)))

# Cholesky decomposition for positive definite matrices
A_pd = np.array([[4, 12, -16], [12, 37, -43], [-16, -43, 98]])
L = np.linalg.cholesky(A_pd)
print("Lower triangular L:\n", L)
print("Verification:", np.allclose(L @ L.T, A_pd))

Eigenvalues and Eigenvectors

Eigendecomposition reveals fundamental properties of linear transformations and appears frequently in dimensionality reduction, stability analysis, and graph algorithms.

# Standard eigenvalue problem
A = np.array([[1, 2], [2, 1]])
eigenvalues, eigenvectors = np.linalg.eig(A)

print("Eigenvalues:", eigenvalues)      # [ 3. -1.]
print("Eigenvectors:\n", eigenvectors)

# Verify Av = λv
for i in range(len(eigenvalues)):
    lam = eigenvalues[i]
    v = eigenvectors[:, i]
    print(f"Av = {A @ v}, λv = {lam * v}")

# For symmetric/Hermitian matrices, use eigh (more efficient)
A_symmetric = np.array([[1, 2, 3], [2, 4, 5], [3, 5, 6]])
eigenvalues_sym, eigenvectors_sym = np.linalg.eigh(A_symmetric)
print("Symmetric eigenvalues:", eigenvalues_sym)

# Matrix power using eigendecomposition
def matrix_power_eigen(A, n):
    eigenvalues, eigenvectors = np.linalg.eig(A)
    D_n = np.diag(eigenvalues ** n)
    return eigenvectors @ D_n @ np.linalg.inv(eigenvectors)

A = np.array([[2, 1], [1, 2]])
A_10 = matrix_power_eigen(A, 10)
print("A^10 via eigen:\n", A_10)
print("A^10 via numpy:\n", np.linalg.matrix_power(A, 10))

Matrix Properties and Norms

Understanding matrix properties helps assess numerical stability and condition of problems.

# Matrix norms
A = np.array([[1, 2, 3], [4, 5, 6]])

frobenius_norm = np.linalg.norm(A)  # Default: Frobenius
print("Frobenius norm:", frobenius_norm)

# Different norm types
l1_norm = np.linalg.norm(A, ord=1)          # Max column sum
l2_norm = np.linalg.norm(A, ord=2)          # Spectral norm
linf_norm = np.linalg.norm(A, ord=np.inf)   # Max row sum
print(f"L1: {l1_norm}, L2: {l2_norm}, L∞: {linf_norm}")

# Condition number (measures numerical stability)
A_square = np.array([[1, 2], [3, 4]])
cond = np.linalg.cond(A_square)
print("Condition number:", cond)

# Ill-conditioned matrix example
A_ill = np.array([[1, 1], [1, 1.0001]])
print("Ill-conditioned:", np.linalg.cond(A_ill))  # Very large

# Matrix rank
A_rank_deficient = np.array([[1, 2, 3], [2, 4, 6], [1, 1, 1]])
rank = np.linalg.matrix_rank(A_rank_deficient)
print("Rank:", rank)  # 2, not full rank

# Determinant and trace
A = np.array([[1, 2], [3, 4]])
det = np.linalg.det(A)
trace = np.trace(A)  # Sum of diagonal elements
print(f"Determinant: {det}, Trace: {trace}")

Practical Applications

Real-world scenarios demonstrate when to apply specific linear algebra operations.

# Principal Component Analysis (PCA) using SVD
def pca_transform(X, n_components):
    """Dimensionality reduction using SVD."""
    # Center the data
    X_centered = X - np.mean(X, axis=0)
    
    # Compute SVD
    U, s, Vt = np.linalg.svd(X_centered, full_matrices=False)
    
    # Transform data
    X_transformed = U[:, :n_components] * s[:n_components]
    
    return X_transformed, Vt[:n_components]

# Example data
X = np.random.randn(100, 5)
X_pca, components = pca_transform(X, n_components=2)
print("Original shape:", X.shape)
print("Transformed shape:", X_pca.shape)

# Solving polynomial regression with normal equations
def polynomial_fit(x, y, degree):
    """Fit polynomial using least squares."""
    # Create Vandermonde matrix
    A = np.vander(x, degree + 1, increasing=True)
    
    # Solve normal equations: A^T A c = A^T y
    coeffs = np.linalg.solve(A.T @ A, A.T @ y)
    
    return coeffs

x = np.linspace(0, 10, 50)
y = 2 * x**2 + 3 * x + 1 + np.random.randn(50) * 5
coeffs = polynomial_fit(x, y, degree=2)
print("Polynomial coefficients:", coeffs)

# Distance matrix computation for clustering
points = np.random.randn(100, 3)
# Efficient distance matrix using broadcasting
diff = points[:, np.newaxis, :] - points[np.newaxis, :, :]
distances = np.linalg.norm(diff, axis=2)
print("Distance matrix shape:", distances.shape)

Performance Considerations

NumPy’s linear algebra operations link to optimized BLAS/LAPACK libraries, but choosing the right function matters.

import time

# Avoid explicit matrix inversion
A = np.random.randn(1000, 1000)
b = np.random.randn(1000)

# Slow: compute inverse explicitly
start = time.time()
x_inv = np.linalg.inv(A) @ b
time_inv = time.time() - start

# Fast: use solve directly
start = time.time()
x_solve = np.linalg.solve(A, b)
time_solve = time.time() - start

print(f"Inverse method: {time_inv:.4f}s")
print(f"Solve method: {time_solve:.4f}s")
print(f"Speedup: {time_inv/time_solve:.2f}x")

# Use specialized functions for matrix structure
A_symmetric = A @ A.T  # Guaranteed symmetric positive definite

# For symmetric systems, consider using Cholesky
L = np.linalg.cholesky(A_symmetric)
y = np.linalg.solve(L, b)
x_cholesky = np.linalg.solve(L.T, y)

The np.linalg module provides industrial-strength implementations suitable for production systems. Always prefer specialized solvers over generic approaches, validate numerical stability through condition numbers, and leverage decomposition methods for both performance and accuracy.