How to Calculate the Inverse of a Matrix in NumPy

Key Insights

Use numpy.linalg.inv() for square, non-singular matrices, but prefer numpy.linalg.solve() when solving linear systems—it’s faster and more numerically stable.
The pseudo-inverse numpy.linalg.pinv() handles singular and non-square matrices where the true inverse doesn’t exist, making it essential for least-squares problems.
Always check for singular matrices before inversion using determinants or try/except blocks—silent numerical errors from near-singular matrices can corrupt your results.

Introduction to Matrix Inversion

Matrix inversion is a fundamental operation in linear algebra that shows up constantly in scientific computing, machine learning, and data analysis. The inverse of a matrix A, denoted A⁻¹, satisfies the property A × A⁻¹ = I, where I is the identity matrix. Think of it as the matrix equivalent of dividing by a number.

You’ll encounter matrix inversion when solving systems of linear equations, computing covariance matrices in statistics, implementing linear regression, or performing coordinate transformations in graphics. Understanding how to compute inverses efficiently—and when to avoid computing them entirely—is essential knowledge for any Python developer working with numerical data.

This article assumes you have NumPy installed (pip install numpy) and understand basic matrix concepts like dimensions and multiplication.

Using numpy.linalg.inv() — The Standard Approach

The numpy.linalg.inv() function is your go-to tool for computing matrix inverses. It accepts a square matrix and returns its inverse using LU decomposition under the hood.

import numpy as np

# Create a simple 2x2 matrix
A = np.array([[4, 7],
              [2, 6]])

# Compute the inverse
A_inv = np.linalg.inv(A)

print("Original matrix A:")
print(A)
print("\nInverse matrix A⁻¹:")
print(A_inv)

# Verify: A @ A_inv should equal the identity matrix
identity = A @ A_inv
print("\nA @ A⁻¹ (should be identity):")
print(identity)

Output:

Original matrix A:
[[4 7]
 [2 6]]

Inverse matrix A⁻¹:
[[ 0.6 -0.7]
 [-0.2  0.4]]

A @ A⁻¹ (should be identity):
[[1.00000000e+00 0.00000000e+00]
 [1.11022302e-16 1.00000000e+00]]

Notice the tiny floating-point error (1.11e-16) instead of exact zeros. This is normal and unavoidable in numerical computing. For practical purposes, use np.allclose() to verify inverses:

# Proper way to verify an inverse
A_3x3 = np.array([[1, 2, 3],
                  [0, 1, 4],
                  [5, 6, 0]])

A_3x3_inv = np.linalg.inv(A_3x3)

# Check if product is close to identity
is_valid = np.allclose(A_3x3 @ A_3x3_inv, np.eye(3))
print(f"Inverse is valid: {is_valid}")  # True

Handling Singular Matrices

Not every matrix has an inverse. A matrix is singular (non-invertible) when its determinant equals zero, which happens when rows or columns are linearly dependent. Attempting to invert a singular matrix raises a LinAlgError.

import numpy as np
from numpy.linalg import LinAlgError

# This matrix is singular - row 2 is 2x row 1
singular_matrix = np.array([[1, 2],
                            [2, 4]])

# Check the determinant first
det = np.linalg.det(singular_matrix)
print(f"Determinant: {det}")  # 0.0 (or very close to 0)

# Safe inversion with exception handling
def safe_inverse(matrix):
    """Compute matrix inverse with proper error handling."""
    try:
        # Check determinant first for better error messages
        det = np.linalg.det(matrix)
        if np.abs(det) < 1e-10:
            print(f"Warning: Matrix is near-singular (det={det:.2e})")
            return None
        
        return np.linalg.inv(matrix)
    
    except LinAlgError as e:
        print(f"Cannot invert matrix: {e}")
        return None

# Test with singular matrix
result = safe_inverse(singular_matrix)
print(f"Inverse result: {result}")  # None

# Test with valid matrix
valid_matrix = np.array([[1, 2],
                         [3, 4]])
result = safe_inverse(valid_matrix)
print(f"Valid inverse computed: {result is not None}")  # True

The determinant check catches near-singular matrices that might technically invert but produce wildly inaccurate results. A determinant close to zero (relative to the matrix scale) signals numerical instability.

Alternative: Moore-Penrose Pseudo-Inverse with numpy.linalg.pinv()

The true inverse only exists for square, non-singular matrices. For everything else—rectangular matrices, singular matrices, or overdetermined systems—use the Moore-Penrose pseudo-inverse.

import numpy as np

# Rectangular matrix (3x2) - no true inverse exists
rectangular = np.array([[1, 2],
                        [3, 4],
                        [5, 6]])

# This would fail:
# np.linalg.inv(rectangular)  # LinAlgError: Last 2 dimensions must be square

# Pseudo-inverse works fine
pseudo_inv = np.linalg.pinv(rectangular)
print(f"Original shape: {rectangular.shape}")      # (3, 2)
print(f"Pseudo-inverse shape: {pseudo_inv.shape}") # (2, 3)

print("\nPseudo-inverse:")
print(pseudo_inv)

The pseudo-inverse has different properties than the true inverse. For a matrix A with pseudo-inverse A⁺:

A × A⁺ × A = A (but A × A⁺ ≠ I in general)
A⁺ minimizes ||Ax - b||² for least-squares problems

# Compare inv() and pinv() on a square matrix
square = np.array([[1, 2],
                   [3, 4]])

true_inv = np.linalg.inv(square)
pseudo_inv = np.linalg.pinv(square)

