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

Question

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

D为我的图像中的每个体素存储3x3对称矩阵的6个独特元素.我需要在矢量化代码中同时找到所有行*cols*deps矩阵的Moore-Penrose伪逆(循环遍历每个图像体素,并且在Python中反转太慢).

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

然而,对于那些单数的矩阵(并且肯定会有一些),我需要Moore-Penrose伪逆.我可以得出一个真实的,奇异的,对称的3x3矩阵的MP的解析公式,但它是一个非常讨厌/冗长的公式,因此会涉及非常大量(元素方面)矩阵算法和相当多的混乱寻找代码.

因此,我想知道是否有一种方法可以同时在数字上同时找到所有这些矩阵的MP伪逆.有没有办法做到这一点？

感激不尽,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

这真的不是那么糟糕,虽然它比简单调用明显(4倍)慢np.linalg.inv,但这当然无法正确处理奇异数组:

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