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

Question

I have a 3D image with row sizes x cols x deps. For each voxel in the image, I calculated a 3x3 real symmetric matrix. They are stored in an array D, which therefore has the form (rows, columns, prints, 6).

D stores 6 unique 3x3 symmetric matrix elements for each voxel in my image. I need to find the Moore-Penrose pseudo version for all rows of * cols * deps matrices at the same time / in vectorized code (looping through each voxel image and inverting is too slow in Python).

Some of these 3x3 symmetric matrices are nonsingular, and I can find their inverse in the vectorized code using the analytical formula for the true inversion of the nonsingular symmetric 3x3 matrix, and I did it.

However, for those matrices that are unique (and there will certainly be some), I need the Moore-Penrose pseudo-reverse. I could get the analytic formula for the MP of a real, singular symmetric 3x3 matrix, but this is really an unpleasant / long formula and therefore will include a VERY large number of (elementary) matrix arithmetic and quite a bit confusing looking code.

Therefore, I would like to know if there is a way to simultaneously find MP pseudo-inverse for all of these matrices immediately numerically. Is there any way to do this?

thankful gf

Answer 1

NumPy 1.8 gufuncs, , . np.linalg.pinv gufunc-ed, np.linalg.svd , - , . , gupinv :

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))

, :

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)

, , , , :

>>> 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

50 * 40 * 30 = 60,000 -:

In [2]: %timeit pi = gu_pinv(b)
1 loops, best of 3: 422 ms per loop

, (4x) , np.linalg.inv, , , :

In [8]: %timeit np.linalg.inv(b)
10 loops, best of 3: 98.8 ms per loop

Answer 2

EDIT: . @Jaime. .

, scipy, pinv (), Moore-Penrose . :

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)

, . pinv, numpy, . . .