4

I have a p-by-p-by-n tensor. I want to extract diagonal element for each p-by-p slice. Are there anyone know how to do this without looping?

Thank you.

1
  • 1
    Why do you want to do it without looping? How big will p and n be?
    – David
    Feb 9, 2015 at 3:23

3 Answers 3

9

Behold the ever mighty and ever powerful bsxfun for vectorizing MATLAB problems to do this task very efficiently using MATLAB's linear indexing -

diags = A(bsxfun(@plus,[1:p+1:p*p]',[0:n-1]*p*p))

Sample run with 4 x 4 x 3 sized input array -

A(:,:,1) =
    0.7094    0.6551    0.9597    0.7513
    0.7547    0.1626    0.3404    0.2551
    0.2760    0.1190    0.5853    0.5060
    0.6797    0.4984    0.2238    0.6991
A(:,:,2) =
    0.8909    0.1493    0.8143    0.1966
    0.9593    0.2575    0.2435    0.2511
    0.5472    0.8407    0.9293    0.6160
    0.1386    0.2543    0.3500    0.4733
A(:,:,3) =
    0.3517    0.9172    0.3804    0.5308
    0.8308    0.2858    0.5678    0.7792
    0.5853    0.7572    0.0759    0.9340
    0.5497    0.7537    0.0540    0.1299
diags =
    0.7094    0.8909    0.3517
    0.1626    0.2575    0.2858
    0.5853    0.9293    0.0759
    0.6991    0.4733    0.1299

Benchmarking

Here's few runtime tests comparing this bsxfun based approach against repmat + eye based approach for big datasizes -

***** Datasize: 500 x 500 x 500 *****
----------------------- With BSXFUN
Elapsed time is 0.008383 seconds.
----------------------- With REPMAT + EYE
Elapsed time is 0.163341 seconds.

***** Datasize: 800 x 800 x 500 *****
----------------------- With BSXFUN
Elapsed time is 0.012977 seconds.
----------------------- With REPMAT + EYE
Elapsed time is 0.402111 seconds.

***** Datasize: 1000 x 1000 x 500 *****
----------------------- With BSXFUN
Elapsed time is 0.017058 seconds.
----------------------- With REPMAT + EYE
Elapsed time is 0.690199 seconds.
5
  • Good pet! Very nice and obedient. +1.
    – rayryeng
    Feb 9, 2015 at 5:05
  • Direct indexing instead of masking. Lovely :) May 16, 2015 at 10:01
  • 2
    @SanthanSalai Yup! Lot less memory requirements and still vectorized, just one of those special cases.
    – Divakar
    May 16, 2015 at 10:06
  • Any body seeing these kind of prob would obviously think of masking but you took a different route ;) May 16, 2015 at 10:12
  • 2
    @SanthanSalai I didn't take a different route, I took my instinctively default route, the bsxfun route! ;)
    – Divakar
    May 16, 2015 at 12:58
6

One suggestion I have is to create a p x p logical identity matrix, replicate this n times in the third dimension, and then use this matrix to access your tensor. Something like this, supposing that your tensor was stored in A:

ind = repmat(logical(eye(p)), [1 1 n]);
out = A(ind);

Example use:

>> p = 5; n = 3;
>> A = reshape(1:75, p, p, n)

A(:,:,1) =

     1     6    11    16    21
     2     7    12    17    22
     3     8    13    18    23
     4     9    14    19    24
     5    10    15    20    25


A(:,:,2) =

    26    31    36    41    46
    27    32    37    42    47
    28    33    38    43    48
    29    34    39    44    49
    30    35    40    45    50


A(:,:,3) =

    51    56    61    66    71
    52    57    62    67    72
    53    58    63    68    73
    54    59    64    69    74
    55    60    65    70    75

>> ind = repmat(logical(eye(p)), [1 1 n]);
>> out = A(ind)

out =

     1
     7
    13
    19
    25
    26
    32
    38
    44
    50
    51
    57
    63
    69
    75

You'll notice that we grab the diagonals of the first slice, followed by the diagonals of the second slice, etc. up until the last slice. These values are all concatenated into a single vector.

1
  • 1
    I used my "pet" here too! ;)
    – Divakar
    Feb 9, 2015 at 4:57
2

Just reading Divakar's answer and trying to understand why he again is roughly 10 times faster than my idea I put together code mixing both, and ended up with code which is faster:

A=reshape(A,[],n);
diags2 = A(1:p+1:p*p,:);

For a 500x500x500 tensor I get 0.008s for Divakar's solution and 0.005s for my solution using Matlab 2013a. Probably plain indexing is the only way to beat bsxfun.

5
  • 1
    Ah involves reshaping to index, nice. Won't that reshape be a bit costly affair?
    – Divakar
    Mar 17, 2016 at 13:23
  • Isn't reshape basically zero cost? I never saw it adding any substantial runtime.
    – Daniel
    Mar 17, 2016 at 13:26
  • Well, reshape isn't exactly costly when you compare with something like permute, but not free for sure. Did you include reshape in the timings?
    – Divakar
    Mar 17, 2016 at 13:28
  • Yes, it was included into my timings. A reshape does not touch the data in memory at all.
    – Daniel
    Mar 17, 2016 at 13:30
  • @Divakar: Had a similar discussion with Ander Biguri yesterday in the chat. Now realizing that even you aren't aware of the speed of reshape, I think I should share that knowledge
    – Daniel
    Mar 17, 2016 at 13:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.