benchmarks

[Paul Kienzle]

Paul Thomas pointed me to these benchmarks.  Anybody want to do
something about them?

        http://www.sciviews.org/other/benchmark.htm

Speed is relative to R.  Note that he didn't test on a variety of 
matrices,
so the result may represent random perturbations in the generated
matrices rather than any real difference between the systems.

============================
I.A   Matlab 0.48 : Octave 2.01

  tic;
    a = abs(randn(1500, 1500)/10);
    b = a';
    a = reshape(b, 750, 3000);
    b = a';
  timing = toc;

Octave's random number generators are slower than they need
to be.  octave-forge helps, but there are even better ones available.
E.g., http://www.jstatsoft.org/v05/i08/

I've written to the author asking for permission to release it under
an GPL compatible license, but didn't get a response.  I haven't
had time to recode it in my own words.  It'll be tricky since the paper
contains the source code.

I.B Matlab 0.42  - Octave 1.22
    b = a.^1000;

Looking at the source, Octave is using indexing for every element.
Try replacing the nested for loop and indexing with a call to
data and directly walking the array.  Bonus, it automatically works
for n-D arrays.

I.C Matlab 0.89  - Octave 7.77
    b = sort(a);

Octave's sort is surprisingly slow.  3x worse than any other package
mentioned.  Anyone know a fast stable sort algorithm?

I.D Matlab 0.73 - Octave 0.35 -  R 0.24
    b = a'*a;

Octave could be faster if it skipped the transpose and used blas
directly with the appropriate transpose code.   Note, the tester
did not use a self-generated atlas, so these results don't mean
much.

I.E Matlab 0.24 - Octave 0.78
    c = a\b';

===================================

I don't understand how O-Matrix is consistently faster
on the following tests.  Aren't they all using the same
underlying libraries?  Or is O-Matrix hand-optimized
for Intel architecture?

II.A Matlab 0.48 - Octave 0.96 - O-Matrix 0.17
    b = fft(a);
Not sure if he used fftw.

II.B Matlab 0.86 - Octave 2.30 - O-Matrix 0.44
    b = eig(a);

II.C Matlab 0.27 - Octave 1.02 - O-Matrix 0.17
   b = det(a);

II.D Matlab 0.33 - Octave 0.21 - O-Matrix 0.22
   b = chol(a);

II.E Matlab 0.23 - Octave 0.47 - O-Matrix 0.11
   b = inv(a);

========================================
III.A Matlab 2.11 - Octave 2.06 - O-Matrix 0.59
    b = (phi.^a - (-phi).^(-a)) / sqrt(5);

III.B Matlab 0.84 - Octave 0.73 - O-Matrix 0.47
    b = ones(a, a)./((1:a)' * ones(1, a) + ones(a, 1) * (0:(a-1)));

III.C Matlab 0.91 - Octave 0.42
    c = gcd2(a, b);                            % gcd2 is a recursive 
function

This is a surprise --- the function is short and presumably
compiled in matlab, yet octave performs as well as the best
of all the packages.

III.D Matlab 0.38 - Octave 4.39 - O-Matrix 0.08
    for j = 1:220
      for k = 1:220
        b(k,j) = abs(j - k) + 1;
      end
    end

It is because of this case and the next that I want to
spit out some Parrot or other VM bytecode from the
octave tree-walker.

III.E Matlab 1.92 - Octave 3.08 - O-Matrix 0.52
    % Calculation of Escoufier's equivalent vectors

Lots of lines of vector code.  Not as bad as I expected.

==================================

We can do something about the first batch pretty easily.

The second batch, I question the methodology (he
should have used many different random matrices
and averaged rather than using the same random
matrix multiple times).  There are hand-optimized
blas for pentium available:

        http://www.cs.utexas.edu/users/flame/goto/

but the license isn't so friendly.  Anyone tried it?

The third batch will be a lot of work --- may something
to work toward for Octave 3.2.

BTW, moving more stuff from octave-forge to octave,
and providing a CRAN-like list of available packages
will help make octave look more complete.  This has
to be linked directly from the octave home page so that
people can find it.

Paul Kienzle
address@hidden