I have a 3D image with dimensions rows x cols x deps. For every voxel in the image, I have computed a 3x3 real symmetric matrix. They are stored in the array D, which therefore has shape (rows, cols, deps, 6).

D stores the 6 unique elements of the 3x3 symmetric matrix for every voxel in my image. I need to find the Moore-Penrose pseudo inverse of all row\*cols\*deps matrices simultaneously/in vectorized code (looping through every image voxel and inverting is far too slow in Python).

Some of these 3x3 symmetric matrices are non-singular, and I can find their inverses, in vectorized code, using the analytical formula for the true inverse of a non-singular 3x3 symmetric matrix, and I've done that.

However, for those matrices that ARE singular (and there are sure to be some) I need the Moore-Penrose pseudo inverse. I could derive an analytical formula for the MP of a real, singular, symmetric 3x3 matrix, but it's a really nasty/lengthy formula, and would therefore involve a VERY large number of (element-wise) matrix arithmetic and quite a bit of confusing looking code.

Hence, I would like to know if there is a way to simultaneously find the MP pseudo inverse for all these matrices at once numerically. Is there a way to do this?

Gratefully,
GF