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.

No comments:

Post a Comment