使用pytorch验证卷积定理
Verify convolution theorem using pytorch
这个定理基本上可以表述如下：
F(f*g) = F(f)xF(g)
我知道这个定理，但我就是无法使用 pytorch 重现结果。
下面是一个可重现的代码：
import torch
import torch.nn.functional as F
# calculate f*g
f = torch.ones((1,1,5,5))
g = torch.tensor(list(range(9))).view(1,1,3,3).float()
conv = F.conv2d(f, g, bias=None, padding=2)
# calculate F(f*g)
F_fg = torch.rfft(conv, signal_ndim=2, onesided=False)
# calculate F x G
f = f.squeeze()
g = g.squeeze()
# need to pad into at least [w1+w2-1, h1+h2-1], which is 7 in our case.
size = f.size(0) + g.size(0) - 1
f_new = torch.zeros((7,7))
g_new = torch.zeros((7,7))
f_new[1:6,1:6] = f
g_new[2:5,2:5] = g
F_f = torch.rfft(f_new, signal_ndim=2, onesided=False)
F_g = torch.rfft(g_new, signal_ndim=2, onesided=False)
FxG = torch.mul(F_f, F_g)
print(FxG - F_fg)
这是 print(FxG - F_fg)
的结果
tensor([[[[[ 0.0000e+00, 0.0000e+00],
[ 4.1426e+02, 1.7270e+02],
[-3.6546e+01, 4.7600e+01],
[-1.0216e+01, -4.1198e+01],
[-1.0216e+01, -2.0223e+00],
[-3.6546e+01, -6.2804e+01],
[ 4.1426e+02, -1.1427e+02]],
...
[[ 4.1063e+02, -2.2347e+02],
[-7.6294e-06, 2.2817e+01],
[-1.9024e+01, -9.0105e+00],
[ 7.1708e+00, -4.1027e+00],
[-2.6739e+00, -1.1121e+01],
[ 8.8471e+00, 7.1710e+00],
[ 4.2528e+01, 9.7559e+01]]]]])
你可以看到差异并不总是 0。
谁能告诉我为什么以及如何正确执行此操作？
谢谢
所以我仔细研究了您到目前为止所做的事情。我在您的代码中确定了三个错误来源。我会尽力在这里充分解决每个问题。
1。复数运算
PyTorch 目前不支持复数乘法 (AFAIK)。 FFT 操作简单地 returns 具有实部和虚部的张量。我们需要显式地编写复数乘法代码，而不是使用 torch.mul 或 * 运算符。
(a + ib) * (c + id) = (a*c - b*d) + i(a*d + b*c)
2。卷积的定义
CNN文献中经常使用的“卷积”的定义，其实和讨论卷积定理时使用的定义是不一样的。我不会详细介绍，但是 theoretical definition 在滑动和相乘之前翻转内核。相反，pytorch、tensorflow、caffe 等中的卷积运算...不会进行这种翻转。
为了解决这个问题，我们可以在应用 FFT 之前简单地翻转 g（水平和垂直）。
3。锚点位置
假定使用卷积定理时的锚点是填充后的左上角g。同样，我不会对此进行详细说明，但这就是数学运算的结果。
第二点和第三点用例子可能更容易理解。假设您使用了以下 g
[1 2 3]
[4 5 6]
[7 8 9]
而不是 g_new 成为
[0 0 0 0 0 0 0]
[0 0 0 0 0 0 0]
[0 0 1 2 3 0 0]
[0 0 4 5 6 0 0]
[0 0 7 8 9 0 0]
[0 0 0 0 0 0 0]
[0 0 0 0 0 0 0]
实际上应该是
[5 4 0 0 0 0 6]
[2 1 0 0 0 0 3]
[0 0 0 0 0 0 0]
[0 0 0 0 0 0 0]
[0 0 0 0 0 0 0]
[0 0 0 0 0 0 0]
[8 7 0 0 0 0 9]
在这里我们垂直和水平翻转内核，然后应用循环移位，使内核的中心位于左上角。
我最终重写了你的大部分代码并对其进行了一些概括。最复杂的操作是正确定义 g_new。我决定使用 meshgrid 和模运算来同时翻转和移动索引。如果这里的某些内容对您没有意义，请发表评论，我会尽力澄清。
import torch
import torch.nn.functional as F
def conv2d_pyt(f, g):
assert len(f.size()) == 2
assert len(g.size()) == 2
f_new = f.unsqueeze(0).unsqueeze(0)
g_new = g.unsqueeze(0).unsqueeze(0)
pad_y = (g.size(0) - 1) // 2
pad_x = (g.size(1) - 1) // 2
fcg = F.conv2d(f_new, g_new, bias=None, padding=(pad_y, pad_x))
return fcg[0, 0, :, :]
def conv2d_fft(f, g):
assert len(f.size()) == 2
assert len(g.size()) == 2
# in general not necessary that inputs are odd shaped but makes life easier
assert f.size(0) % 2 == 1
assert f.size(1) % 2 == 1
assert g.size(0) % 2 == 1
assert g.size(1) % 2 == 1
size_y = f.size(0) + g.size(0) - 1
size_x = f.size(1) + g.size(1) - 1
f_new = torch.zeros((size_y, size_x))
g_new = torch.zeros((size_y, size_x))
# copy f to center
f_pad_y = (f_new.size(0) - f.size(0)) // 2
f_pad_x = (f_new.size(1) - f.size(1)) // 2
f_new[f_pad_y:-f_pad_y, f_pad_x:-f_pad_x] = f
# anchor of g is 0,0 (flip g and wrap circular)
g_center_y = g.size(0) // 2
g_center_x = g.size(1) // 2
g_y, g_x = torch.meshgrid(torch.arange(g.size(0)), torch.arange(g.size(1)))
g_new_y = (g_y.flip(0) - g_center_y) % g_new.size(0)
g_new_x = (g_x.flip(1) - g_center_x) % g_new.size(1)
g_new[g_new_y, g_new_x] = g[g_y, g_x]
# take fft of both f and g
F_f = torch.rfft(f_new, signal_ndim=2, onesided=False)
F_g = torch.rfft(g_new, signal_ndim=2, onesided=False)
# complex multiply
FxG_real = F_f[:, :, 0] * F_g[:, :, 0] - F_f[:, :, 1] * F_g[:, :, 1]
FxG_imag = F_f[:, :, 0] * F_g[:, :, 1] + F_f[:, :, 1] * F_g[:, :, 0]
FxG = torch.stack([FxG_real, FxG_imag], dim=2)
# inverse fft
fcg = torch.irfft(FxG, signal_ndim=2, onesided=False)
# crop center before returning
return fcg[f_pad_y:-f_pad_y, f_pad_x:-f_pad_x]
# calculate f*g
f = torch.randn(11, 7)
g = torch.randn(5, 3)
fcg_pyt = conv2d_pyt(f, g)
fcg_fft = conv2d_fft(f, g)
avg_diff = torch.mean(torch.abs(fcg_pyt - fcg_fft)).item()
print('Average difference:', avg_diff)
这给了我
Average difference: 4.6866085767760524e-07
这非常接近于零。我们没有得到完全零的原因仅仅是由于浮点错误。
这个定理基本上可以表述如下：
F(f*g) = F(f)xF(g)
我知道这个定理，但我就是无法使用 pytorch 重现结果。
下面是一个可重现的代码：
import torch
import torch.nn.functional as F
# calculate f*g
f = torch.ones((1,1,5,5))
g = torch.tensor(list(range(9))).view(1,1,3,3).float()
conv = F.conv2d(f, g, bias=None, padding=2)
# calculate F(f*g)
F_fg = torch.rfft(conv, signal_ndim=2, onesided=False)
# calculate F x G
f = f.squeeze()
g = g.squeeze()
# need to pad into at least [w1+w2-1, h1+h2-1], which is 7 in our case.
size = f.size(0) + g.size(0) - 1
f_new = torch.zeros((7,7))
g_new = torch.zeros((7,7))
f_new[1:6,1:6] = f
g_new[2:5,2:5] = g
F_f = torch.rfft(f_new, signal_ndim=2, onesided=False)
F_g = torch.rfft(g_new, signal_ndim=2, onesided=False)
FxG = torch.mul(F_f, F_g)
print(FxG - F_fg)
这是 print(FxG - F_fg)
的结果tensor([[[[[ 0.0000e+00, 0.0000e+00],
[ 4.1426e+02, 1.7270e+02],
[-3.6546e+01, 4.7600e+01],
[-1.0216e+01, -4.1198e+01],
[-1.0216e+01, -2.0223e+00],
[-3.6546e+01, -6.2804e+01],
[ 4.1426e+02, -1.1427e+02]],
...
[[ 4.1063e+02, -2.2347e+02],
[-7.6294e-06, 2.2817e+01],
[-1.9024e+01, -9.0105e+00],
[ 7.1708e+00, -4.1027e+00],
[-2.6739e+00, -1.1121e+01],
[ 8.8471e+00, 7.1710e+00],
[ 4.2528e+01, 9.7559e+01]]]]])
你可以看到差异并不总是 0。
谁能告诉我为什么以及如何正确执行此操作？
谢谢
所以我仔细研究了您到目前为止所做的事情。我在您的代码中确定了三个错误来源。我会尽力在这里充分解决每个问题。
1。复数运算
PyTorch 目前不支持复数乘法 (AFAIK)。 FFT 操作简单地 returns 具有实部和虚部的张量。我们需要显式地编写复数乘法代码，而不是使用 torch.mul 或 * 运算符。
(a + ib) * (c + id) = (a*c - b*d) + i(a*d + b*c)
2。卷积的定义
CNN文献中经常使用的“卷积”的定义，其实和讨论卷积定理时使用的定义是不一样的。我不会详细介绍，但是 theoretical definition 在滑动和相乘之前翻转内核。相反，pytorch、tensorflow、caffe 等中的卷积运算...不会进行这种翻转。
为了解决这个问题，我们可以在应用 FFT 之前简单地翻转 g（水平和垂直）。
3。锚点位置
假定使用卷积定理时的锚点是填充后的左上角g。同样，我不会对此进行详细说明，但这就是数学运算的结果。
第二点和第三点用例子可能更容易理解。假设您使用了以下 g
[1 2 3]
[4 5 6]
[7 8 9]
而不是 g_new 成为
[0 0 0 0 0 0 0]
[0 0 0 0 0 0 0]
[0 0 1 2 3 0 0]
[0 0 4 5 6 0 0]
[0 0 7 8 9 0 0]
[0 0 0 0 0 0 0]
[0 0 0 0 0 0 0]
实际上应该是
[5 4 0 0 0 0 6]
[2 1 0 0 0 0 3]
[0 0 0 0 0 0 0]
[0 0 0 0 0 0 0]
[0 0 0 0 0 0 0]
[0 0 0 0 0 0 0]
[8 7 0 0 0 0 9]
在这里我们垂直和水平翻转内核，然后应用循环移位，使内核的中心位于左上角。
我最终重写了你的大部分代码并对其进行了一些概括。最复杂的操作是正确定义 g_new。我决定使用 meshgrid 和模运算来同时翻转和移动索引。如果这里的某些内容对您没有意义，请发表评论，我会尽力澄清。
import torch
import torch.nn.functional as F
def conv2d_pyt(f, g):
assert len(f.size()) == 2
assert len(g.size()) == 2
f_new = f.unsqueeze(0).unsqueeze(0)
g_new = g.unsqueeze(0).unsqueeze(0)
pad_y = (g.size(0) - 1) // 2
pad_x = (g.size(1) - 1) // 2
fcg = F.conv2d(f_new, g_new, bias=None, padding=(pad_y, pad_x))
return fcg[0, 0, :, :]
def conv2d_fft(f, g):
assert len(f.size()) == 2
assert len(g.size()) == 2
# in general not necessary that inputs are odd shaped but makes life easier
assert f.size(0) % 2 == 1
assert f.size(1) % 2 == 1
assert g.size(0) % 2 == 1
assert g.size(1) % 2 == 1
size_y = f.size(0) + g.size(0) - 1
size_x = f.size(1) + g.size(1) - 1
f_new = torch.zeros((size_y, size_x))
g_new = torch.zeros((size_y, size_x))
# copy f to center
f_pad_y = (f_new.size(0) - f.size(0)) // 2
f_pad_x = (f_new.size(1) - f.size(1)) // 2
f_new[f_pad_y:-f_pad_y, f_pad_x:-f_pad_x] = f
# anchor of g is 0,0 (flip g and wrap circular)
g_center_y = g.size(0) // 2
g_center_x = g.size(1) // 2
g_y, g_x = torch.meshgrid(torch.arange(g.size(0)), torch.arange(g.size(1)))
g_new_y = (g_y.flip(0) - g_center_y) % g_new.size(0)
g_new_x = (g_x.flip(1) - g_center_x) % g_new.size(1)
g_new[g_new_y, g_new_x] = g[g_y, g_x]
# take fft of both f and g
F_f = torch.rfft(f_new, signal_ndim=2, onesided=False)
F_g = torch.rfft(g_new, signal_ndim=2, onesided=False)
# complex multiply
FxG_real = F_f[:, :, 0] * F_g[:, :, 0] - F_f[:, :, 1] * F_g[:, :, 1]
FxG_imag = F_f[:, :, 0] * F_g[:, :, 1] + F_f[:, :, 1] * F_g[:, :, 0]
FxG = torch.stack([FxG_real, FxG_imag], dim=2)
# inverse fft
fcg = torch.irfft(FxG, signal_ndim=2, onesided=False)
# crop center before returning
return fcg[f_pad_y:-f_pad_y, f_pad_x:-f_pad_x]
# calculate f*g
f = torch.randn(11, 7)
g = torch.randn(5, 3)
fcg_pyt = conv2d_pyt(f, g)
fcg_fft = conv2d_fft(f, g)
avg_diff = torch.mean(torch.abs(fcg_pyt - fcg_fft)).item()
print('Average difference:', avg_diff)
这给了我
Average difference: 4.6866085767760524e-07
这非常接近于零。我们没有得到完全零的原因仅仅是由于浮点错误。