极大的加权平均

Posted 2023-03-11

【中文标题】极大的加权平均

【英文标题】：Extremely large weighted average

【发布时间】：2011-12-14 14:09:43

【问题描述】：

我正在使用 64 位 matlab 和 32g 的 RAM（只是让你知道）。

我有一个包含 130 万个数字（整数）的文件（向量）。我想制作另一个相同长度的向量，其中每个点都是整个第一个向量的加权平均值，由与该位置的反距离加权（实际上是位置 ^-0.1，而不是 ^-1，但用于示例目的） .我不能使用 matlab 的“过滤器”功能，因为它只能平均当前点之前的东西，对吧？为了更清楚地解释，这里以 3 个元素为例

data = [ 2 6 9 ]
weights = [ 1 1/2 1/3; 1/2 1 1/2; 1/3 1/2 1 ]
results=data*weights= [ 8 11.5 12.666 ]
i.e.
8 = 2*1 + 6*1/2 + 9*1/3
11.5 = 2*1/2 + 6*1 + 9*1/2
12.666 = 2*1/3 + 6*1/2 + 9*1

所以新向量中的每个点都是整个第一个向量的加权平均值，权重为 1/（与该位置的距离+1）。

我可以为每个点重新制作权重向量，然后逐个元素计算结果向量，但这需要 for 循环的 130 万次迭代，每个迭代包含 130 万次乘法。我宁愿使用直接矩阵乘法，将 1x1.3mil 乘以 1.3milx1.3mil，理论上可行，但我无法加载那么大的矩阵。

然后我尝试使用 shell 脚本制作矩阵并在 matlab 中对其进行索引，因此一次只调用矩阵的相关列，但这也需要很长时间。

我不必在 matlab 中执行此操作，因此人们对使用如此大的数字并获得平均值提出任何建议将不胜感激。由于我使用的是 ^-0.1 而不是 ^-1 的权重，因此它不会下降得那么快 - 与原始点的权重 1 相比，第 100 万个点的权重仍然为 0.25，所以我不能只是削减它当它变大时关闭。

希望这已经足够清楚了吗？

这是下面答案的代码（所以可以格式化？）：

data = load('/Users/mmanary/Documents/test/insertion.txt');
data=data.';
total=length(data);
x=1:total;
datapad=[zeros(1,total) data];
weights = ([(total+1):-1:2 1:total]).^(-.4);
weights = weights/sum(weights);
Fdata = fft(datapad);
Fweights = fft(weights);
Fresults = Fdata .* Fweights;
results = ifft(Fresults);
results = results(1:total);
plot(x,results)

【问题讨论】：

对于不熟悉 FFT 方法来解决此类问题的任何人（请参阅下面的答案），我强烈推荐免费在线 Scientist and Engineer's Guide to Digital Signal Processing。即使对于像我这样的非 DSP 人来说，了解基础知识也是无价的。使用 FFT 与蛮力卷积时的加速几乎就像魔术一样。 :)

【参考方案1】：

做到这一点的唯一明智方法是使用 FFT 卷积，作为 filter 函数和类似函数的基础。手动操作很容易：

% Simulate some data
n = 10^6;
x = randi(10,1,n);
xpad = [zeros(1,n) x];

% Setup smoothing kernel
k = 1 ./ [(n+1):-1:2 1:n];

% FFT convolution
Fx = fft(xpad);
Fk = fft(k);

Fxk = Fx .* Fk;

xk = ifft(Fxk);
xk = xk(1:n);

n=10^6 需要 不到半秒！

【讨论】：

而且对于太长的信号，FFT 可能会被分解成帧。

@MicahManary 太棒了！祝你的项目好运。

@JohnColby 所以卷积工作，但我如何除以总权重（以获得实际平均值）。边缘处的点需要除以比中间点小得多的数字。我的代码现在有问题（因此可以格式化）。

我相信只需扩展您的内核（即weights），使其总和为 1。

@JohnColby 当然 - 谢谢。上面编辑的答案代码可以工作！

【参考方案2】：

我不能使用 matlab 的 'filter' 函数，因为它只能取平均值 当前点之前的事情，对吧？

这是不正确的。您始终可以从数据或过滤后的数据中添加样本（即添加或删除零）。由于使用filter 过滤（顺便说一下，您也可以使用conv）是一个线性动作，它不会改变结果（就像添加和删除零，什么都不做，然后过滤。然后线性允许你交换顺序以添加样本 -> 过滤 -> 删除样本）。

无论如何，在您的示例中，您可以将平均内核设为：

weights = 1 ./ [3 2 1 2 3]; % this kernel introduces a delay of 2 samples

然后简单地说：

result =  filter(w,1,[data, zeros(1,3)]); % or conv (data, w)
% removing the delay introduced by the kernel
result = result (3:end-1);

【参考方案3】：

蛮力方法可能对您有用，但需要进行一些小的优化。

创建权重的 ^-0.1 操作将比计算加权均值的 + 和 * 操作花费更长的时间，但是您可以在所有百万加权均值操作中重复使用权重。算法变为：

创建一个权重向量，其中包含任何计算所需的所有权重： weights = (-n:n).^-0.1

对于向量中的每个元素：

索引weights 向量的相关部分以将当前元素视为“中心”。

对权重部分和整个向量执行加权平均。这可以通过快速矢量点乘和标量除法来完成。

主循环n^2 加减。 n 等于 130 万，即 3.4 万亿次操作。现代 3GHz CPU 的单核每秒可以进行 60 亿次加法/乘法运算，因此大约需要 10 分钟。再加上索引weights 向量和开销的时间，我仍然估计您可以在半小时内完成。

【参考方案4】：

您只考虑了 2 个选项： 将 1.3M*1.3M 矩阵与向量相乘一次或将 2 个 1.3M 向量相乘 1.3M 次。

但是您可以将权重矩阵划分为任意数量的子矩阵，然后将 n*1.3M 矩阵与向量乘以 1.3M/n 次。

我假设当迭代次数最少并且 n 会创建适合您的内存的最大子矩阵时，最快将是这样，而不会使您的计算机开始将页面交换到您的硬盘驱动器。

根据您的内存大小，您应该从 n=5000 开始。

您还可以使用 parfor 使其更快（n 除以处理器数量）。

【参考方案5】：

这可能不是最好的方法，但如果有大量内存，您绝对可以并行化该过程。

您可以构建由原始矩阵的条目组成的稀疏矩阵，其值为 i^(-1)（其中 i = 1 .. 1.3 million），将它们与原始向量相乘，然后将所有结果相加。

因此，对于您的示例，产品本质上是：

a = rand(3,1);
b1 = [1 0 0;
      0 1 0;
      0 0 1];
b2 = [0 1 0;
      1 0 1;
      0 1 0] / 2;
b3 = [0 0 1;
      0 0 0;
      1 0 0] / 3;

c = sparse(b1) * a + sparse(b2) * a + sparse(b3) * a;

当然，您不会以这种方式构造稀疏矩阵。如果您想减少内部循环的迭代次数，则可以在每个矩阵中拥有多个 i。

查看 MATLAB 中的 parfor 循环：http://www.mathworks.com/help/toolbox/distcomp/parfor.html