Python, одновременное псевдообращение множества сингулярных симметричных матриц 3x3

NumPy 1.8 включает элементы линейной алгебры, которые делают именно то, что вам нужно. В то время как 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

Что на самом деле не так уж и плохо, хотя это заметно (в 4 раза) медленнее, чем простой вызов np.linalg.inv, но это, конечно, не позволяет должным образом обрабатывать отдельные массивы:

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

РЕДАКТИРОВАТЬ: см. Ответ @Jaime. Только обсуждение в комментариях к этому ответу сейчас полезно и только для конкретной проблемы.

Вы можете сделать эту матрицу по матрице, используя scipy, который предоставляет pinv (link) для вычисления псевдообратного преобразования Мура-Пенроуза. Пример ниже:

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, которая предположительно быстрее и отличается. См. это сообщение.