以完全矢量化的方式用前一个元素逐行或逐列替换矩阵中的零（或 NAN）答案

【问题标题】：Replacing zeros (or NANs) in a matrix with the previous element row-wise or column-wise in a fully vectorized way以完全矢量化的方式用前一个元素逐行或逐列替换矩阵中的零（或 NAN）
【发布时间】：2014-01-30 20:08:33
【问题描述】：

我需要用前一个元素逐行替换矩阵中的零（或 NaN），所以基本上我需要这个矩阵 X

[0,1,2,2,1,0;  
5,6,3,0,0,2;  
0,0,1,1,0,1]

变成这样：

[0,1,2,2,1,1;  
5,6,3,3,3,2;  
0,0,1,1,1,1],

请注意，如果第一行元素为零，它将保持不变。

我知道这已经以矢量化方式解决了单个行或列向量，这是最好的方法之一：

id = find(X);         
X(id(2:end)) = diff(X(id));       
Y = cumsum(X)

问题在于 Matlab/Octave 中矩阵的索引是连续的，并且按列递增，因此它适用于单行或单列，但不能应用相同的确切概念，但需要用多行进行修改，因为每行raw/column 从新开始，必须被视为独立的。我已经尽力了，用谷歌搜索了整个谷歌，但找不到出路。如果我在循环中应用同样的想法，它会变得太慢，因为我的矩阵至少包含 3000 行。谁能帮我解决这个问题？

【问题讨论】：

好的，我明白你在转置上的意思了。
然而必须再次编辑它以插入 2 个初始零的大小写，并强调它们不应继承前一行的最后一个值这一事实。很抱歉持续编辑，但随着我们的进行，特殊情况和更深入的理解会出现......

标签： matlab loops matrix vectorization

【解决方案1】：

在每一行中隔离零的特殊情况

您可以使用find 的双输出版本来定位除第一列之外的所有列中的零和NaN，然后使用linear indexing 用它们的行前值填充这些条目：

[ii jj] = find( (X(:,2:end)==0) | isnan(X(:,2:end)) );
X(ii+jj*size(X,1)) = X(ii+(jj-1)*size(X,1));

一般情况（每行允许连续的零）

X(isnan(X)) = 0; %// handle NaN's and zeros in a unified way
aux = repmat(2.^(1:size(X,2)), size(X,1), 1) .* ...
    [ones(size(X,1),1) logical(X(:,2:end))]; %// positive powers of 2 or 0
