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

Question

我有一个三维图像，其尺寸为行×列×深度。对于图像中的每个体素，我计算了一个3x3实对称矩阵。它们存储在数组D中，因此其形状为（行，列，深度，6）。

D存储了图像中每个体素的3x3对称矩阵的6个唯一元素。我需要同时/向量化地找到所有行*列*深度矩阵的Moore-Penrose伪逆（在Python中循环遍历每个图像体素并求逆太慢了）。

其中一些3x3对称矩阵是非奇异的，我可以使用真实逆矩阵的解析公式找到它们的逆矩阵，并且我已经做到了这一点。

但是，对于那些奇异的矩阵（肯定会有一些），我需要Moore-Penrose伪逆。我可以推导出实数奇异对称3x3矩阵的MP伪逆的解析公式，但它是一个非常麻烦/冗长的公式，因此将涉及非常多的（逐元素）矩阵运算和相当多看起来令人困惑的代码。

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

感激不尽， GF

- NLi10Me

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

这其实也不算太糟糕，虽然速度明显慢了4倍，但是比起直接调用np.linalg.inv要更可靠，因为后者在处理奇异数组时会发生错误：

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

- Jaime

---

非常酷，这是新版本的一个惊人特性！带有新广义ufunc列表的ipython笔记本在此处。 - gg349

有趣的是，这个笔记本来自我在圣地亚哥Python用户组的演讲！我从未想过除了那些参加演讲的人之外，会有其他人发现它。互联网真是一个小地方... - Jaime

谢谢你的回答，Jamie。这看起来是一些强大的代码...但我还没有完全理解它。我对Python相对较新。 gufunc是GPU实现吗？我相信我可以查找这个问题，但想问一下是否有其他人不知道。当然，在422毫秒内取60,000个任意3x3实/对称矩阵的伪逆是非常棒的。 - NLi10Me

与GPU无关。您可能需要阅读有关通用函数（即ufunc）的文档，此处。然后，您可能希望阅读有关广义通用函数（即gufunc）的内容，此处。基本上，ufunc将在标量上定义的操作（如加法），并对其输入数组的所有元素逐个运行它。Gufunc使用在数组上定义的操作执行类似的操作，因此单个调用会在许多堆叠的数组上运行相同的操作。 - Jaime

Answer 2

编辑：请查看@Jaime的回答。现在只有对特定问题的讨论在本答案的评论中是有用的。

您可以使用scipy逐个矩阵地进行操作，它提供了pinv(link)来计算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)

虽然矩阵是对称的，但这并没有被利用。您可能还想尝试numpy提供的pinv版本，据说更快且不同。请参见this帖子。

- gg349

---

我知道scipy.linalg.pinv，但我的问题是如何避免逐个矩阵进行操作。对所有行列依赖矩阵进行循环太慢了，不适用于我的应用程序。 - NLi10Me

抱歉，我有这样的印象，因为它们中只有一些是单数。你对代码进行了一些分析吗？你确定调用 pinv 是你的问题吗？通常有多少矩阵是奇异的？我怀疑很少（？）。如果是这种情况，试图向量化 pinv 就是过早的优化。 - gg349

我认为你是对的。奇异矩阵只占了一小部分。我会尝试将它们隔离出来，仅循环处理这些矩阵，并使用pinv。也许这样就足够快了。最初，我在所有体素上都有很多计算（不仅仅是调用pinv），所以我不认为是pinv特别拖慢了速度。大多数其他操作已经向量化，这是唯一的麻烦制造者。 - NLi10Me

好的，请随时更新。您也可以尝试使用 numpy 版本，可能会更快或更慢，请参见上面的链接。 - gg349

我确实使用了numpy版本。由于奇异矩阵的百分比非常小，它足够快。非常感谢您的帮助。然而，原始问题仍然存在，有没有一种方法可以在Python中同时反转一组小矩阵？ - NLi10Me

我想另一个“非Python”解决方案是并行处理，但我正在寻求Python内部的解决方案。 - NLi10Me