用先前的非零值替换向量中的所有零
Replace all zeros in vector by previous non-zero value
Matlab/Octave算法示例:
input vector: [ 1 0 2 0 7 7 7 0 5 0 0 0 9 ]
output vector: [ 1 1 2 2 7 7 7 7 5 5 5 5 9 ]
算法非常简单:遍历向量并用最后一个非零值替换所有零。这看起来微不足道,当用一个缓慢的 for (i=1:length) 循环完成并能够引用前一个元素 (i-1) 时也是如此,但看起来不可能以快速矢量化形式表达。
我尝试了 merge() 和 shift() 但它只适用于第一次出现的零,而不是任意数量的零。
是否可以在 Octave/Matlab 中以矢量化形式完成,或者必须使用 C 才能在大量数据上具有足够的性能?
我有 another similar slow for-loop algorithm to speed up 并且似乎通常不可能以矢量化形式引用以前的值,例如 SQL lag()
或 group by
或 loop (i-1)
很容易做到。但是 Octave/Matlab 循环非常慢。
有没有人找到这个一般问题的解决方案,或者出于基本的 Octave/Matlab 设计原因,这是徒劳的吗?
性能基准:
解决方案 1(慢循环)
in = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1 ,100000);
out = in;
tic
for i=2:length(out)
if (out(i)==0)
out(i)=out(i-1);
end
end
toc
[in(1:20); out(1:20)] % test to show side by side if ok
经过的时间是 15.047 秒。
Dan 的解决方案 2(快约 80 倍)
in = V = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1 ,100000);
tic;
d = double(diff([0,V])>0);
d(find(d(2:end))+1) = find(diff([0,~V])==-1) - find(diff([0,~V])==1);
out = V(cumsum(~~V+d)-1);
toc;
[in(1:20); out(1:20)] % shows it works ok
经过的时间是 0.188167 秒。
15.047 / 0.188167 = 79.97 倍改进
GameOfThrows 的解决方案 3(快约 115 倍)
in = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1 ,100000);
a = in;
tic;
pada = [a,888];
b = pada(pada >0);
bb = b(:,1:end-1);
c = find (pada==0);
d = find(pada>0);
len = d(2:end) - (d(1:end-1));
t = accumarray(cumsum([1,len])',1);
out = bb(cumsum(t(1:end-1)));
toc;
经过的时间是 0.130558 秒。
15.047 / 0.130558 = 115.25 倍改进
魔法 Luis Mendo 的解决方案 4 (快约 250 倍)
in = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] , 1, 100000);
tic;
u = nonzeros(in);
out = u(cumsum(in~=0)).';
toc;
经过的时间是 0.0597501 秒。
15.047 / 0.0597501 = 251.83 倍改进
(更新 2019/03/13)使用 MATLAB R2017a 的时间安排:
Slow loop: 0.010862 seconds.
Dan: 0.072561 seconds.
GameOfThrows: 0.066282 seconds.
Luis Mendo: 0.032257 seconds.
fillmissing: 0.053366 seconds.
所以我们再次得出相同的结论:MATLAB 中的循环不再慢!
另请参阅:
Trivial/impossible algorithm challenge in Octave/Matlab Part II: iterations memory
矢量运算通常假设各个项目是独立的。如果您对较早的项目有依赖性,那么循环是最好的方法。
关于 matlab 的一些额外背景知识:在 matlab 中,操作通常更快,不是因为特定的矢量操作,而是因为矢量操作只是在本机 C++ 代码中执行循环而不是通过解释器
我认为是可以的,让我们从基础开始,你想捕获数字大于0的地方:
a = [ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] %//Load in Vector
pada = [a,888]; %//Pad A with a random number at the end to help in case the vector ends with a 0
b = pada(find(pada >0)); %//Find where number if bigger than 0
bb = b(:,1:end-1); %//numbers that are bigger than 0
c = find (pada==0); %//Index where numbers are 0
d = find(pada>0); %//Index where numbers are greater than 0
length = d(2:end) - (d(1:end-1)); %//calculate number of repeats needed for each 0 trailing gap.
%//R = [cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), bb, length,'uniformoutput',0))]; %//Repeat the value
----------EDIT---------
%// Accumarray and cumsum method, although not as nice as Dan's 1 liner
t = accumarray(cumsum([1,length])',1);
R = bb(cumsum(t(1:end-1)));
注意:我使用了 arrayfun
,但您可以使用 accumarray
,因为 well.I 认为这表明可以并行执行此操作?
R =
第 1 列到第 10 列
1 1 2 2 7 7 7 7 5 5
第 11 列到第 13 列
5 5 9
测试:
a = [ 1 0 2 0 7 7 7 0 5 0 0 0 9 0 0 0 ]
R =
第 1 列到第 10 列
1 1 2 2 7 7 7 7 5 5
第 11 至 16 列
5 5 9 9 9 9
表现:
a = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1,10000); %//Double of 130,000
Arrayfun Method : Elapsed time is 6.840973 seconds.
AccumArray Method : Elapsed time is 2.097432 seconds.
我认为是矢量化的解决方案。适用于您的示例:
V = [1 0 2 0 7 7 7 0 5 0 0 0 9]
%// This is where the numbers you will repeat lie. You have to cast to a double otherwise later when you try assign numbers to it it caps them at logical 1s
d = double(diff([0,V])>0)
%// find(diff([0,~V])==-1) - find(diff([0,~V])==1) is the length of each zero cluster
d(find(d(2:end))+1) = find(diff([0,~V])==-1) - find(diff([0,~V])==1)
%// ~~V is the same as V ~= 0
V(cumsum(~~V+d)-1)
下面的简单方法可以满足您的需求,而且速度可能非常快:
in = [1 0 2 0 7 7 7 0 5 0 0 0 9];
t = cumsum(in~=0);
u = nonzeros(in);
out = u(t).';
这是另一个解决方案,使用 linear interpolation with previous neighbor lookup。
我认为它也相当快,因为只有查找和索引,没有计算:
in = [1 0 2 0 7 7 7 0 5 0 0 0 9]
mask = logical(in);
idx = 1:numel(in);
in(~mask) = interp1(idx(mask),in(mask),idx(~mask),'previous');
%// out = in
说明
您需要创建一个索引向量:
idx = 1:numel(in) $// = 1 2 3 4 5 ...
还有一个逻辑掩码,屏蔽所有非零值:
mask = logical(in);
这样你就可以得到用于插值的网格点idx(mask)
和网格数据in(mask)
。查询点 idx(~mask)
是零数据的索引。查询数据 in(~mask)
然后 "calculated" 通过 next previous neighbor 插值,所以它基本上在网格中查看前一个网格点的值是多少。正是你想要的。不幸的是,所涉及的函数对于所有可想到的情况都有巨大的开销,这就是为什么它仍然比 Luis Mendo 的答案慢,尽管不涉及算术计算。
此外,可以稍微减少 interp1
的开销:
F = griddedInterpolant(idx(mask),in(mask),'previous');
in(~mask) = F(idx(~mask));
但是没有太大的影响
in = %// = out
1 1 2 2 7 7 7 7 5 5 5 5 9
基准
0.699347403200000 %// thewaywewalk
1.329058123200000 %// GameOfThrows
0.408333643200000 %// LuisMendo
1.585014923200000 %// Dan
代码
function [t] = bench()
in = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1 ,100000);
% functions to compare
fcns = {
@() thewaywewalk(in);
@() GameOfThrows(in);
@() LuisMendo(in);
@() Dan(in);
};
% timeit
t = zeros(4,1);
for ii = 1:10;
t = t + cellfun(@timeit, fcns);
end
format long
end
function in = thewaywewalk(in)
mask = logical(in);
idx = 1:numel(in);
in(~mask) = interp1(idx(mask),in(mask),idx(~mask),'previous');
end
function out = GameOfThrows(a)
pada = [a,888];
b = pada(find(pada >0));
bb = b(:,1:end-1);
c = find (pada==0);
d = find(pada>0);
length = d(2:end) - (d(1:end-1));
t = accumarray(cumsum([1,length])',1);
out = bb(cumsum(t(1:end-1)));
end
function out = LuisMendo(in)
t = cumsum(in~=0);
u = nonzeros(in);
out = u(t).';
end
function out = Dan(V)
d = double(diff([0,V])>0);
d(find(d(2:end))+1) = find(diff([0,~V])==-1) - find(diff([0,~V])==1);
out = V(cumsum(~~V+d)-1);
end
MATLAB R2016b 中的新功能:fillmissing
,它完全符合问题中的描述:
in = [ 1 0 2 0 7 7 7 0 5 0 0 0 9 ];
in(in==0) = NaN;
out = fillmissing(in,'previous');
[this duplicate question 中发现的新功能]。
Matlab/Octave算法示例:
input vector: [ 1 0 2 0 7 7 7 0 5 0 0 0 9 ]
output vector: [ 1 1 2 2 7 7 7 7 5 5 5 5 9 ]
算法非常简单:遍历向量并用最后一个非零值替换所有零。这看起来微不足道,当用一个缓慢的 for (i=1:length) 循环完成并能够引用前一个元素 (i-1) 时也是如此,但看起来不可能以快速矢量化形式表达。 我尝试了 merge() 和 shift() 但它只适用于第一次出现的零,而不是任意数量的零。
是否可以在 Octave/Matlab 中以矢量化形式完成,或者必须使用 C 才能在大量数据上具有足够的性能?
我有 another similar slow for-loop algorithm to speed up 并且似乎通常不可能以矢量化形式引用以前的值,例如 SQL lag()
或 group by
或 loop (i-1)
很容易做到。但是 Octave/Matlab 循环非常慢。
有没有人找到这个一般问题的解决方案,或者出于基本的 Octave/Matlab 设计原因,这是徒劳的吗?
性能基准:
解决方案 1(慢循环)
in = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1 ,100000);
out = in;
tic
for i=2:length(out)
if (out(i)==0)
out(i)=out(i-1);
end
end
toc
[in(1:20); out(1:20)] % test to show side by side if ok
经过的时间是 15.047 秒。
Dan 的解决方案 2(快约 80 倍)
in = V = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1 ,100000);
tic;
d = double(diff([0,V])>0);
d(find(d(2:end))+1) = find(diff([0,~V])==-1) - find(diff([0,~V])==1);
out = V(cumsum(~~V+d)-1);
toc;
[in(1:20); out(1:20)] % shows it works ok
经过的时间是 0.188167 秒。
15.047 / 0.188167 = 79.97 倍改进
GameOfThrows 的解决方案 3(快约 115 倍)
in = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1 ,100000);
a = in;
tic;
pada = [a,888];
b = pada(pada >0);
bb = b(:,1:end-1);
c = find (pada==0);
d = find(pada>0);
len = d(2:end) - (d(1:end-1));
t = accumarray(cumsum([1,len])',1);
out = bb(cumsum(t(1:end-1)));
toc;
经过的时间是 0.130558 秒。
15.047 / 0.130558 = 115.25 倍改进
魔法 Luis Mendo 的解决方案 4 (快约 250 倍)
in = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] , 1, 100000);
tic;
u = nonzeros(in);
out = u(cumsum(in~=0)).';
toc;
经过的时间是 0.0597501 秒。
15.047 / 0.0597501 = 251.83 倍改进
(更新 2019/03/13)使用 MATLAB R2017a 的时间安排:
Slow loop: 0.010862 seconds.
Dan: 0.072561 seconds.
GameOfThrows: 0.066282 seconds.
Luis Mendo: 0.032257 seconds.
fillmissing: 0.053366 seconds.
所以我们再次得出相同的结论:MATLAB 中的循环不再慢!
另请参阅: Trivial/impossible algorithm challenge in Octave/Matlab Part II: iterations memory
矢量运算通常假设各个项目是独立的。如果您对较早的项目有依赖性,那么循环是最好的方法。
关于 matlab 的一些额外背景知识:在 matlab 中,操作通常更快,不是因为特定的矢量操作,而是因为矢量操作只是在本机 C++ 代码中执行循环而不是通过解释器
我认为是可以的,让我们从基础开始,你想捕获数字大于0的地方:
a = [ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] %//Load in Vector
pada = [a,888]; %//Pad A with a random number at the end to help in case the vector ends with a 0
b = pada(find(pada >0)); %//Find where number if bigger than 0
bb = b(:,1:end-1); %//numbers that are bigger than 0
c = find (pada==0); %//Index where numbers are 0
d = find(pada>0); %//Index where numbers are greater than 0
length = d(2:end) - (d(1:end-1)); %//calculate number of repeats needed for each 0 trailing gap.
%//R = [cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), bb, length,'uniformoutput',0))]; %//Repeat the value
----------EDIT---------
%// Accumarray and cumsum method, although not as nice as Dan's 1 liner
t = accumarray(cumsum([1,length])',1);
R = bb(cumsum(t(1:end-1)));
注意:我使用了 arrayfun
,但您可以使用 accumarray
,因为 well.I 认为这表明可以并行执行此操作?
R =
第 1 列到第 10 列
1 1 2 2 7 7 7 7 5 5
第 11 列到第 13 列
5 5 9
测试:
a = [ 1 0 2 0 7 7 7 0 5 0 0 0 9 0 0 0 ]
R =
第 1 列到第 10 列
1 1 2 2 7 7 7 7 5 5
第 11 至 16 列
5 5 9 9 9 9
表现:
a = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1,10000); %//Double of 130,000
Arrayfun Method : Elapsed time is 6.840973 seconds.
AccumArray Method : Elapsed time is 2.097432 seconds.
我认为是矢量化的解决方案。适用于您的示例:
V = [1 0 2 0 7 7 7 0 5 0 0 0 9]
%// This is where the numbers you will repeat lie. You have to cast to a double otherwise later when you try assign numbers to it it caps them at logical 1s
d = double(diff([0,V])>0)
%// find(diff([0,~V])==-1) - find(diff([0,~V])==1) is the length of each zero cluster
d(find(d(2:end))+1) = find(diff([0,~V])==-1) - find(diff([0,~V])==1)
%// ~~V is the same as V ~= 0
V(cumsum(~~V+d)-1)
下面的简单方法可以满足您的需求,而且速度可能非常快:
in = [1 0 2 0 7 7 7 0 5 0 0 0 9];
t = cumsum(in~=0);
u = nonzeros(in);
out = u(t).';
这是另一个解决方案,使用 linear interpolation with previous neighbor lookup。
我认为它也相当快,因为只有查找和索引,没有计算:
in = [1 0 2 0 7 7 7 0 5 0 0 0 9]
mask = logical(in);
idx = 1:numel(in);
in(~mask) = interp1(idx(mask),in(mask),idx(~mask),'previous');
%// out = in
说明
您需要创建一个索引向量:
idx = 1:numel(in) $// = 1 2 3 4 5 ...
还有一个逻辑掩码,屏蔽所有非零值:
mask = logical(in);
这样你就可以得到用于插值的网格点idx(mask)
和网格数据in(mask)
。查询点 idx(~mask)
是零数据的索引。查询数据 in(~mask)
然后 "calculated" 通过 next previous neighbor 插值,所以它基本上在网格中查看前一个网格点的值是多少。正是你想要的。不幸的是,所涉及的函数对于所有可想到的情况都有巨大的开销,这就是为什么它仍然比 Luis Mendo 的答案慢,尽管不涉及算术计算。
此外,可以稍微减少 interp1
的开销:
F = griddedInterpolant(idx(mask),in(mask),'previous');
in(~mask) = F(idx(~mask));
但是没有太大的影响
in = %// = out
1 1 2 2 7 7 7 7 5 5 5 5 9
基准
0.699347403200000 %// thewaywewalk
1.329058123200000 %// GameOfThrows
0.408333643200000 %// LuisMendo
1.585014923200000 %// Dan
代码
function [t] = bench()
in = repmat([ 1 0 2 0 7 7 7 0 5 0 0 0 9 ] ,1 ,100000);
% functions to compare
fcns = {
@() thewaywewalk(in);
@() GameOfThrows(in);
@() LuisMendo(in);
@() Dan(in);
};
% timeit
t = zeros(4,1);
for ii = 1:10;
t = t + cellfun(@timeit, fcns);
end
format long
end
function in = thewaywewalk(in)
mask = logical(in);
idx = 1:numel(in);
in(~mask) = interp1(idx(mask),in(mask),idx(~mask),'previous');
end
function out = GameOfThrows(a)
pada = [a,888];
b = pada(find(pada >0));
bb = b(:,1:end-1);
c = find (pada==0);
d = find(pada>0);
length = d(2:end) - (d(1:end-1));
t = accumarray(cumsum([1,length])',1);
out = bb(cumsum(t(1:end-1)));
end
function out = LuisMendo(in)
t = cumsum(in~=0);
u = nonzeros(in);
out = u(t).';
end
function out = Dan(V)
d = double(diff([0,V])>0);
d(find(d(2:end))+1) = find(diff([0,~V])==-1) - find(diff([0,~V])==1);
out = V(cumsum(~~V+d)-1);
end
MATLAB R2016b 中的新功能:fillmissing
,它完全符合问题中的描述:
in = [ 1 0 2 0 7 7 7 0 5 0 0 0 9 ];
in(in==0) = NaN;
out = fillmissing(in,'previous');
[this duplicate question 中发现的新功能]。