将二维高斯核分离成两个一维核
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
计算了一个高斯核,并通过查看核的秩来检查它是否可以分离。
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