print("True inverse:")
print(true_inv)
print("\nPseudo-inverse:")
print(pseudo_inv)
print(f"\nAre they equal? {np.allclose(true_inv, pseudo_inv)}")  # True

# For singular matrices, only pinv() works
singular = np.array([[1, 2],
                     [2, 4]])

pseudo_inv_singular = np.linalg.pinv(singular)
print("\nPseudo-inverse of singular matrix:")
print(pseudo_inv_singular)

Use pinv() when you’re unsure about matrix properties or when working with data matrices that might be rank-deficient.

Performance Considerations

Matrix inversion has O(n³) time complexity—doubling the matrix size increases computation time roughly eightfold. More importantly, explicit inversion is often unnecessary and wasteful.

When solving Ax = b for x, you might think to compute x = A⁻¹b. Don’t. Use numpy.linalg.solve() instead:

import numpy as np
import time

def benchmark_methods(n, num_trials=100):
    """Compare inversion vs solve for linear systems."""
    # Create random system Ax = b
    np.random.seed(42)
    A = np.random.randn(n, n)
    b = np.random.randn(n)
    
    # Method 1: Explicit inversion
    start = time.perf_counter()
    for _ in range(num_trials):
        A_inv = np.linalg.inv(A)
        x_inv = A_inv @ b
    time_inv = time.perf_counter() - start
    
    # Method 2: Direct solve
    start = time.perf_counter()
    for _ in range(num_trials):
        x_solve = np.linalg.solve(A, b)
    time_solve = time.perf_counter() - start
    
    # Verify both give same answer
    assert np.allclose(x_inv, x_solve)
    
    return time_inv, time_solve

# Benchmark different sizes
for n in [50, 100, 200, 500]:
    t_inv, t_solve = benchmark_methods(n)
    speedup = t_inv / t_solve
    print(f"n={n:3d}: inv+matmul={t_inv:.4f}s, solve={t_solve:.4f}s, speedup={speedup:.1f}x")

Typical output:

n= 50: inv+matmul=0.0089s, solve=0.0051s, speedup=1.7x
n=100: inv+matmul=0.0198s, solve=0.0082s, speedup=2.4x
n=200: inv+matmul=0.0712s, solve=0.0241s, speedup=3.0x
n=500: inv+matmul=0.5765s, solve=0.1842s, speedup=3.1x

The solve() function uses LU decomposition without explicitly forming the inverse, saving both time and memory while improving numerical stability.

Practical Application: Solving Linear Systems

Let’s tie everything together with a real example: fitting a polynomial to data points using linear algebra.

import numpy as np

# Data points to fit
x_data = np.array([0, 1, 2, 3, 4])
y_data = np.array([1, 2.5, 3.8, 8.1, 15.2])

# Fit a quadratic: y = a + bx + cx²
# This creates the Vandermonde matrix
degree = 2
A = np.vander(x_data, degree + 1, increasing=True)
print("Design matrix A:")
print(A)

# Method 1: Using explicit inverse (NOT recommended)
# For overdetermined systems, use normal equations: A^T A x = A^T b
ATA = A.T @ A
ATb = A.T @ y_data

coeffs_inv = np.linalg.inv(ATA) @ ATb
print(f"\nCoefficients via inverse: {coeffs_inv}")

# Method 2: Using solve (better)
coeffs_solve = np.linalg.solve(ATA, ATb)
print(f"Coefficients via solve:  {coeffs_solve}")

# Method 3: Using lstsq (best for this problem)
coeffs_lstsq, residuals, rank, s = np.linalg.lstsq(A, y_data, rcond=None)
print(f"Coefficients via lstsq:  {coeffs_lstsq}")

# All methods give the same answer
print(f"\nAll methods agree: {np.allclose(coeffs_inv, coeffs_solve)}")

# Predict values
x_test = np.linspace(0, 4, 50)
y_pred = coeffs_lstsq[0] + coeffs_lstsq[1]*x_test + coeffs_lstsq[2]*x_test**2
print(f"\nFitted polynomial: y = {coeffs_lstsq[0]:.2f} + {coeffs_lstsq[1]:.2f}x + {coeffs_lstsq[2]:.2f}x²")

For least-squares problems like this, np.linalg.lstsq() is the right choice—it handles rank-deficient matrices and provides useful diagnostics.

Summary and Best Practices

Here’s your quick reference for matrix inversion in NumPy:

Situation	Function	Notes
Square, non-singular matrix	`np.linalg.inv()`	Verify with determinant check
Solving Ax = b	`np.linalg.solve()`	Faster and more stable than inv()
Non-square or singular matrix	`np.linalg.pinv()`	Returns pseudo-inverse
Least-squares problems	`np.linalg.lstsq()`	Best for overdetermined systems

Common pitfalls to avoid:

Computing explicit inverses for linear systems. Use solve() instead—it’s 2-3x faster and more accurate.
Ignoring numerical stability. Near-singular matrices can produce garbage results without raising errors. Always check condition numbers with np.linalg.cond() for critical applications.
Assuming exact results. Floating-point arithmetic means A @ A_inv won’t exactly equal I. Use np.allclose() for comparisons.
Using inv() on rectangular matrices. It won’t work. Use pinv() or lstsq() depending on your goal.

Matrix inversion is a powerful tool, but it’s often the wrong tool. When you reach for inv(), ask yourself: “Am I solving a linear system?” If yes, use solve(). Your code will be faster, more accurate, and more robust.