Skip to main content

Python NumPy: How to Calculate Matrix Inverse in NumPy

The matrix inverse A⁻¹ satisfies A × A⁻¹ = I (identity matrix). NumPy provides efficient methods for computing inverses, though understanding when to use alternatives is equally important.

Basic Inverse Calculation​

Use np.linalg.inv() to compute the inverse of a square matrix:

import numpy as np

A = np.array([
    [1, 2],
    [3, 4]
])

A_inv = np.linalg.inv(A)
print("Inverse:")
print(A_inv)

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

Output:

Inverse:
[[-2.   1. ]
 [ 1.5 -0.5]]

A @ A_inv:
[[1.0000000e+00 0.0000000e+00]
 [8.8817842e-16 1.0000000e+00]

Verification with np.allclose()​

Due to floating-point precision, use np.allclose() for verification:

import numpy as np

A = np.array([[1, 2], [3, 4]])
A_inv = np.linalg.inv(A)

# Check if product equals identity
is_identity = np.allclose(A @ A_inv, np.eye(2))
print(f"A @ A_inv ≈ I: {is_identity}")  # True

# Also verify A_inv @ A = I
is_identity_2 = np.allclose(A_inv @ A, np.eye(2))
print(f"A_inv @ A ≈ I: {is_identity_2}")  # True

Output:

A @ A_inv ≈ I: True
A_inv @ A ≈ I: True

Handling Singular Matrices​

A singular matrix (determinant = 0) has no inverse:

import numpy as np

# Singular matrix: second row is 2× first row
singular = np.array([
    [1, 2],
    [2, 4]
])

# Check determinant first
det = np.linalg.det(singular)
print(f"Determinant: {det}")  # 0.0

np.linalg.inv(singular) # Error: Singular matrix

Output:

Determinant: 0.0
numpy.linalg.LinAlgError: Singular matrix

Safe Inverse Function​

import numpy as np

def safe_inverse(matrix, tolerance=1e-10):
    """Compute inverse with singularity check."""
    det = np.linalg.det(matrix)
    
    if np.abs(det) < tolerance:
        print(f"Warning: Matrix is singular or near-singular (det={det:.2e})")
        return None
    
    return np.linalg.inv(matrix)

# Test
A = np.array([[1, 2], [3, 4]])
print(safe_inverse(A))

B = np.array([[1, 2], [2, 4]])
print(safe_inverse(B))  # Returns None with warning

Output:

[[-2.   1. ]
 [ 1.5 -0.5]]
Warning: Matrix is singular or near-singular (det=0.00e+00)
None

Don't Use inv() to Solve Linear Systems​

For solving Ax = b, use np.linalg.solve() instead of computing the inverse:

import numpy as np

A = np.array([[1, 2], [3, 4]])
b = np.array([5, 11])

# ❌ Slow and less accurate
x_slow = np.linalg.inv(A) @ b
print(f"Using inv(): {x_slow}")

# ✅ Faster and more stable
x_fast = np.linalg.solve(A, b)
print(f"Using solve(): {x_fast}")

# Verify both solutions
print(f"\nVerification (A @ x = b):")
print(f"inv method: {A @ x_slow}")
print(f"solve method: {A @ x_fast}")

Output:

Using inv(): [1. 2.]
Using solve(): [1. 2.]

Verification (A @ x = b):
inv method: [ 5. 11.]
solve method: [ 5. 11.]
warning

Computing the inverse explicitly is slower and introduces more numerical error than using solve(). Only compute the inverse when you specifically need the inverse matrix itself.

Performance Comparison​

import numpy as np
import timeit

n = 500
A = np.random.randn(n, n)
b = np.random.randn(n)

# Ensure A is invertible
A = A + n * np.eye(n)

# Method 1: Using inverse
def solve_with_inv():
    return np.linalg.inv(A) @ b

# Method 2: Using solve
def solve_direct():
    return np.linalg.solve(A, b)

inv_time = timeit.timeit(solve_with_inv, number=10)
solve_time = timeit.timeit(solve_direct, number=10)

print(f"inv() method: {inv_time:.4f}s")
print(f"solve() method: {solve_time:.4f}s")
print(f"solve() is {inv_time/solve_time:.1f}x faster")

Output:

inv() method: 17.3293s
solve() method: 0.1554s
solve() is 111.5x faster

Pseudo-Inverse for Non-Square Matrices​

The Moore-Penrose pseudo-inverse works for any matrix shape:

import numpy as np

# Non-square matrix (3×2)
A = np.array([
    [1, 2],
    [3, 4],
    [5, 6]
])

# Regular inverse fails
try:
    np.linalg.inv(A)
except np.linalg.LinAlgError as e:
    print(f"inv() error: {e}")

# Pseudo-inverse works
A_pinv = np.linalg.pinv(A)
print(f"Shape of A: {A.shape}")             # (3, 2)
print(f"Shape of pinv(A): {A_pinv.shape}")  # (2, 3)
print(A_pinv)

Output:

inv() error: Last 2 dimensions of the array must be square
Shape of A: (3, 2)
Shape of pinv(A): (2, 3)
[[-1.33333333 -0.33333333  0.66666667]
 [ 1.08333333  0.33333333 -0.41666667]]

Pseudo-Inverse for Least Squares​

import numpy as np

# Overdetermined system (more equations than unknowns)
A = np.array([
    [1, 1],
    [1, 2],
    [1, 3]
])
b = np.array([1, 2, 2])

# Least squares solution using pseudo-inverse
x = np.linalg.pinv(A) @ b
print(f"Least squares solution: {x}")

# Equivalent using lstsq
x_lstsq, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
print(f"lstsq solution: {x_lstsq}")

Output:

Least squares solution: [0.66666667 0.5       ]
lstsq solution: [0.66666667 0.5       ]

Inverse of Special Matrices​

Diagonal Matrix​

import numpy as np

# Diagonal matrix: inverse is reciprocal of diagonal
D = np.diag([2, 4, 5])
D_inv = np.linalg.inv(D)

print("D:")
print(D)
print("\nD_inv (1/diagonal):")
print(D_inv)

Output:

D:
[[2 0 0]
 [0 4 0]
 [0 0 5]]

D_inv (1/diagonal):
[[0.5  0.   0.  ]
 [0.   0.25 0.  ]
 [0.   0.   0.2 ]]

Orthogonal Matrix​

import numpy as np

# Orthogonal matrix: inverse equals transpose
Q, _ = np.linalg.qr(np.random.randn(3, 3))

Q_inv = np.linalg.inv(Q)
Q_T = Q.T

print("Q_inv ≈ Q.T:", np.allclose(Q_inv, Q_T))  # True

# For orthogonal matrices, use transpose instead of inv()

Output:

Q_inv ≈ Q.T: True
tip

For orthogonal matrices, the inverse equals the transpose. Use Q.T instead of np.linalg.inv(Q) for better performance.

Checking Invertibility​

Multiple ways to check if a matrix is invertible:

import numpy as np

def check_invertibility(matrix):
    """Check if a matrix is invertible using multiple methods."""
    n = matrix.shape[0]
    
    # Method 1: Determinant
    det = np.linalg.det(matrix)
    print(f"Determinant: {det:.6f}")
    
    # Method 2: Rank
    rank = np.linalg.matrix_rank(matrix)
    print(f"Rank: {rank} (need {n} for invertibility)")
    
    # Method 3: Condition number
    cond = np.linalg.cond(matrix)
    print(f"Condition number: {cond:.2e}")
    
    # Verdict
    is_invertible = not np.isclose(det, 0) and rank == n
    print(f"Invertible: {is_invertible}")
    
    return is_invertible

# Test
A = np.array([[1, 2], [3, 4]])
check_invertibility(A)

print()

B = np.array([[1, 2], [2, 4]])
check_invertibility(B)

Output:

Determinant: -2.000000
Rank: 2 (need 2 for invertibility)
Condition number: 1.49e+01
Invertible: True

Determinant: 0.000000
Rank: 1 (need 2 for invertibility)
Condition number: 2.52e+16
Invertible: False

Batch Inverse (Multiple Matrices)​

Invert multiple matrices efficiently:

import numpy as np

# Stack of 3 matrices, each 2×2
matrices = np.array([
    [[1, 2], [3, 4]],
    [[2, 0], [0, 2]],
    [[1, 1], [0, 1]]
])

print(f"Shape: {matrices.shape}")  # (3, 2, 2)

# Invert all at once
inverses = np.linalg.inv(matrices)
print(f"Inverses shape: {inverses.shape}")  # (3, 2, 2)

# Verify first matrix
print(np.allclose(matrices[0] @ inverses[0], np.eye(2)))  # True

Output:

Shape: (3, 2, 2)
Inverses shape: (3, 2, 2)
True

Summary​

MethodUse CaseMatrix Shape
np.linalg.inv(A)True inverse neededSquare, invertible
np.linalg.pinv(A)Pseudo-inverse, least squaresAny shape
np.linalg.solve(A, b)Solving Ax = bSquare
A.TOrthogonal matrix inverseOrthogonal
CheckPurpose
np.linalg.det(A) ≠ 0Invertibility
np.linalg.cond(A)Numerical stability
np.linalg.matrix_rank(A)Full rank check

Best Practice: Use np.linalg.solve() for linear systems-it's faster and more numerically stable. Only compute the explicit inverse when you need to apply it multiple times or analyze the inverse matrix itself.

Table of Contents