How to Find the Null Space of a Matrix in Python

Key Insights

The null space of a matrix A contains all vectors x where Ax = 0, revealing linear dependencies and constraint relationships in your data
Use scipy.linalg.null_space() for production code—it handles numerical precision automatically and returns orthonormal basis vectors
Always verify your null space vectors by computing Ax and checking the result is near-zero within floating-point tolerance (typically 1e-10)

Introduction to the Null Space

The null space (or kernel) of a matrix A is the set of all vectors x that satisfy Ax = 0. While this sounds abstract, it’s fundamental to understanding linear systems, data dependencies, and constraint satisfaction problems.

In practical terms, null space vectors represent combinations of your variables that produce zero output. If you’re modeling physical constraints, these vectors show you which combinations satisfy homogeneous constraints. In data science, they reveal linear dependencies between features. In computer graphics, they help solve inverse kinematics and constraint systems.

The dimension of the null space (called nullity) relates directly to the rank through the rank-nullity theorem: rank(A) + nullity(A) = number of columns. A full-rank matrix has only the trivial null space (the zero vector), while a rank-deficient matrix has a non-trivial null space that reveals redundancies in your system.

Mathematical Foundation

Before diving into code, understand that finding the null space involves identifying which linear combinations of columns sum to zero. Row reduction to Reduced Row Echelon Form (RREF) is the classical approach taught in linear algebra courses, but it’s numerically unstable for computer implementation.

The rank-nullity theorem tells us that for an m×n matrix A:

rank(A) + nullity(A) = n
If rank(A) = n, the null space contains only the zero vector
If rank(A) < n, there are n - rank(A) linearly independent null space vectors

Let’s start with a simple example where we can visualize the null space:

import numpy as np
from scipy.linalg import null_space

# Create a simple 3x3 rank-deficient matrix
# Third row is sum of first two rows
A = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [5, 7, 9]
])

print(f"Matrix A:\n{A}")
print(f"\nRank: {np.linalg.matrix_rank(A)}")
print(f"Expected nullity: {A.shape[1] - np.linalg.matrix_rank(A)}")

# Find null space
ns = null_space(A)
print(f"\nNull space basis:\n{ns}")

# Verify: A @ ns should be approximately zero
result = A @ ns
print(f"\nA @ null_space:\n{result}")
print(f"Maximum absolute value: {np.max(np.abs(result))}")

This matrix has rank 2, so we expect a 1-dimensional null space. The null space vector shows us the linear dependency: the third row equals the sum of the first two.

Finding Null Space with NumPy

While NumPy doesn’t have a dedicated null space function, you can extract it from Singular Value Decomposition (SVD). The SVD decomposes A = UΣV^T, where the columns of V corresponding to zero (or near-zero) singular values form the null space basis.

Here’s a complete implementation:

import numpy as np

def null_space_svd(A, rcond=None):
    """
    Compute null space using SVD.
    
    Parameters:
    -----------
    A : array_like
        Input matrix
    rcond : float, optional
        Relative condition number. Singular values smaller than
        rcond * largest_singular_value are considered zero.
    """
    # Perform SVD
    u, s, vh = np.linalg.svd(A, full_matrices=True)
    
    # Determine tolerance for zero singular values
    if rcond is None:
        rcond = np.finfo(A.dtype).eps * max(A.shape)
    
    tol = np.max(s) * rcond
    
    # Find number of non-zero singular values
    num_nonzero = np.sum(s > tol)
    
    # Null space is spanned by rows of vh after the rank
    null_space = vh[num_nonzero:].T.conj()
    
    return null_space

# Test with a 4x6 matrix
A = np.array([
    [1, 2, 3, 4, 5, 6],
    [2, 4, 6, 8, 10, 12],
    [1, 1, 1, 1, 1, 1],
    [0, 1, 2, 3, 4, 5]
])

print(f"Matrix shape: {A.shape}")
print(f"Rank: {np.linalg.matrix_rank(A)}")

ns = null_space_svd(A)
print(f"\nNull space shape: {ns.shape}")
print(f"Null space basis vectors:\n{ns}")

# Verify orthonormality
print(f"\nOrthonormality check (should be identity):\n{ns.T @ ns}")

The SVD approach is numerically stable and gives you orthonormal basis vectors. The rcond parameter controls the threshold for treating singular values as zero—critical for numerical stability with floating-point arithmetic.

Using SciPy’s Null Space Function

SciPy provides scipy.linalg.null_space(), which implements the SVD approach with sensible defaults. This is your go-to function for production code:

from scipy.linalg import null_space
import numpy as np

# Same 4x6 matrix
A = np.array([
    [1, 2, 3, 4, 5, 6],
    [2, 4, 6, 8, 10, 12],
    [1, 1, 1, 1, 1, 1],
    [0, 1, 2, 3, 4, 5]
])

# Find null space with default tolerance
ns_scipy = null_space(A)

print(f"Null space dimensions: {ns_scipy.shape}")
print(f"Null space vectors:\n{ns_scipy}")

# Custom tolerance for stricter zero detection
ns_strict = null_space(A, rcond=1e-10)
print(f"\nWith stricter tolerance: {ns_strict.shape}")

