What is the most efficient way to extract an index from each column of a matrix without using a loop?

I would like to be able to extract a vector 'y' from an n-by-n matrix 'A', taking one element from each column of A, with the vector 'x' specifying the elements to be extracted. More specifically, I want to extract y(i)= A(x(i), i) for i= 1, 2, …, without using a loop. Currently, the only mechanism that I know of for doing this is the following:

A= magic(5);

x= [3 3 3 4 1]; 

y= diag(A(x,1:5))

The disadvantage of the above is that the full matrix corresponding to the outer product of x and 1:n is generated as an intermediate result. In the above example this is not a concern, but in large scale problems this is very inefficient. I would like to know if there is a better way to go about this.

Answer 1

[MathWorks Support Team]

In this problem, performance improvements could be gained by the use of linear indexing, rather than DIAG or a loop to reference the row indices of matrix columns.

On some platforms, the following code snippet yields a more-than-halving of time taken to execute the matrix reference done using DIAG in the problem statement, while also having the benefit of being easier to read as code.

As an example, the following script illustrates the implementation of linear indexing.

A= magic(5);

% The matrix A can also be considered the same as this output when using linear % indexing:

%      A=reshape(A,[1 25]); 

x= [3 3 3 4 1]; 

xl = [(0:5:24)] + x;  % Linear indices

disp('Starting linear indexing timing...')

tic

for i=1:1e5

      yr = A(xl);

end

tResh = toc

A comparison of the two index referencing methods yields better performance in the case of linear indexing:

     Starting DIAG timing...

     tDiag =

         0.2301

     Starting linear indexing timing...

     tResh =

         0.0910