Obscure Matlab features #2: linear indexing, sub2ind, ind2sub

Matlab supports 3 forms of indexing matrices. There is the familiar range indexing that selects a submatrix within a larger matrix:

A = randn(10);
 A(1:2,1:4) %a 2x4 matrix extracted from the top of A

This is also called sub-indexing. There is also logical indexing:

A  = randn(10);
 A(A>1.5)  = Inf; %Transform values greater than 1.5 to Inf.

These two methods are fairly expressive, yet there are a certain number of operations which are awkward or impossible to perform with these. Let’s say that we have a matrix A and we want to set its off diagonal (not the regular \ one, the one shaped like /) to all ones. You might think this requires the use of a loop or some special function within Matlab. Yet this operation is easily done in a single step with linear indexing.

Background: a Matlab matrix is internally stored as a flag defining its shape, followed by sequential numeric data organized such that the first index of the matrix grows fastest and the last grows slowest. Thus, in a 2d matrix, the elements are ordered in a И pattern. If a matrix has size (10,20), the element (1,1) has linear index 1, the element (10,1) has index 10, (1,2) has index 11, and so forth.

To target the off-diagonal of a 10×10 matrix, we write down the paired indices of the off diagonal:

(10,1), (9,2), (8,3), …

And translate them into linear indices. For a 2D matrix of size MxN, the element (i,j) has linear index equal to i + (j-1)*M. So you could select the off diagonal using:

A((10:-1:1) + (0:9)*10)

This is not very readable, however. The functions sub2ind and ind2sub can translate linear indices to subindices and vice-versa. So, to set the off-diagonal of a matrix to all ones, we can use:

A = randn(10);
 A(sub2ind(size(A),size(A,1):-1:1,1:size(A,2))) = 1;

This come in handy in other situations. Say that you have a vector of parameters for a model that represent spatial gains. So for a 16×16 receptive field you have represented the data as a 256×1 vector w. Say that you want to create a matrix D that computes horizonal spatial derivatives. Here’s one way to do it that involves no for loop:

inds = sub2ind([16,16],reshape(repmat((1:16),15*2,1),2,16*15),repmat([1:15;2:16],1,16));
 D = zeros(15*16,16*16);
 D(sub2ind([15*16,16*16],ones(2,1)*(1:15*16),inds)) = repmat([1,-1]',1,15*16);

Unreadable, but will net you a mighty amount of geek points! Now for something a little more useful. Say that you want to find the location of the maximum element in a matrix. You might use something like:

A = randn(10);
 [rowmax,bestrows] = max(A);
 [~,bestcol] = max(bestrows)
 %The max element is at (bestrows(bestcol),bestcol)

It’s cleaner however to use the fact that the (:) operator reshapes a vector or matrix A to size (numel(A),1), obtain the max position as a linear index, and then translate to sub-indices:

A = randn(10);
 [~,bestpos] = max(A(:));
 [bestrow,bestcol] = ind2sub(size(A),bestpos);

I should also note that this works regardless of how many dimensions A has. So for example the largest element in a 5-dimensional matrix can be determined like so:

A = randn(5,4,3,4,6);
 [~,bestpos] = max(A(:));
 [mi,mj,mk,ml,mm] = ind2sub(size(A),bestpos);

sub2ind and ind2sub thus can make it quite a bit easier to work with multidimensional matrices.

Leave a comment

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s