How to Compute the Pseudoinverse in Python

Key Insights

The Moore-Penrose pseudoinverse generalizes matrix inversion to non-square and singular matrices, making it essential for solving overdetermined and underdetermined linear systems.
NumPy’s linalg.pinv() is the standard implementation for computing pseudoinverses in Python, using SVD decomposition with configurable tolerance for numerical stability.
While convenient for solving linear systems, direct pseudoinverse computation is O(n³) and numerically expensive—consider QR decomposition or iterative methods for large-scale problems.

Introduction to the Pseudoinverse

The Moore-Penrose pseudoinverse extends the concept of matrix inversion to matrices that don’t have a regular inverse. While a regular inverse exists only for square, non-singular matrices, the pseudoinverse is defined for any matrix—rectangular, singular, or otherwise.

Mathematically, for a matrix A, its pseudoinverse A⁺ satisfies four conditions: AA⁺A = A, A⁺AA⁺ = A⁺, (AA⁺)ᵀ = AA⁺, and (A⁺A)ᵀ = A⁺A. When A is invertible, A⁺ equals A⁻¹.

You need the pseudoinverse when working with:

Non-square matrices (more rows than columns or vice versa)
Singular square matrices (determinant equals zero)
Overdetermined systems (more equations than unknowns)
Underdetermined systems (fewer equations than unknowns)

Here’s what happens when you try to invert a singular matrix:

import numpy as np

# Create a singular matrix (determinant = 0)
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

print(f"Determinant: {np.linalg.det(A):.10f}")

try:
    A_inv = np.linalg.inv(A)
except np.linalg.LinAlgError as e:
    print(f"Error: {e}")
    
# Pseudoinverse works fine
A_pinv = np.linalg.pinv(A)
print(f"\nPseudoinverse computed successfully")
print(f"Shape: {A_pinv.shape}")

Computing Pseudoinverse with NumPy

NumPy provides numpy.linalg.pinv() as the standard function for computing the Moore-Penrose pseudoinverse. It uses Singular Value Decomposition (SVD) internally, which is numerically stable and works for any matrix.

The function signature is: numpy.linalg.pinv(a, rcond=1e-15, hermitian=False). The rcond parameter sets the cutoff ratio for small singular values—any singular value smaller than rcond * largest_singular_value is treated as zero.

import numpy as np

# Create a rectangular matrix
A = np.array([[1, 2],
              [3, 4],
              [5, 6]])

# Compute pseudoinverse
A_pinv = np.linalg.pinv(A)

print(f"Original matrix A:\n{A}")
print(f"\nPseudoinverse A+:\n{A_pinv}")
print(f"\nShape A: {A.shape}, Shape A+: {A_pinv.shape}")

# Verify the property: A @ A+ @ A = A
verification = A @ A_pinv @ A
print(f"\nVerification (A @ A+ @ A):\n{verification}")
print(f"Close to original? {np.allclose(A, verification)}")

# Also verify: A+ @ A @ A+ = A+
verification2 = A_pinv @ A @ A_pinv
print(f"\nA+ @ A @ A+ equals A+? {np.allclose(A_pinv, verification2)}")

Computing Pseudoinverse with SciPy

SciPy offers an alternative implementation via scipy.linalg.pinv(). The main difference is the check_finite parameter, which can be set to False to skip input validation for better performance when you’re certain your data is well-formed.

import numpy as np
from scipy import linalg as scipy_linalg
from numpy import linalg as numpy_linalg
import time

# Create a test matrix
np.random.seed(42)
A = np.random.randn(100, 50)

# NumPy version
start = time.perf_counter()
A_pinv_numpy = numpy_linalg.pinv(A)
numpy_time = time.perf_counter() - start

# SciPy version with input checking
start = time.perf_counter()
A_pinv_scipy = scipy_linalg.pinv(A)
scipy_time = time.perf_counter() - start

# SciPy version without input checking
start = time.perf_counter()
A_pinv_scipy_fast = scipy_linalg.pinv(A, check_finite=False)
scipy_fast_time = time.perf_counter() - start

print(f"NumPy time: {numpy_time*1000:.3f} ms")
print(f"SciPy time (with check): {scipy_time*1000:.3f} ms")
print(f"SciPy time (no check): {scipy_fast_time*1000:.3f} ms")
print(f"\nResults match? {np.allclose(A_pinv_numpy, A_pinv_scipy)}")

For most applications, NumPy’s implementation is sufficient and more commonly used. SciPy’s version is useful when you need the additional control or are already using SciPy for other linear algebra operations.

Practical Application: Solving Linear Systems

The pseudoinverse excels at solving linear systems Ax = b when A isn’t invertible. For overdetermined systems (more equations than unknowns), x = A⁺b gives the least-squares solution that minimizes ||Ax - b||.

import numpy as np

# Overdetermined system: 5 equations, 3 unknowns
# Represents fitting a plane to 5 points
A = np.array([[1, 1, 1],
              [1, 2, 1],
              [1, 3, 1],
              [1, 4, 1],
              [1, 5, 1]])

