3

I have a 7*1 vector a = (1:7).'. I want to form a matrix A of size 4*4 from vector a such that the elements of a form the anti-diagonals of matrix A as follows:

A = [1 2 3 4;
     2 3 4 5;
     3 4 5 6;
     4 5 6 7]

I would like this to work for a general a, not just when the elements are consecutive integers.

I appreciate any help.

Wolfie
  • 27,562
  • 7
  • 28
  • 55
sn at
  • 61
  • 7

2 Answers2

6

Setting up the indexing

Adding the two outputs of meshgrid can give the indices:

[x, y] = meshgrid(1:4, 0:3);
x + y;
% ans = [1     2     3     4
%        2     3     4     5
%        3     4     5     6
%        4     5     6     7];

If a was just as in your example you could stop there. Alternatively, use this to index a general vector a. For comparison, I'll use the same example input as rahnema1 did for their method:

a = [4 6 2 7 3 5 1];
[x, y] = meshgrid(1:4, 0:3);
A = a(x + y);    
% A = [4     6     2     7
%      6     2     7     3
%      2     7     3     5
%      7     3     5     1] 

There are many ways you could create the indices instead of using meshgrid, see the benchmarking functions below for some exampels!


Benchmarking and seven different methods.

Here are some timings for running different methods, including methods using cumsum, repmat, hankel and a simple for loop. This benchmark was done in Matlab 2015b, so takes advantage of Matlab optimisations etc. which the Octave benchmarks in rahnema1's answer might not do. I also use the timeit function which is more robust than tic/toc because it does multiple trials etc.

function benchie()
    n = 10000;    % (large) square matrix size
    a = 1:2*n-1;  % array of correct size, could be anything this long
    f1 = @() m1(a,n);  disp(['bsxfun:   ', num2str(timeit(f1))]);
    f2 = @() m2(a,n);  disp(['cumsum:   ', num2str(timeit(f2))]);
    f3 = @() m3(a,n);  disp(['meshgrid: ', num2str(timeit(f3))]);
    f4 = @() m4(a,n);  disp(['repmat:   ', num2str(timeit(f4))]);
    f5 = @() m5(a,n);  disp(['for loop: ', num2str(timeit(f5))]);
    f6 = @() m6(a,n);  disp(['hankel1:  ', num2str(timeit(f6))]); 
    f7 = @() m7(a,n);  disp(['hankel2:  ', num2str(timeit(f7))]); 
end
% Use bsxfun to do broadcasting of addition
function m1(a,n); A = a(bsxfun(@plus, (1:n), (0:n-1).')); end
% Use cumsum to do cumulative vertical addition to create indices
function m2(a,n); A = a(cumsum([(1:n); ones(n-1,n)])); end
% Add the two meshgrid outputs to get indices
function m3(a,n); [x, y] = meshgrid(1:n, 0:n-1); A = a(x + y); end
% Use repmat twice to replicate the meshgrid results, for equivalent one liner
function m4(a,n); A = a(repmat((1:n)',1,n) + repmat(0:n-1,n,1)); end
% Use a simple for loop. Initialise A and assign values to each row in turn
function m5(a,n); A = zeros(n); for ii = 1:n; A(:,ii) = a(ii:ii+n-1); end; end
% Create a Hankel matrix (constant along anti-diagonals) for indexing
function m6(a,n); A = a(hankel(1:n,n:2*n-1)); end
% Create a Hankel matrix directly from elements
function m7(a,n); A = hankel(a(1:n),a(n:2*n-1)); end

Output:

bsxfun:   1.4397 sec
cumsum:   2.0563 sec
meshgrid: 2.0169 sec
repmat:   1.8598 sec
for loop: 0.4953 sec % MUCH quicker!
hankel1:  2.6154 sec
hankel2:  1.4235 sec

So you are best off using rahnema1's suggestion of bsxfun or direct generation of a hankel matrix if you want a one liner, here is a brilliant StackOverflow answer which explains some of bsxfun's advantages: In Matlab, when is it optimal to use bsxfun?

However, the for loop is more than twice as quick! Conclusion: Matlab has lots of neat ways to achieve things like this, sometimes a simple for loop with some appropriate pre-allocation and Matlab's internal optimisations can be the quickest.

Wolfie
  • 27,562
  • 7
  • 28
  • 55
  • @snat I recommend trying the above benchmark to see which is quickest on your machine and version of Matlab! Mileage may vary – Wolfie Jun 23 '17 at 14:28
5

You can use hankel:

n= 4;
A= hankel(a(1:n),a(n:2*n-1))

Other solution(expansion/bsxfun):

In MATLAB r2016b /Octave It can be created as:

A = a((1:4)+(0:3).')

In pre r2016b you can use bsxfun:

A = a(bsxfun(@plus,1:4, (0:3).'))

Example input/output

a = [4 6 2 7 3 5 1]

A =

   4   6   2   7
   6   2   7   3
   2   7   3   5
   7   3   5   1

Using the benchmark provided by @Wolfie tested in Octave:

 _____________________________________
|Method   |memory peak(MB)|timing(Sec)|
|=========|===============|===========|
|bsxfun   |2030           |1.50       |
|meshgrid |3556           |2.43       |
|repmat   |2411           |2.64       |
|hankel   |886            |0.43       |
|for loop |886            |0.82       |
rahnema1
  • 15,264
  • 3
  • 15
  • 27