将二维高斯核分离成两个一维核

separate 2D gaussian kernel into two 1D kernels

计算了一个高斯核,并通过查看核的秩来检查它是否可以分离。

kernel = gaussian_kernel(kernel_size,sigma)
print(kernel)

[[ 0.01054991  0.02267864  0.0292689   0.02267864  0.01054991]
 [ 0.02267864  0.04875119  0.06291796  0.04875119  0.02267864]
 [ 0.0292689   0.06291796  0.0812015   0.06291796  0.0292689 ]
 [ 0.02267864  0.04875119  0.06291796  0.04875119  0.02267864]
 [ 0.01054991  0.02267864  0.0292689   0.02267864  0.01054991]]

 rank = np.linalg.matrix_rank(kernel)

if rank == 1:
    print('The Kernel is separable')
else:
    print('The kernel is not separable')

现在我认为分离是不正确的。我按以下方式进行:

 u,s,v = np.linalg.svd(kernel)
 k1 = (u[:,0] * np.sqrt(s[0]))[np.newaxis].T
 k2 = v[:,0] * np.sqrt(s[0]) 

然后我将上面的两个内核相乘,得到了原来的内核。但是我没看懂。

if not  np.all(k1 * k2 == kernel):
    print('k1 * k2 is not equal to kernel')

我假设我尝试使用 svd 和进一步进行的分离是不正确的。一些解释会有所帮助。

矩阵等级 1 意味着所有行要么为零,要么在缩放之前相同,列也是如此。它们也达到等于这两个因子的比例。 因此,您可以使用

之类的方法恢复它们
I,J = np.unravel_index(np.abs(kernel).argmax(), kernel.shape)
f1 = np.nansum(kernel / (kernel[None,:,J]@kernel),1,keepdims=True)
f2 = np.nansum(kernel / (kernel@kernel[I,:,None]),0,keepdims=True)
scaling = np.sqrt(np.abs(kernel).sum()/np.abs(f1*f2).sum())
f1 *= scaling * np.sign(f1[I,0]) * np.sign(kernel[I,J])
f2 *= scaling * np.sign(f2[0,J])

请注意,大部分的复杂性来自于我试图对尽可能多的数据进行平均。一个更简单但我认为在数值上不太稳定的方法是

I,J = np.unravel_index(np.abs(kernel).argmax(), kernel.shape)
f1 = kernel[:,J,None]
f2 = kernel[None,I,:] / kernel[I,J]

当然,一旦你的索引正确,你的方法也可以工作:

k1 = u[:,0,None] * np.sqrt(s[0])
k2 = v[None,0,:] * np.sqrt(s[0])

np.allclose(kernel, k1*k2)
# True