MATLAB：按特定行减去矩阵子集

Question

这是我想使用的矩阵子集的示例：

1 3 5

2 3 6

1 1 1

3 5 4

5 5 5

8 8 0

这个矩阵实际上是 3000 x 3。

对于前 3 行，我希望用这三行中的第一行减去每一行。

对于后 3 行，我希望用这三行中的第一行减去每一行，依此类推。

因此，输出矩阵将如下所示：

0 0 0

1 0 1

0 -2 -4

0 0 0

2 0 1

5 3 -4

MATLAB 中的哪些代码会为我完成这项工作？

Answer 1

您可以只使用索引来做到这一点：

a(:) = a(:) - a(3*floor((0:numel(a)-1)/3)+1).';

当然，上面的3 可以替换为任何其他数字。即使该数字不除以行数，它也有效。

Answer 2

您也可以通过使用mat2cell、cellfun，然后使用cell2mat 来完全矢量化。假设我们的矩阵存储在A，尝试：

numBlocks = size(A,1) / 3;
B = mat2cell(A, 3*ones(1,numBlocks), 3);
C = cellfun(@(x) x - x([1 1 1], :), B, 'UniformOutput', false);
D = cell2mat(C); %//Output

第一行计算出我们需要多少个 3 x 3 块。这是假设行数是 3 的倍数。第二行使用mat2cell 分解每个 3 x 3 块并将它们放入单独的单元格中。然后第三行使用cellfun，这样对于我们的元胞数组（这是一个 3 x 3 矩阵）中的每个元胞，它获取 3 x 3 矩阵的每一行并用第一行减去自身。这与@David 所做的非常相似，只是我没有使用repmat 来最小化开销。第四行然后获取这些矩阵中的每一个并将它们堆叠回来，以便我们最终得到最终的矩阵。

示例（这是使用您帖子中定义的矩阵）：

A = [1 3 5; 2 3 6; 1 1 1; 3 5 4; 5 5 5; 8 8 0];
numBlocks = size(A,1) / 3;
B = mat2cell(A, 3*ones(1, numBlocks), 3);
C = cellfun(@(x) x - x([1 1 1], :), B, 'UniformOutput', false);
D = cell2mat(C);

输出：

D =

 0     0     0
 1     0     1
 0    -2    -4
 0     0     0
 2     0     1
 5     3    -4

事后看来，我认为@David 在性能提升方面是正确的。除非这段代码重复多次，否则我认为for 循环会更有效率。不管怎样，我想提供另一种选择。很酷的运动！

编辑：时间和大小测试

由于我们之前的讨论，我决定进行时间和尺寸测试。这些测试是在具有 16 GB RAM 的 Intel i7-4770 @ 3.40 GHz CPU 上执行的，在 Windows 7 Ultimate 上使用 MATLAB R2014a。基本上，我做了以下事情：

• 测试 #1 - 将随机种子生成器设置为 1 以获得可重复性。我写了一个循环 10000 次的循环。对于循环中的每次迭代，我生成一个随机整数 3000 x 3 矩阵，然后执行此处描述的每个方法。我记下每种方法在 10000 次循环后完成所需的时间。计时结果为：

  1. 大卫的方法：0.092129 seconds
  2. rayryeng的方法：1.9828 seconds
  3. natan的方法：0.20097 seconds
  4. natan 的bsxfun 方法：0.10972 seconds
  5. Divakar 的bsxfun 方法：0.0689 seconds

因此，Divakar 的方法是最快的，其次是 David 的 for 循环方法，紧随其后的是 natan 的 bsxfun 方法，其次是 natan 的原始 kron 方法，其次是树懒（又名我的）。

• 测试 #2 - 我决定看看随着矩阵大小的增加，它会变得多快。设置如下。我进行了 1000 次迭代，在每次迭代中，我每次将矩阵行的大小增加 3000。因此，迭代 1 由 3000 x 3 矩阵组成，下一次迭代由 6000 x 3 矩阵组成，依此类推。随机种子再次设置为 1。在每次迭代中，都会记录完成代码所花费的时间。为了确保公平，在处理代码开始之前，每次迭代都会清除变量。因此，这是一个stem 图，它向您展示了每种矩阵大小的时间。我对绘图进行了子集化，以便它显示从 200000 x 3 到 300000 x 3 的时间。请注意，水平轴记录了每次迭代的行数。第一个词干用于 3000 行，下一个词干用于 6000 行，依此类推。列在 3 处保持不变（当然）。