# Verify
verification = A @ ns_scipy
print(f"\nVerification (A @ null_space):")
print(f"Max absolute error: {np.max(np.abs(verification))}")

The rcond parameter works identically to the SVD implementation above. For most applications, the default value works well, but you may need to adjust it for ill-conditioned matrices or when you know the expected numerical precision.

Verification and Validation

Always verify your null space computation. Floating-point errors mean Ax won’t be exactly zero, but it should be very close:

import numpy as np
from scipy.linalg import null_space

def verify_null_space(A, ns, tol=1e-10):
    """
    Verify that ns contains valid null space vectors for matrix A.
    
    Returns True if all checks pass, raises AssertionError otherwise.
    """
    # Check 1: A @ ns should be approximately zero
    result = A @ ns
    max_error = np.max(np.abs(result))
    
    print(f"Maximum |A @ ns|: {max_error:.2e}")
    assert max_error < tol, f"Null space verification failed: {max_error} >= {tol}"
    
    # Check 2: Null space vectors should be orthonormal
    gram = ns.T @ ns
    identity = np.eye(ns.shape[1])
    ortho_error = np.max(np.abs(gram - identity))
    
    print(f"Orthonormality error: {ortho_error:.2e}")
    assert ortho_error < tol, f"Vectors not orthonormal: {ortho_error} >= {tol}"
    
    # Check 3: Dimension check via rank-nullity theorem
    expected_nullity = A.shape[1] - np.linalg.matrix_rank(A)
    actual_nullity = ns.shape[1]
    
    print(f"Expected nullity: {expected_nullity}, Actual: {actual_nullity}")
    assert expected_nullity == actual_nullity, "Nullity doesn't match rank-nullity theorem"
    
    return True

# Example usage
A = np.random.rand(5, 8)
A[3] = A[0] + A[1]  # Create dependency
A[4] = 2 * A[2]     # Another dependency

ns = null_space(A)
verify_null_space(A, ns)
print("✓ All verification checks passed")

This verification function should be part of your testing suite whenever you compute null spaces, especially with real-world data where numerical issues are common.

Practical Applications

Here’s a real-world example: solving a constraint satisfaction problem in robotics or computer graphics. Suppose you have a system where certain joint angles must satisfy linear constraints:

import numpy as np
from scipy.linalg import null_space

# Example: 3D point cloud with known distance constraints
# We want to find configurations that satisfy: sum of coordinates = 0
# and x + y - z = 0

# Constraint matrix (each row is a constraint)
C = np.array([
    [1, 1, 1],      # x + y + z = 0
    [1, 1, -1]      # x + y - z = 0
])

# Find null space (valid configurations)
valid_configs = null_space(C)

print(f"Constraint-satisfying configurations:\n{valid_configs}")

# Generate specific solutions
# Any linear combination of null space vectors satisfies constraints
solution1 = valid_configs @ np.array([1.0])
solution2 = valid_configs @ np.array([2.5])

print(f"\nSolution 1: {solution1}")
print(f"Verification: {C @ solution1}")

print(f"\nSolution 2: {solution2}")
print(f"Verification: {C @ solution2}")

This pattern appears in optimization problems, computer vision (epipolar geometry), and data analysis (finding feature combinations that explain zero variance).

Performance Considerations and Edge Cases

Different matrix types require different strategies:

import numpy as np
from scipy.linalg import null_space
from scipy.sparse.linalg import svds
import time

# Full-rank matrix (trivial null space)
A_full_rank = np.random.rand(5, 5)
ns_full = null_space(A_full_rank)
print(f"Full-rank matrix null space shape: {ns_full.shape}")
# Should be (5, 0) - empty null space

# Rank-deficient matrix
A_deficient = np.random.rand(5, 10)
A_deficient[:, 5:] = A_deficient[:, :5]  # Copy first 5 columns

start = time.time()
ns_dense = null_space(A_deficient)
dense_time = time.time() - start

print(f"\nDense method time: {dense_time:.4f}s")
print(f"Null space dimension: {ns_dense.shape[1]}")

# For large sparse matrices, consider sparse methods
from scipy.sparse import csr_matrix

# Note: scipy.sparse doesn't have null_space, but you can use SVD
A_sparse = csr_matrix(A_deficient)
# For sparse matrices, use iterative methods or convert to dense
# if the matrix is small enough

# Numerical stability test
A_ill_conditioned = np.array([
    [1, 1, 1],
    [1, 1.0001, 1],
    [1, 1, 1.0001]
])

ns_default = null_space(A_ill_conditioned)
ns_strict = null_space(A_ill_conditioned, rcond=1e-3)

print(f"\nIll-conditioned matrix:")
print(f"Default tolerance null space: {ns_default.shape}")
print(f"Strict tolerance null space: {ns_strict.shape}")

For large sparse matrices, converting to dense for null space computation may be impractical. Consider whether you actually need the full null space or just need to solve Ax = 0 for specific cases, where iterative methods like conjugate gradient might be more appropriate.

The key takeaway: scipy.linalg.null_space() is robust and efficient for most use cases. Use SVD-based methods for numerical stability, always verify your results, and adjust tolerance parameters when working with ill-conditioned matrices. Understanding the null space gives you powerful tools for solving constraint systems, identifying redundancies, and understanding the structure of linear transformations in your data.