b = np.array([2.1, 3.9, 6.2, 7.8, 10.1])

# Solve using pseudoinverse
A_pinv = np.linalg.pinv(A)
x = A_pinv @ b

print(f"Solution x: {x}")

# Verify the solution
b_predicted = A @ x
residual = b - b_predicted

print(f"\nOriginal b: {b}")
print(f"Predicted b: {b_predicted}")
print(f"Residual: {residual}")
print(f"Residual norm: {np.linalg.norm(residual):.6f}")

# Compare with numpy's least squares solver
x_lstsq, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
print(f"\nSolutions match? {np.allclose(x, x_lstsq)}")

Practical Application: Linear Regression

Linear regression can be solved elegantly using the pseudoinverse. The normal equations β = (XᵀX)⁻¹Xᵀy simplify to β = X⁺y when using the pseudoinverse.

import numpy as np
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt

# Generate synthetic data
np.random.seed(42)
X = np.random.randn(100, 1) * 10
y = 3 * X.squeeze() + 7 + np.random.randn(100) * 5

# Add bias term (column of ones)
X_with_bias = np.column_stack([np.ones(len(X)), X])

# Solve using pseudoinverse
X_pinv = np.linalg.pinv(X_with_bias)
beta = X_pinv @ y

print(f"Coefficients (pseudoinverse): intercept={beta[0]:.3f}, slope={beta[1]:.3f}")

# Compare with sklearn
model = LinearRegression()
model.fit(X, y)

print(f"Coefficients (sklearn): intercept={model.intercept_:.3f}, slope={model.coef_[0]:.3f}")
print(f"\nCoefficients match? {np.allclose([model.intercept_, model.coef_[0]], beta)}")

# Make predictions
y_pred_pinv = X_with_bias @ beta
y_pred_sklearn = model.predict(X)

# Calculate R² score manually
ss_res = np.sum((y - y_pred_pinv) ** 2)
ss_tot = np.sum((y - np.mean(y)) ** 2)
r2 = 1 - (ss_res / ss_tot)

print(f"R² score: {r2:.4f}")

Performance and Numerical Stability

Computing the pseudoinverse has O(mn²) complexity for an m×n matrix where m ≥ n, dominated by the SVD computation. For large matrices, this becomes computationally expensive.

The rcond parameter controls numerical stability. Small singular values amplify noise, so they’re zeroed out. The default 1e-15 works for most cases, but ill-conditioned matrices may need adjustment.

import numpy as np
import time

# Benchmark different matrix sizes
sizes = [100, 200, 500, 1000]
times = []

for n in sizes:
    A = np.random.randn(n, n)
    
    start = time.perf_counter()
    A_pinv = np.linalg.pinv(A)
    elapsed = time.perf_counter() - start
    
    times.append(elapsed)
    print(f"Size {n}x{n}: {elapsed:.4f} seconds")

# Demonstrate ill-conditioned matrix
print("\n--- Ill-conditioned Matrix Example ---")
# Create matrix with very small singular values
U = np.random.randn(5, 5)
S = np.diag([1000, 100, 10, 0.01, 0.0001])  # Wide range of singular values
V = np.random.randn(5, 5)
A_ill = U @ S @ V.T

# Different rcond values
for rcond in [1e-15, 1e-10, 1e-5, 1e-2]:
    A_pinv = np.linalg.pinv(A_ill, rcond=rcond)
    reconstruction = A_ill @ A_pinv @ A_ill
    error = np.linalg.norm(A_ill - reconstruction)
    print(f"rcond={rcond:.0e}: reconstruction error = {error:.6e}")

# Alternative: Use SVD directly for more control
print("\n--- Direct SVD Approach ---")
U, s, Vt = np.linalg.svd(A_ill, full_matrices=False)
print(f"Singular values: {s}")

# Manually construct pseudoinverse with custom threshold
threshold = 1e-5
s_inv = np.where(s > threshold, 1/s, 0)
A_pinv_manual = Vt.T @ np.diag(s_inv) @ U.T
print(f"Manual pseudoinverse matches pinv? {np.allclose(A_pinv_manual, np.linalg.pinv(A_ill, rcond=1e-5))}")

For large-scale problems, consider alternatives:

QR decomposition: More efficient for least-squares problems
Iterative methods: LSQR or LSMR for sparse matrices
Direct solvers: np.linalg.lstsq() is optimized for solving systems

Conclusion

The pseudoinverse is an indispensable tool for solving linear systems beyond the constraints of regular matrix inversion. NumPy’s linalg.pinv() provides a reliable, SVD-based implementation suitable for most applications.

Use the pseudoinverse when you need to:

Solve overdetermined or underdetermined linear systems
Implement linear regression or other statistical models from scratch
Work with singular or rectangular matrices

However, remember that pseudoinverse computation is expensive. For production code with large matrices, profile your application and consider specialized solvers. The rcond parameter deserves attention when dealing with ill-conditioned matrices—don’t accept the default blindly.

The pseudoinverse elegantly generalizes matrix inversion, but like all numerical methods, understanding its limitations and computational characteristics ensures you apply it appropriately.