我使用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)。
我可以重新制作每个点的权重向量,然后逐个元素地计算结果向量,但这需要130万次for循环迭代,每个迭代包含130万次乘法。我宁愿使用直接矩阵乘法,将1x1.3mil乘以1.3milx1.3mil,这在理论上是有效的,但我无法加载大的矩阵。
然后我尝试使用shell脚本制作矩阵并在matlab中对其进行索引,因此一次只调用矩阵的相关列,但这也需要很长时间。
我没有必要在matlab中这样做,所以人们对利用如此大的数字和获得平均值的任何建议都将不胜感激。因为我使用的是^ -0.1的重量,而不是^ -1,所以它不会快速下降 - 百万分之一点仍然加权为0.25,而原始点加权为1,所以我不能只是削减它因为它变大了。
希望这很清楚吗?
以下是答案的代码(因此可以格式化?):
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卷积,作为支撑 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
!
这可能不是最好的方法,但是有了大量的内存,你肯定可以并行化这个过程。
您可以构造由原始矩阵的条目组成的稀疏矩阵,这些条目具有值 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
在每个矩阵中。
看看 parfor
在MATLAB中循环: http://www.mathworks.com/help/toolbox/distcomp/parfor.html
我不能使用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);
您只考虑了两个选项:
将1.3M * 1.3M矩阵与载体相乘一次或将2个1.3M载体乘以1.3M倍。
但是你可以将权重矩阵除以你想要的多个子矩阵,并将n * 1.3M矩阵乘以1.3M / n倍的向量。
我认为最快的是迭代次数最少而n会产生适合你内存的最大子矩阵,而不会让你的计算机开始将页面交换到你的硬盘。
你的内存大小应该从n = 5000开始。
你也可以通过使用parfor来加快速度 n
除以处理器数量)。
蛮力方式可能对您有用,只需对其进行一次小优化。
创建权重的^ -0.1操作将比计算加权平均值的+和*操作花费更长的时间,但是您可以在所有百万加权平均操作中重复使用权重。算法变为:
主循环呢 ñ^ 2添加和减法。同 ñ 相当于130万,即3.4万亿次运营。现代3GHz CPU的单核可以说每秒60亿次加法/乘法,因此大约需要10分钟。添加时间索引 weights
矢量和管理费用,我仍然估计你可以在半小时内到达。