Python，许多 3x3、奇异、对称、矩阵的同时伪反转

Question

我有一个尺寸为行 x 列 x 深度的 3D 图像。对于图像中的每个体素，我计算了一个 3x3 实对称矩阵。它们存储在数组 D 中，因此具有形状 (rows, cols, deps, 6)。

D 为图像中的每个体素存储 3x3 对称矩阵的 6 个唯一元素。我需要同时/在矢量化代码中找到所有 row*cols*deps 矩阵的 Moore-Penrose 伪逆(遍历每个图像体素并在 Python 中求逆太慢)。

其中一些 3x3 对称矩阵是非奇异矩阵，我可以在向量化代码中找到它们的逆矩阵，使用非奇异 3x3 对称矩阵的真逆的解析公式，我已经做到了。

但是，对于那些奇异的矩阵(并且肯定有一些)，我需要 Moore-Penrose 伪逆。我可以推导出一个实数、奇异、对称的 3x3 矩阵的 MP 的分析公式，但它是一个非常讨厌/冗长的公式，因此会涉及非常多的(元素方面的)矩阵算术和相当多的困惑寻找代码。

因此，我想知道是否有一种方法可以在数值上一次同时找到所有这些矩阵的 MP 伪逆。有办法做到这一点吗？

感激，广发

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

这确实没那么糟糕，虽然它明显比简单地调用 np.linalg.inv 慢 (4x)，但这当然无法正确处理奇异数组:

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