NumPy 1.8 included linear algebra gufuncs, which do exactly what you are after. While np.linalg.pinv
is not gufunc-ed, np.linalg.svd
is, and behind the scenes that is the function that gets called. So you can define your own gupinv
function, based on the source code of the original function, as follows:
def gu_pinv(a, rcond=1e-15):
a = np.asarray(a)
swap = np.arange(a.ndim)
swap[[-2, -1]] = swap[[-1, -2]]
u, s, v = np.linalg.svd(a)
cutoff = np.maximum.reduce(s, axis=-1, keepdims=True) * rcond
mask = s > cutoff
s[mask] = 1. / s[mask]
s[~mask] = 0
return np.einsum('...uv,...vw->...uw',
np.transpose(v, swap) * s[..., None, :],
np.transpose(u, swap))
And you can now do things like:
a = np.random.rand(50, 40, 30, 6)
b = np.empty(a.shape[:-1] + (3, 3), dtype=a.dtype)
# Expand the unique items into a full symmetrical matrix
b[..., 0, :] = a[..., :3]
b[..., 1:, 0] = a[..., 1:3]
b[..., 1, 1:] = a[..., 3:5]
b[..., 2, 1:] = a[..., 4:]
# make matrix at [1, 2, 3] singular
b[1, 2, 3, 2] = b[1, 2, 3, 0] + b[1, 2, 3, 1]
# Find all the pseudo-inverses
pi = gu_pinv(b)
And of course the results are correct, both for singular and non-singular matrices:
>>> np.allclose(pi[0, 0, 0], np.linalg.pinv(b[0, 0, 0]))
True
>>> np.allclose(pi[1, 2, 3], np.linalg.pinv(b[1, 2, 3]))
True
And for this example, with 50 * 40 * 30 = 60,000
pseudo-inverses calculated:
In [2]: %timeit pi = gu_pinv(b)
1 loops, best of 3: 422 ms per loop
Which is really not that bad, although it is noticeably (4x) slower than simply calling np.linalg.inv
, but this of course fails to properly handle the singular arrays:
In [8]: %timeit np.linalg.inv(b)
10 loops, best of 3: 98.8 ms per loop