Python, simultaneous pseudo-inversion of many 3x3, singular, symmetric, matrices

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

EDIT: See @Jaime's answer. Only the discussion in the comments to this answer is useful now, and only for the specific problem at hand.

You can do this matrix by matrix, using scipy, that provides pinv (link) to calculate the Moore-Penrose pseudo inverse. An example follows:

from scipy.linalg import det,eig,pinv
from numpy.random import randint
#generate a random singular matrix M first
while True:
    M = randint(0,10,9).reshape(3,3)
    if det(M)==0:
        break
M = M.astype(float)
#this is the method you need
MPpseudoinverse = pinv(M)

This does not exploit the fact that the matrix is symmetric though. You may also want to try the version of pinv exposed by numpy, that is supposedely faster, and different. See this post.