计算离散卷积的 "product" 的有效方法

Question

我正在寻找一种优雅的方法来计算离散卷积的 "product"，而不是总和。

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

在这种情况下我们可以使用：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

Answer 1

这里有一个避免循环的方法：

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

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

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

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

Answer 2

`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) 宽滑动 window 沿 f 的乘积。然后乘以 g.

的元素的乘积

注意： 当 f 有很多元素或极值（大或小）时，您可能运行会遇到准确性问题或 underflow/overflow 执行累积乘积时。一种可能的缓解方法是在累积乘积之前对 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)-wide 滑动 windows 的内容。将每列的乘积乘以 g 的元素的乘积即可得到答案。但是，对于大向量 f and/or g 这可能涉及大量额外计算并使用大量内存。

对结果进行计时：

我使用 f = rand(1,1000000); 和 [=38 测试了 6 个解决方案（你问题中的循环， (conv(log) and movsum(log)), the bsxfun solution from 中的 2 个解决方案，以及我的 cumprod 和 hankel 解决方案） =] 并平均超过 40 次迭代。这是我得到的（运行ning 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

Answer 3

另一个解决方案，部分灵感来自相对相同的问题

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.

计算离散卷积的 "product" 的有效方法

Efficient way to calculate the "product" of a discrete convolution

arrays

matlab

multiplication

convolution

sliding-window

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

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

对结果进行计时：

计算离散卷积的 "product" 的有效方法

Efficient way to calculate the "product" of a discrete convolution

arrays

matlab

multiplication

convolution

sliding-window

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

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

对结果进行计时：

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

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