11

我正在使用 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)
4

5 回答 5

11

唯一明智的方法是使用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！

于 2011-10-25T07:40:47.503 回答
3

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

您可以构建由原始矩阵的条目组成的稀疏矩阵，这些条目具有值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's。

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

于 2011-10-25T05:37:46.340 回答
2

我不能使用 matlab 的“过滤器”功能，因为它只能平均当前点之前的东西，对吧？

这是不正确的。您始终可以从您的数据或过滤的数据中添加样本（即添加或删除零）。由于过滤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);
于 2011-10-25T07:35:35.173 回答
0

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

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

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

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

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

于 2011-10-25T07:35:51.897 回答
0

蛮力的方式可能会为你工作，在混合中进行一个小的优化。

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

  • 使用任何计算所需的所有权重创建一个权重向量： weights = (-n:n).^-0.1

  • 对于向量中的每个元素：

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

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

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

于 2011-10-25T07:43:18.460 回答