用前一个非零值替换向量中的所有零

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 倍改进

Magical 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

Answer 1

向量运算通常假设各个项目是独立的。如果您依赖于较早的项目，那么循环是最好的方法。

关于 matlab 的一些额外背景：在 matlab 中，运算通常更快，不是因为专门的向量运算，而是因为向量运算只是在本机 C++ 代码中而不是通过解释器进行循环

Answer 2

我觉得是可以的，我们先从基础开始，你想捕捉数字大于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。我认为这表明可以并行执行此操作？

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.

Answer 3

我认为是一个矢量化的解决方案。适用于您的示例：

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)

Answer 4

以下简单的方法可以满足您的需求，并且可能非常快：

in = [1 0 2 0 7 7 7 0 5 0 0 0 9];
t = cumsum(in~=0);
u = nonzeros(in);
out = u(t).';

Answer 5

这是另一个解决方案，使用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) 然后通过下一个前一个邻居插值“计算”，所以它基本上在网格中查看前一个网格点的值是多少。正是你想要的。不幸的是，所涉及的函数对于所有可以想到的情况都有巨大的开销，这就是为什么它仍然比 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

Answer 6

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 中发现的这项新功能]。