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

Question

以下是我想要使用的矩阵子集的示例:

1 3 5

2 3 6

1 1 1

3 5 4

5 5 5

8 8 0

该矩阵实际上是3000 x 3.

对于前3行,我希望用这三行中的第一行减去这些行中的每一行.

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

因此,输出矩阵将如下所示:

0 0 0

1 0 1

0 -2 -4

0 0 0

2 0 1

5 3 -4

MATLAB中的哪些代码会为我做这个？

Answer 1

你也可以通过这样做完全矢量化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×3矩阵),它采用3×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. 迪瓦卡的bsxfun方法:0.0689 seconds

因此,Divakar的方法是最快的,其次是David的for循环方法,紧接着是natan的bsxfun方法,接着是natan的原始kron方法,然后是树懒(又名我的).

• 测试#2 - 我决定看看当你增加矩阵的大小时会有多快.设置如下.我做了1000次迭代,并且在每次迭代时,我每次都将矩阵行的大小增加3000.因此,迭代1由3000×3矩阵组成,下一次迭代由6000×3矩阵组成,依此类推.随机种子再次设置为1.在每次迭代时,记录完成代码所花费的时间.为了确保公平性,在处理代码开始之前,每次迭代都清除变量.因此,这是一个stem图表,显示每种矩阵大小的时间.我对绘图进行了子集化,使其显示从200000 x 3到300000 x 3的时序.请注意,水平轴记录每次迭代时的行数.第一个主干为3000行,下一个为6000行,依此类推.列(当然)的列保持不变.

我无法解释整个图表中的随机峰值....可能归因于RAM中发生的事情.但是,我非常确定我在每次迭代时清除变量以确保没有偏差.无论如何,迪瓦卡和大卫都是紧密联系在一起的.接下来是natan的bsxfun方法,然后是natan的kron方法,最后是我的方法.有趣的是看看Divakar的bsxfun方法和大卫的for方法在时间上是如何并排的.

• 测试#3 - 我重复了我为测试#2做的事情,但是使用natan的建议,我决定采用对数比例.我做了6次迭代,从3000 x 3矩阵开始,然后将行增加10倍.因此,第二次迭代有30000 x 3,第三次迭代有300000 x 3,依此类推,直到最后一次迭代,即3e8 x 3.

我在水平轴上绘制了半对数刻度,而垂直轴仍然是线性刻度.同样,横轴表示矩阵中的行数.

我更改了垂直限制,因此我们可以看到大多数方法.我的方法表现很差,以至于它会将其他时间压缩到图表的下端.因此,我改变了观看限制,使我的方法脱离了图片.基本上,在这里验证了测试#2中看到的内容.

Answer 2

这是实现这个的另一种方式,与natan的bsxfun实现bsxfun略有不同-

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

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,并且这在大数组时会发生变化(比率<1表示kron在给定大小的时候花费的时间更少数组).这是在Matlab 2012a win7(i5机器)完成的