计算离散卷积“乘积”的有效方法

Posted 2023-02-14

【中文标题】计算离散卷积“乘积”的有效方法

【英文标题】：Efficient way to calculate the "product" of a discrete convolution

【发布时间】：2017-05-13 00:22:15

【问题描述】：

我正在寻找一种优雅的方法来计算离散卷积的“乘积”，而不是求和。

这是一个离散卷积的公式：

在这种情况下我们可以使用：conv(x,y)

现在我想实现这些操作

当然我可以使用循环，但我正在寻找一个技巧来线性化这个操作。

示例：

f = [2 4 3 9 7 1]
g = [3 2 1]

dist = length(g)-1;

for ii = 1:length(f)-dist
    x(ii) = prod(f(ii:ii+dist).*g)
end

x =

144    648   1134    378

【问题讨论】：

如果您的输入 (f,g) 始终包含整数，则无论算法如何，您都可以通过将它们转换为 uint16 等来获得一些速度（假设您可以排除溢出）。

【参考方案1】：

cumprod解决方案：（非常高效）

>> pf = cumprod(f);
>> x = prod(g).*pf(numel(g):end)./[1 pf(1:(end-numel(g)))]

x =

         144         648        1134         378

这首先使用cumprod 获取f 的累积乘积。通过将每个元素除以它之前的 3 个元素的累积乘积，我们得到每个 numel(g)-wide 滑动窗口沿 f 的乘积。然后乘以g的元素的乘积。

注意：当f 有许多元素或极值（大或小）时，您在执行累积乘积时可能会遇到准确性问题或下溢/上溢。缓解这种情况的一种潜在方法是在累积乘积之前对f 应用缩放，然后再撤消它：

c = ...set a scaling factor...
pf = cumprod(f./c);
x = prod(c.*g).*pf(numel(g):end)./[1 pf(1:(end-numel(g)))];

c 的选择可能类似于 mean(abs(f)) 或 max(abs(f))，以便缩放后的 f 提供更好的有界结果（即更接近 1 的值）。这不会明显改变下面的计时结果。

---

hankel 解决方案：（效率不高但仍然很有趣）

>> x = prod(g).*prod(hankel(f(1:numel(g)), f(numel(g):end)))

x =

         144         648        1134         378

对hankel 的调用会创建一个矩阵，其中每一列都包含numel(g) 宽滑动窗口之一的内容。将每列的乘积乘以g 元素的乘积，即可得出答案。但是，对于大型向量 f 和/或 g，这可能涉及大量额外计算并占用大量内存。

---

计时结果：

我测试了 6 个解决方案（您问题中的循环，rahnema 的 2 个解决方案（conv(log) 和 movsum(log)），Luis 的 bsxfun 解决方案，以及我的 cumprod 和 hankel 解决方案) 使用 f = rand(1,1000000); 和 g = rand(1,100); 并平均超过 40 次迭代。这是我得到的（运行 Windows 7 x64、16 GB RAM、MATLAB R2016b）：

solution    | avg. time (s)
------------+---------------
loop        |       1.10671
conv(log)   |       0.04984
movsum(log) |       0.03736
bsxfun      |       1.20472
cumprod     |       0.01469
hankel      |       1.17704

【讨论】：

我认为如果f 包含诸如4 5 10 而不是rand 之类的值，cumprod(f) 的结果将不会那么可靠

@rahnema1：你是对的。我谈到了如何补偿。

感谢您的回答。 cumprod 方法确实很高效，但是如果 f 包含 0，这个方法就会失败，而且 IMO 实现这样的解决方案有点危险。

【参考方案2】：

另一个解决方案部分受Dev-iL answer 启发，解决了相对相同的问题

exp(sum(log(g))+conv(log(f),[1 1 1],'valid'))

或

exp(sum(log(g))+movsum(log(f),numel(g),'Endpoints', 'discard'))

自从exp(sum(log(x))) = prod(x)

但这里不是一个向量，而是两个向量 f 和 g。

所需的公式可以重新表述为：

八度音阶：

f= rand(1,1000000);
g= rand(1,100);

disp('----------EXP(LOG)------')
tic
    exp(sum(log(g))+conv(log(f),ones(1,numel(g))));
toc

disp('----------BSXFUN------')
tic
    ind = bsxfun(@plus, 0:numel(f)-numel(g), (1:numel(g)).');
    x = prod(bsxfun(@times, f(ind), g(:)),1);
toc
disp('----------HANKEL------')
tic
    prod(g)*prod(hankel(f(1:numel(g)), f(numel(g):end)));
toc

disp('----------CUMPROD-----')
tic
    pf = cumprod(f);
    x = prod(g)*pf(numel(g):end)./[1 pf(1:(end-numel(g)))];
toc

结果：

----------EXP(LOG)------%rahnema1
Elapsed time is 0.211445 seconds.
----------BSXFUN--------%Luis Mendo
Elapsed time is 1.94182 seconds.
----------HANKEL--------%gnovice
Elapsed time is 1.46593 seconds.
----------CUMPROD-------%gnovice
Elapsed time is 0.00748992 seconds.

【讨论】：

要匹配问题中示例的结果，您需要将..., 'valid') 添加到conv 调用，并将..., 'Endpoints', 'discard'); 添加到movsum 调用。

【参考方案3】：

这是一种避免循环的方法：

ind = bsxfun(@plus, 0:numel(f)-numel(g), (1:numel(g)).');
x = prod(bsxfun(@times, f(ind), g(:)), 1);

它的工作原理如下。 ind 表示索引的滑动窗口，每一列对应一个位移。例如，如果g 的大小为3，则矩阵ind 将是

1 2 3 4 ...
2 3 4 5 ...
3 4 5 6 ...

这用于索引f。结果乘以g（通过广播）作为一列。最后计算每列元素的乘积。