col = floor(log2(cumsum(aux,2))); %// col index
ind = bsxfun(@plus, (col-1)*size(X,1), (1:size(X,1)).'); %'// linear index
Y = X(ind);

诀窍是利用矩阵aux，如果X的对应条目为0且其列号大于1，则该矩阵包含0；否则包含 2 提高到列号。因此，将cumsum 逐行应用于此矩阵，取log2 并向下舍入（矩阵col）为每一行给出最右边的非零条目的列索引，直到当前条目（所以这是一种按行的“累积最大值”函数。）只剩下从列号转换为线性索引（使用bsxfun；也可以使用sub2ind）并使用它来索引X。

这仅对中等大小的X 有效。对于大尺寸，代码使用的 2 的幂很快接近 realmax，并导致不正确的索引。

例子：

X =
     0     1     2     2     1     0     0
     5     6     3     0     0     2     3
     1     1     1     1     0     1     1

给予

>> Y
Y =
     0     1     2     2     1     1     1
     5     6     3     3     3     2     3
     1     1     1     1     1     1     1

【讨论】：

嗨 Luis，谢谢，恐怕我的例子不合适，我现在修好了，你伟大的代码只有在有一个零要替换的情况下才有效，但当有更多时它不会填充其他零：5,6,3,0,0,2；应该变成5,6,3,3,3,2；抱歉使用了一个不恰当的例子。
该解决方案仅适用于最大 n*1023 的矩阵。对于 1024 或更大，aux 包含 inf
路易斯令人印象深刻，不幸的是有一个棘手的问题，如果我生成一个像 X = randint(1000,1000) 这样均匀分布的 1 和 0 的矩阵，它可以工作，但如果我生成 X = randint(1000 ,1000,100) 所以 0 的频率较低，它大多会失败，实际上在 ind 中你可以找到大于 1000000 (1000*1000) 的值。
该死！丹尼尔，你是对的！除了我提到的问题之外，当超过 1024 时，它会在索引中获得 NAN。天哪。看起来这很棘手。
@DanielR 是的。 2^1024 非常接近 realmax

【解决方案2】：

您可以将自己的解决方案概括如下：

Y = X.';                                       %'// Make a transposed copy of X
Y(isnan(Y)) = 0;
idx = find([ones(1, size(X, 1)); Y(2:end, :)]);
Y(idx(2:end)) = diff(Y(idx));
Y = reshape(cumsum(Y(:)), [], size(X, 1)).';   %'// Reshape back into a matrix

它的工作原理是将输入数据视为一个长向量，应用原始解决方案，然后将结果重新整形为矩阵。第一列始终被视为非零，因此这些值不会在整个行中传播。另请注意，原始矩阵已转置，以便以行优先顺序将其转换为向量。

【讨论】：

嗨 Eitan，与 Louis 的评论相同，我可能使用了不正确的示例，代码也应该填充连续的零，我现在更正了，抱歉，感谢您的回答。
明白了。请看一下修改后的答案。
好吧，也许我想通了：它是 Y！
我认为这个解决方案的一个限制是它会跨行传播值。例如，如果一行中的第一个和第二个值为零，则第二个值将获得最新的非零值（来自不同的行），而第一个值显式恢复为零。
@Charles 是的，它有这个确切的问题。谢谢你的展示，我太热情了。

【解决方案3】：

Eitan 的答案的修改版本，以避免跨行传播值：

Y = X'; %'
tf = Y > 0;
tf(1,:) = true;
idx = find(tf);
Y(idx(2:end)) = diff(Y(idx));
Y = reshape(cumsum(Y(:)),fliplr(size(X)))';

【讨论】：

显然，这正是我所做的（没有注意到您发布了第二个解决方案）。请记住，您的虽然容易受到 NaN 的影响。
是的，NaN（或者同样是 Inf 或 -Inf）会为 cumsum 搞砸。如果您的值具有非常不同的大小，浮点错误甚至可能是一个问题。我对这类事情的首选解决方案实际上是我的另一个答案，它使用索引方法。这会将 Infs 处理为常规值，可以应用于非数字类型，支持限制填充距离等。
也就是说 cumsum 方法可能更适合 OP。
很好，对我来说看起来没有问题，而且速度很快：在 3000X3000 矩阵上平均 0.4 秒（在我 4 岁的 i7 上），大约 50% 零。谢谢
@Charles 我一直坚持使用“cumsum 方法”，因为它是由 OP 在问题开始时提出的。

【解决方案4】：

x=[0,1,2,2,1,0;
5,6,3,0,1,2;
1,1,1,1,0,1];
%Do it column by column is easier
x=x';
rm=0;
while 1
    %fields to replace
    l=(x==0);
    %do nothing for the first row/column
    l(1,:)=0;
    rm2=sum(sum(l));
    if rm2==rm
        %nothing to do
        break;
    else
        rm=rm2;
    end
    %replace zeros
    x(l) = x(find(l)-1);
end
x=x';

【讨论】：

这是一个近乎矢量化的解决方案。该循环只需要处理零行，如果您最多有 5 个连续的零，它将迭代 5 次。检查性能，这是一个非常快速的解决方案。
Daniel，谢谢！，我已经尝试过您的解决方案，它实际上速度非常快，比在 1000X1000 矩阵上应用矢量化原始逐行解决方案快大约 100 (!) 倍，在2000X2000。
嗨@Daniel R 我注意到，当一行开头有 2 个或更多连续零时，代码会卡住，对吗？
@aldo：您使用的是我的代码的当前版本吗？写完我的答案后不久，我注意到了一个类似的问题并对其进行了更新。变量rm 和rm2 是新的。
您的解决方案是最快的，只有很少数量的零

【解决方案5】：

我有一个用于填充 NaN 的类似问题的函数。这可能会被进一步缩减或加速 - 它是从具有更多功能（向前/向后填充、最大距离等）的预先存在的代码中提取的。

X = [
    0 1 2 2 1 0
    5 6 3 0 0 2
    1 1 1 1 0 1
    0 0 4 5 3 9
];

X(X == 0) = NaN;
Y = nanfill(X,2);
Y(isnan(Y)) = 0

function y = nanfill(x,dim)
if nargin < 2, dim = 1; end
if dim == 2, y = nanfill(x',1)'; return; end
i = find(~isnan(x(:)));
j = 1:size(x,1):numel(x);
j = j(ones(size(x,1),1),:);
ix = max(rep([1; i],diff([1; i; numel(x) + 1])),j(:));
y = reshape(x(ix),size(x));

function y = rep(x,times)
i = find(times);
if length(i) < length(times), x = x(i); times = times(i); end
i = cumsum([1; times(:)]);
j = zeros(i(end)-1,1);
j(i(1:end-1)) = 1;
y = x(cumsum(j));

【讨论】：

这可能比这个问题需要的更复杂 - 请参阅我的其他答案。