Compute the Moore-Penrose pseudoinverse of a stack of matrices in Python

The Moore-Penrose pseudoinverse is a generalization of the matrix inverse for non-square or singular matrices. In NumPy, you can compute the pseudoinverse of a stack of matrices using numpy.linalg.pinv(), which uses singular value decomposition (SVD) internally.

Syntax

numpy.linalg.pinv(a, rcond=1e-15, hermitian=False)

Parameters

The function accepts the following parameters:

• a ? Matrix or stack of matrices to be pseudo-inverted
• rcond ? Cutoff for small singular values. Values ? rcond × largest_singular_value are set to zero
• hermitian ? If True, assumes the matrix is Hermitian for more efficient computation

Example

Let's compute the pseudoinverse of a stack of 2×2 matrices ?

import numpy as np

# Create a stack of three 2x2 matrices
arr = np.array([
    [[1, 2], [3, 4]], 
    [[1, 2], [2, 1]], 
    [[1, 3], [3, 1]]
])

print("Original stack of matrices:")
print(arr)
print("\nShape:", arr.shape)
print("Dimensions:", arr.ndim)

# Compute the Moore-Penrose pseudoinverse
pseudoinverse = np.linalg.pinv(arr)
print("\nPseudoinverse:")
print(pseudoinverse)

Original stack of matrices:
[[[1 2]
  [3 4]]

 [[1 2]
  [2 1]]

 [[1 3]
  [3 1]]]

Shape: (3, 2, 2)
Dimensions: 3

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

 [[-0.33333333  0.66666667]
  [ 0.66666667 -0.33333333]]

 [[-0.125  0.375]
  [ 0.375 -0.125]]]

Verification

We can verify that A × A? × A = A (where A? is the pseudoinverse) ?

import numpy as np

# Original matrices
arr = np.array([[[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]]])

# Compute pseudoinverse
pinv_arr = np.linalg.pinv(arr)

# Verify: A @ A+ @ A should equal A
verification = arr @ pinv_arr @ arr

print("Original matrices:")
print(arr[0])  # First matrix
print("\nVerification (A @ A+ @ A):")
print(verification[0])  # Should match original
print("\nAre they equal?", np.allclose(arr[0], verification[0]))

Original matrices:
[[1 2]
 [3 4]]

Verification (A @ A+ @ A):
[[1. 2.]
 [3. 4.]]

Are they equal? True

Using rcond Parameter

The rcond parameter controls numerical precision by setting small singular values to zero ?

import numpy as np

# Create a nearly singular matrix
arr = np.array([[[1, 2], [2, 4.0001]]])

# Default rcond
pinv1 = np.linalg.pinv(arr)
print("Default rcond result:")
print(pinv1[0])

# Higher rcond (more aggressive cutoff)
pinv2 = np.linalg.pinv(arr, rcond=1e-3)
print("\nWith rcond=1e-3:")
print(pinv2[0])

Default rcond result:
[[ 2.49996875e+04 -1.24993750e+04]
 [-4.99987500e+04  2.49993750e+04]]

With rcond=1e-3:
[[0.2 0.2]
 [0.1 0.1]]

Conclusion

Use numpy.linalg.pinv() to compute the Moore-Penrose pseudoinverse of matrix stacks. The rcond parameter helps handle numerical stability, while hermitian=True can speed up computation for symmetric matrices.