strange performance behavior - microbenchmark

There are two implementations of the same algorithm:

function B = alg1( B, R, alpha )
% Fast computation of L1 distance transform
K = numel(B);
% forward pass
for k=2:K
    B(k) = min( B(k), B(k-1) + alpha * (R(k) - R(k-1)));
end
% backward pass
for k=K-1:-1:1
    B(k) = min( B(k), B(k+1) + alpha * (R(k+1) - R(k)));
end
end

and

function B = alg2( B, R, alpha )
% Fast computation of L1 distance transform
alphaRdiff = alpha*diff(R);
K = numel(B);
% forward pass
for k=2:K
    B(k) = min( B(k), B(k-1) + alphaRdiff(k-1));
end
% backward pass
for k=K-1:-1:1
    B(k) = min( B(k), B(k+1) + alphaRdiff(k));
end
end

These two algorithms differs only by elimination of term

alpha * (R(k) - R(k-1))) 

from the both for-loop.

And by scaled differences pre-computation of vector R

alphaRdiff = alpha*diff(R);

So, the second algorithm alg2 should be faster, because only half of multiplications alpha*(R(k)-R(k-1)) is performed.

But speed of both algorithm is still nearly same (or unmodified alg 1 is even faster), see

N = 1e5;
B = rand(1,N);
R = rand(1,N);

tic;A = alg1( B, R, 1/3 );toc
tic;A_= alg2( B, R, 1/3 );toc
isequal(A,A_)

Elapsed time is 0.037654 seconds.
Elapsed time is 0.029562 seconds.

ans =
logical
1

N = 1e8;
tic;A = alg1( B, R, 1/3 );toc
tic;A_= alg2( B, R, 1/3 );toc

Elapsed time is 1.444549 seconds.
Elapsed time is 1.563630 seconds.

What is the explanation of this strange behavior?

Answer 1

[Bruno Luong]

My guess is in the first code, the intermediate result does have to be stored in memory. It can be in processor register or in the first level of cache. So even the operation is performed twice it is still faster.

I made an algo_3 & 4 that replace min by if, and they seem even faster.

N = 1e8;
B = rand(1,N);
R = rand(1,N);
tic;A = alg1( B, R, 1/3 );toc % Elapsed time is 0.819863 seconds.
tic;A= alg2( B, R, 1/3 );toc  % Elapsed time is 0.852902 seconds.
tic;A= alg3( B, R, 1/3 );toc  % Elapsed time is 0.619154 seconds.
tic;A= alg4( B, R, 1/3 );toc  % Elapsed time is 0.564323 seconds.

function B = alg1( B, R, alpha )
% Fast computation of L1 distance transform
K = numel(B);
% forward pass
for k=2:K
    B(k) = min( B(k), B(k-1) + alpha * (R(k) - R(k-1)));
end
% backward pass
for k=K-1:-1:1
    B(k) = min( B(k), B(k+1) + alpha * (R(k+1) - R(k)));
end
end

function B = alg2( B, R, alpha )
% Fast computation of L1 distance transform
alphaRdiff = alpha*diff(R);
K = numel(B);
% forward pass
for k=2:K
    B(k) = min( B(k), B(k-1) + alphaRdiff(k-1));
end
% backward pass
for k=K-1:-1:1
    B(k) = min( B(k), B(k+1) + alphaRdiff(k));
end
end

function B = alg3( B, R, alpha )
% Fast computation of L1 distance transform
alphaRdiff = alpha*diff(R,1,2);
K = numel(B);
% forward pass
for k=2:K
    C = B(k-1) + alphaRdiff(k-1);
    if C < B(k)
        B(k) = C;
    end
end
% backward pass
for k=K-1:-1:1
    C = B(k+1) + alphaRdiff(k);
    if C < B(k)
        B(k) = C;
    end
end
end

function B = alg4( B, R, alpha )
K = numel(B);
% forward pass
for k=2:K
    C = B(k-1) + alpha * (R(k) - R(k-1));
    if C < B(k)
        B(k) = C;
    end
end
% backward pass
for k=K-1:-1:1
    C = B(k+1) + alpha * (R(k+1) - R(k));
    if C < B(k)
        B(k) = C;
    end
end
end

[Michal]

Thanks for answer! These results significantly change the MATLAB programming best practice. Now is really hard to decide how to implement even simple algorithms.

Your best solution (alg4) eliminate all built-in functions...??!!

[Bruno Luong]

"Your best solution (alg4) eliminate all built-in functions"

10 year ago this is not the case.

The rule of thumbs now is for-loop with basic arithmetics is fast and competitive.