## Saturday, February 28, 2009

### Vectorize Matrix Slice Extraction

There was an interesting thread on the Octave mailing lists recently about vectorizing matrix slice extraction. That is avoiding the dreaded for-loop in something like:

%% create a random sparse matrix
n = 10000; %% matrix dimension
nz = 300; %% number of non-zero elements

ia = floor(rand(nz,1)*n); %% dimension 1 indices
ja = floor(rand(nz,1)*n); %% dimension 2 indices
A(:,:) = sparse(ia, ja, rand(nz,1), n, n, n);

%% taking slices with a for loop, slow, poor scaling to large matrices:
for jj = nz:-1:1
slow_columns(:,jj) = A(:,ja(jj));
endfor

%% example using built-in sub2ind, fast, good scaling to large matrices:
jj = sub2ind(size(A)(1),ja);
fast_columns = A(:,jj);

This is a simplified example, the one in the thread shows how you can access a 5-dimensional matrix with 4 indexes. Knowing the built-in functions that come with your high-level language is well worth it, even if it means you have to do weird things to stay in the high-level language and still get reasonable speed. This allows you to avoid slow looping and use the quick compiled functions that someone else has already spent time and effort developing.