Python，同时伪逆多个3x3，奇异，对称矩阵

Question

我有一张三维图像，尺寸是行 x 列 x 深度。对于图像中的每一个体素（也就是三维像素），我计算了一个3x3的实对称矩阵。这些矩阵存储在一个叫D的数组里，所以D的形状是（行，列，深度，6）。

D存储了每个体素对应的3x3对称矩阵的6个独特元素。我需要同时找到所有行*列*深度矩阵的Moore-Penrose伪逆，或者说用向量化的代码来处理（因为逐个体素循环计算反转在Python中太慢了）。

其中一些3x3的对称矩阵是非奇异的，我可以使用非奇异3x3对称矩阵的真实反转的公式，通过向量化的代码找到它们的逆，我已经做到这一点了。

但是，对于那些是奇异的矩阵（肯定会有一些），我需要Moore-Penrose伪逆。我可以推导出一个实奇异对称3x3矩阵的MP的解析公式，但这个公式非常复杂且冗长，因此会涉及大量的（逐元素）矩阵运算和看起来很混乱的代码。

所以，我想知道有没有办法同时以数值方式找到所有这些矩阵的MP伪逆。有没有这样的办法呢？

感激不尽，GF

Answer 1

编辑：请查看@Jaime的回答。现在只有这个回答下的评论讨论对当前具体问题有用。

你可以通过逐个矩阵来处理这个问题，使用 scipy 库，它提供了 pinv 方法（链接），可以计算摩尔-彭若斯伪逆。下面是一个例子：

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)

不过，这个方法并没有利用矩阵是对称的这个特点。你也可以尝试使用 numpy 提供的 pinv 版本，据说它更快，而且有所不同。可以查看 这篇文章。

Answer 2

NumPy 1.8版本加入了线性代数的gufuncs，这些功能正好满足你的需求。虽然np.linalg.pinv并不是gufunc，但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

这其实还不错，虽然它的速度明显比直接调用np.linalg.inv慢了4倍，但后者无法正确处理奇异数组：

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