MATLAB: Subtracting matrix subsets by specific rows

Question

Here is an example of a subset of the matrix I would like to use:

1 3 5

2 3 6

1 1 1

3 5 4

5 5 5

8 8 0

This matrix is in fact 3000 x 3.

For the first 3 rows, I wish to subtract each of these rows with the first row of these three.

For the second 3 rows, I wish to subtract each of these rows with the first of these three, and so on.

As such, the output matrix will look like:

0 0 0

1 0 1

0 -2 -4

0 0 0

2 0 1

5 3 -4

What code in MATLAB will do this for me?

Answer 1

a slightly shorter and vectorized way will be (if a is your matrix) :

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

let's test:

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

Edit

Here's a bsxfun solution (less elegant, but hopefully faster):

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

Edit 2

Ok, this got me curios, as I know bsxfun starts to be less efficient for bigger array sizes. So I tried to check using timeit my two solutions (because they are one liners it's easy). And here it is:

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);

So I didn't test for the for loop and Divkar's bsxfun, but you can see that for arrays smaller than 3e4 kron is better than bsxfun, and this changes at larger arrays (ratio of <1 means kron took less time given the size of the array). This was done at Matlab 2012a win7 (i5 machine)

Answer 2

You could also do this completely vectorized by using mat2cell, cellfun, then cell2mat. Assuming our matrix is stored in A, try:

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

The first line figures out how many 3 x 3 blocks we need. This is assuming that the number of rows is a multiple of 3. The second line uses mat2cell to decompose each 3 x 3 block and places them into individual cells. The third line then uses cellfun so that for each cell in our cell array (which is a 3 x 3 matrix), it takes each row of the 3 x 3 matrix and subtracts itself with the first row. This is very much like what @David did, except I didn't use repmat to minimize overhead. The fourth line then takes each of these matrices and stacks them back so that we get our final matrix in the end.

Example (this is using the matrix that was defined in your post):

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);

Output:

D =

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

In hindsight, I think @David is right with respect to performance gains. Unless this code is repeated many times, I think the for loop will be more efficient. Either way, I wanted to provide another alternative. Cool exercise!

Edit: Timing and Size Tests

Because of our discussion earlier, I have decided to do timing and size tests. These tests were performed on an Intel i7-4770 @ 3.40 GHz CPU with 16 GB of RAM, using MATLAB R2014a on Windows 7 Ultimate. Basically, I did the following:

• Test #1 - Set the random seed generator to 1 for reproducibility. I wrote a loop that cycled 10000 times. For each iteration in the loop, I generated a random integer 3000 x 3 matrix, then performed each of the methods that were described here. I took note of how long it took for each method to complete after 10000 cycles. The timing results are:

  1. David's method: 0.092129 seconds
  2. rayryeng's method: 1.9828 seconds
  3. natan's method: 0.20097 seconds
  4. natan's bsxfun method: 0.10972 seconds
  5. Divakar's bsxfun method: 0.0689 seconds

As such, Divakar's method is the fastest, followed by David's for loop method, followed closely by natan's bsxfun method, followed by natan's original kron method, followed by the sloth (a.k.a mine).

• Test #2 - I decided to see how fast this would get as you increase the size of the matrix. The set up was as follows. I did 1000 iterations, and at each iteration, I increase the size of the matrix rows by 3000 each time. As such, iteration 1 consisted of a 3000 x 3 matrix, the next iteration consisted of a 6000 x 3 matrix and so on. The random seed was set to 1 again. At each iteration, the time taken to complete the code was taken a note of. To ensure fairness, the variables were cleared at each iteration before the processing code began. As such, here is a stem plot that shows you the timing for each size of matrix. I subsetted the plot so that it displays timings from 200000 x 3 to 300000 x 3. Take note that the horizontal axis records the number of rows at each iteration. The first stem is for 3000 rows, the next is for 6000 rows and so on. The columns remain the same at 3 (of course).

I can't explain the random spikes throughout the graph.... probably attributed to something happening in RAM. However, I'm very sure I'm clearing the variables at each iteration to ensure no bias. In any case, Divakar and David are closely tied. Next comes natan's bsxfun method, then natan's kron method, followed last by mine. Interesting to see how Divakar's bsxfun method and David's for method are side-by-side in timing.

• Test #3 - I repeated what I did for Test #2, but using natan's suggestion, I decided to go on a logarithmic scale. I did 6 iterations, starting at a 3000 x 3 matrix, and increasing the rows by 10 fold after. As such, the second iteration had 30000 x 3, the third iteration had 300000 x 3 and so on, up until the last iteration, which is 3e8 x 3.

I have plotted on a semi-logarithmic scale on the horizontal axis, while the vertical axis is still a linear scale. Again, the horizontal axis describes the number of rows in the matrix.

I changed the vertical limits so we can see most of the methods. My method is so poor performing that it would squash the other timings towards the lower end of the graph. As such, I changed the viewing limits to take my method out of the picture. Essentially what was seen in Test #2 is verified here.

Answer 3

Here's another way to implement this with bsxfun, slightly different from natan's bsxfun implementation -

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

Answer 4

Simple for loop. This does each 3x3 block separately.

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

Whilst it may be possible to vectorise this, I don't think the performance gains would be worth it, unless you will do this many times. For a 3000x3 matrix it doesn't take long at all.

Edit: In fact this seems to be pretty fast. I think that's because Matlab's JIT compilation can speed up simple for loops well.

Answer 5

You can do it using just indexing:

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

Of course, the 3 above can be replaced by any other number. It works even if that number doesn't divide the number of rows.