我无法解释整个图表中的随机峰值......可能归因于 RAM 中发生的某些事情。但是，我很确定我会在每次迭代中清除变量以确保没有偏差。无论如何，Divakar 和 David 关系密切。接下来是 natan 的 bsxfun 方法，然后是 natan 的 kron 方法，最后是我的方法。有趣的是，看看 Divakar 的 bsxfun 方法和 David 的 for 方法在时间上是如何并排的。

• 测试 #3 - 我重复了测试 #2 的操作，但根据 natan 的建议，我决定采用对数刻度。我进行了 6 次迭代，从 3000 x 3 矩阵开始，然后将行数增加 10 倍。因此，第二次迭代有 30000 x 3，第三次迭代有 300000 x 3，依此类推，直到最后一次迭代，即 3e8 x 3。

我在水平轴上绘制了半对数刻度，而垂直轴仍然是线性刻度。同样，水平轴描述了矩阵中的行数。

我更改了垂直限制，以便我们可以看到大多数方法。我的方法性能很差，以至于它会将其他时间压向图表的下端。因此，我更改了查看限制以将我的方法排除在图片之外。基本上在测试 #2 中看到的内容在这里得到验证。

Answer 3

一种略短且矢量化的方式将是（如果a 是您的矩阵）：

b=a-kron(a(1:3:end,:),ones(3,1));

让我们测试一下：

a=[1 3 5
   2 3 6
   1 1 1
   3 5 4
   5 5 5
   8 8 0]

a-kron(a(1:3:end,:),ones(3,1))

ans =
 0     0     0
 1     0     1
 0    -2    -4
 0     0     0
 2     0     1
 5     3    -4

编辑

这是一个 bsxfun 解决方案（不太优雅，但希望更快）：

a-reshape(bsxfun(@times,ones(1,3),permute(a(1:3:end,:),[2 3 1])),3,[])'

ans =

 0     0     0
 1     0     1
 0    -2    -4
 0     0     0
 2     0     1
 5     3    -4

编辑 2

好的，这让我感到好奇，因为我知道 bsxfun 对于更大的数组大小开始效率降低。所以我尝试使用timeit 检查我的两个解决方案（因为它们是一个衬垫，这很容易）。这里是：

range=3*round(logspace(1,6,200));
for n=1:numel(range)
    a=rand(range(n),3);
    f=@()a-kron(a(1:3:end,:),ones(3,1));
    g=@() a-reshape(bsxfun(@times,ones(1,3),permute(a(1:3:end,:),[2 3 1])),3,[])';
    t1(n)=timeit(f);
    t2(n)=timeit(g);
end
semilogx(range,t1./t2);

所以我没有测试 for 循环和 Divkar 的 bsxfun，但是你可以看到对于小于 3e4 的数组，kron 比 bsxfun 更好，并且这在更大的数组中发生了变化（比率

Answer 4

这是使用bsxfun 实现此功能的另一种方法，与natan's bsxfun implementation 略有不同-

t1 = reshape(a,3,[]); %// a is the input matrix
out = reshape(bsxfun(@minus,t1,t1(1,:)),[],3); %// Desired output

Answer 5

简单的for 循环。这会分别执行每个 3x3 块。

A=randi(5,9,3)
B=A(1:3:end,:)
for i=1:length(A(:,1))/3
    D(3*i-2:3*i,:)=A(3*i-2:3*i,:)-repmat(B(i,:),3,1)
end
D

虽然有可能对此进行矢量化，但我认为性能提升并不值得，除非您会多次这样做。对于 3000x3 矩阵，它根本不需要很长时间。

编辑：事实上这似乎相当快。我认为这是因为 Matlab 的 JIT 编译可以很好地加速简单的for 循环。