performance difference between subscript indexing and linear indexing

Question

I have a 2D matrix in MATLAB and I use two different ways to access its elements. One is based on subscript indexing and the other is based on linear indexing. I test both methods by following code:

N = 512; it = 400; im = zeros(N);
%// linear indexing
[ind_x,ind_y] = ndgrid(1:2:N,1:2:N);
index = sub2ind(size(im),ind_x,ind_y);

tic
for i=1:it
    im(index) = im(index) + 1;
end
toc %//cost 0.45 seconds on my machine (MATLAB2015b, Thinkpad T410)

%// subscript indexing
x = 1:2:N;
y = 1:2:N;

tic
for i=1:it
    im(x,y) = im(x,y) +1;
end
toc %// cost 0.12 seconds on my machine(MATLAB2015b, Thinkpad T410)

%//someone pointed that double or uint32 might an issue, so we turn both into uint32

%//uint32 for linear indexing
index = uint32(index);
tic
for i=1:it
    im(index) = im(index) +1;
end
toc%// cost 0.25 seconds on my machine(MATLAB2015b, Thinkpad T410)

%//uint32 for the subscript indexing
x = uint32(1:2:N);
y = uint32(1:2:N);
tic
for i=1:it
    im(x,y) = im(x,y) +1;
end
toc%// cost 0.11 seconds on my machine(MATLAB2015b, Thinkpad T410)

%% /*********************comparison with others*****************/
%//third way of indexing, loops
tic
for i=1:it
    for j=1:2:N
        for k=1:2:N
            im(j,k) = im(j,k)+1;
        end
    end
end
toc%// cost 0.74 seconds on my machine(MATLAB2015b, Thinkpad T410)

It seems that directly using subscript indexing is faster than the linear indexing obtained from sub2ind. Does anyone know why? I thought they were almost the same.

Answer 1

The intuition

As Daniel mentioned in his answer, the linear index takes up more space in RAM while the subscripts are much smaller.

For the subscripted indexing, internally, Matlab will not create the linear index, but it will use a (double) compiled loop to cycle through all elements.

The subscripted version on the other hand will have to loop through all the linear indices passed from outside, which will require more reads from memory, thus will take longer.

Claims

• Linear indexing is faster
• ...as long as the total number of indices is the same

Timings

From the timings we see a direct confirmation for the first claim and we can infer the second with some additional testing (below).

LOOPED
      subs assignment: 0.2878s
    linear assignment: 0.0812s

VECTORIZED
      subs assignment: 0.0302s
    linear assignment: 0.0862s

First claim

We can test it with loops. The number of subref operations is the same but the linear index points directly to the element of interest while subscripts, internally, need to be converted.

The functions of interest:

function B = subscriptedIndexing(A,row,col)
n = numel(row);
B = zeros(n);
for r = 1:n
    for c = 1:n
        B(r,c) = A(row(r),col(c));
    end
end
end

function B = linearIndexing(A,index)
B = zeros(size(index));
for ii = 1:numel(index)
    B(ii) = A(index(ii));
end
end

Second claim

This claim is an inference from the observed difference in speed when using the vectorized approach.

First, the vectorized approach (as opposed to the looped) speeds up the subscripted assignment while linear indexing is slightly slower (probably not statistically significant).

Second, the only difference in the two indexing methods comes from the size of the indices/subscripts. We want to isolate this as the only possible cause of the difference in the timings. One other major player could be JIT optimization.

The testing functions:

function B = subscriptedIndexingVect(A,row,col)
n = numel(row);
B = zeros(n);
B = A(row,col);
end

function B = linearIndexingVect(A,index)
B = zeros(size(index));
B = A(index);
end

NOTE: I keep the superfluous preallocation of B, to keep the vectorized and looped approaches comparable. In other words, differences in timings should only come from indexing and the internal implementation of the loops.

All tests are run with:

function testFun(N)
A             = magic(N);
row           = 1:2:N;
col           = 1:2:N;
[ind_x,ind_y] = ndgrid(row,col);
index         = sub2ind(size(A),ind_x,ind_y);

% isequal(linearIndexing(A,index), subscriptedIndexing(A,row,col))
% isequal(linearIndexingVect(A,index), subscriptedIndexingVect(A,row,col))

fprintf('<strong>LOOPED</strong>\n')
fprintf('      subs assignment: %.4fs\n',  timeit(@()subscriptedIndexing(A,row,col)))
fprintf('    linear assignment: %.4fs\n\n',timeit(@()linearIndexing(A,index)))
fprintf('<strong>VECTORIZED</strong>\n')
fprintf('      subs assignment: %.4fs\n',  timeit(@()subscriptedIndexingVect(A,row,col)))
fprintf('    linear assignment: %.4fs\n',  timeit(@()linearIndexingVect(A,index)))
end

Turning JIT on/off has NO impact:

feature accel off
testFun(5e3)
...

VECTORIZED
      subs assignment: 0.0303s
    linear assignment: 0.0873s

feature accel on
testFun(5e3)
...

VECTORIZED
      subs assignment: 0.0303s
    linear assignment: 0.0871s

This excludes that subscripted assignment's superior speed comes from JIT optimization which leaves us with the only plausible cause, number of RAM accesses. It is true that the final matrix has the same number of elements. However, the linear assignment has to retrieve all elements of the index in order to fetch the numbers.

SETUP

Tested on Win7 64 with MATLAB R2015b. Prior versions of Matlab will provide different results due to recent changes in Matlab's execution engine

In fact, turning JIT off in Matlab R2014a affects timings, but only for the loops (expected result):

feature accel off
testFun(5e3)

LOOPED
      subs assignment: 7.8915s
    linear assignment: 6.4418s

VECTORIZED
      subs assignment: 0.0295s
    linear assignment: 0.0878s

This again confirms that the difference in timings between linear and sibscripted assignment should come from the number of RAM accesses, since JIT does not play a role in the vectorized approach.

Answer 2

It does not really surprise me that the subscript indexing is much faster here. If you take a look at your input data, the index is much smaller in this case. For the subscript indexing case you have 512 elements while for the linear indexing case you have 65536 elements.

When you apply your example to a vector instead, you will notice that there is no difference between both methods.

Here is the slightly modified code I used to evaluate different matrix sizes:

it = 400; im = zeros(512*512,1);
x = 1:2:size(im,1);
y = 1:2:size(im,2);
%// linear indexing
[ind_x,ind_y] = ndgrid(x,y);
index = sub2ind(size(im),ind_x,ind_y);

tic
for i=1:it
    im(index) = im(index) + 1;
end
toc 

%// subscript indexing


tic
for i=1:it
    im(x,y) = im(x,y) +1;
end
toc