CODEWITHSUNDEEP
performancematlabjulia
user25004
user25004

Reputation: 2088

Strange observation about timing comparison between Julia and Matlab

The goal of this experiment is to compare speed of Matlab and Julia with a small piece of code below.

First the Matlab code:

>> t = 5000; n = 10000;  x = 1:t*n;
>> x = reshape(x, t, n);
>> tic(); y1 = sum(x(:) .* x(:)); toc()
Elapsed time is 0.229563 seconds.
>> y1

    y1 =
         4.1667e+22
>> tic(); y2 = trace(x * x'); toc()
Elapsed time is 15.332694 seconds.

>> y2
y2 =
       4.1667e+22

Versus in Julia

julia> t = 5000; n = 10000; x = 1: t*n;    
julia> x = reshape(x, t, n);
julia> tic(); y1 = sum(x[:].* x[:]); toc();
elapsed time: 1.235170533 seconds
julia> y1
-4526945843202100544
julia> tic();y2 = trace(x*x'); toc();

The second one did not finish the job in more than 1 minutes. So what is the matter here with Julia? This piece of code happen to run both slower and overflow in Julia? Is there any problem in my style? I think one reason to convert from Matlab to Julia is the speed, and I used to think that Julia handles big number arithmetic operations by default. Now, it looks like these are not correct. Can someone explain?

Upvotes: 2

Views: 406

Answers (1)

Simon Byrne
Simon Byrne

Reputation: 7874

There are a couple of things going on here.

Firstly, unlike Matlab, the x is an array of machine integers, not floating point values. This appears to be the main difference in speed, as it is unable to use BLAS routines for linear algebra.

You need to do either

x = 1.0:t*n

or explicitly convert it via

x = float(x)

This is also the reason for a different answer: Julia uses machine arithmetic, which for integers will wrap around on overflow (hence the negative number). You won't have this problem with floating point values.

(Julia does have arbitrary-precision integers, but doesn't promote by default: instead you would need to promote them yourself via big(n))

This will give you some speed up, but you can do better. Julia does require a slightly different programming style than Matlab. See the links that Isaiah provided in the comments above.

In regards to your specific examples

sum(x[:].* x[:])

This is slow as it creates 3 intermediate vectors (2 copies of x[:], though hopefully this will change in future, and their product).

Similarly, 

trace(x*x')

computes an intermediate matrix (most of which is not used in the trace).

My suggestion would be to use

dot(vec(x),vec(x))

vec(x) just reshapes x into a vector (so no copying), and dot is the usual sum-product. You could also "roll-your-own":

function test(x)
    s = zero(eltype(x)) # this prevents type-instability
    for xi in x
        s += xi*xi
    end
    s
end

test(x)

(this needs to be in a function for the JIT compiler to work its magic). It should be reasonably fast, though probably still not as fast as dot, which uses BLAS calls underneath.

Upvotes: 8

Related Questions