[Developers] Fwd: [New post] Julia, Matlab, and C
John Sibert
sibert at hawaii.edu
Mon Sep 17 10:18:38 PDT 2012
-------- Original Message --------
Subject: [New post] Julia, Matlab, and C
Date: Mon, 17 Sep 2012 02:37:37 +0000
From: Justin Domkes Weblog <comment-reply at wordpress.com>
Reply-To: Justin Domkes Weblog
<comment+_dvt62oavm4btqopdkeujz at comment.wordpress.com>
To: sibert at hawaii.edu
WordPress.com
justindomke posted: "Julia is a new language in the same arena as Matlab
or R. I've had failed attempts to quit the Matlab addiction in the past,
making me generally quite conservative about new platforms. However,
I've recently been particularly annoyed by Matlab's slow spee"
Respond to this post by replying above this line
New post on *Justin Domke's Weblog*
<http://justindomke.wordpress.com/author/justindomke/>
Julia, Matlab, and C
<http://justindomke.wordpress.com/2012/09/17/julia-matlab-and-c/>
by justindomke <http://justindomke.wordpress.com/author/justindomke/>
Julia <http://julialang.org/> is a new language in the same arena as
Matlab or R. I've had failed attempts
<http://justindomke.wordpress.com/2008/11/29/mandelbrot-in-scala/> to
quit the Matlab addiction in the past, making me generally quite
conservative about new platforms. However, I've recently been
particularly annoyed by Matlab's slow speed, evil license manager
errors, restrictions on parallel processes, C++ .mex file pain, etc.,
and so I decided to check it out. It seems inevitable that Matlab will
eventually displaced by /something/. The question is: is that something
Julia?
"We want a language that’s open source, with a liberal license. We
want the speed of C with the dynamism of Ruby. We want a language
that’s homoiconic, with true macros like Lisp, but with obvious,
familiar mathematical notation like Matlab. We want something as
usable for general programming as Python, as easy for statistics as
R, as natural for string processing as Perl, as powerful for linear
algebra as Matlab, as good at gluing programs together as the shell.
Something that is dirt simple to learn, yet keeps the most serious
hackers happy. We want it interactive and we want it compiled."
Essentially, the goal seems to be a faster, freer Matlab that treats
users like adults (macros!) and /doesn't require writing any .mex files
in C++ or Fortan/. Sounds too good to be true? I decided to try it out.
My comparisons are to Matlab (R2012a) and C (gcc 4.2 with -O2 and
-fno-builtin to prevent compile-time computations) all on a recent
MacBook Pro (2.3 GHz Intel Core i7 w/ 8GB RAM).
Installation was trivial: I just grabbed a pre-compiled binary
<https://github.com/JuliaLang/julia/downloads> and started it up.
I deliberately used naive algorithms, since I am just testing raw speed.
It should be a fair comparison, as long as the algorithm is constant.
Please let me know about any bugs, though.
Click here to skip to the results <#results>.
First benchmark: Fibonnaci
% Matlab
function f=fib(n)
if n <= 2
f=1.0;
else
f=fib(n-1)+fib(n-2);
end
end
// C
double fib(int n){
if(n<=2)
return(1.0);
else
return(fib(n-2)+fib(n-1));
}
% julia
function fib(n)
if n <= 2
1.0
else
fib(n-1)+fib(n-2);
end
end
Clarity is basically a tie. Running them for n=30 we get:
time in matlab (fib): 14.344231
time in c (fib): 0.005795
time in julia (fib): 0.281221
Second benchmark: Matrix Multiplication
This is a test of the naive O(N^3) matrix multiplication algorithm.
% matlab
function C=mmult(A,B,C)
[M,N] = size(A);
for i=1:M
for j=1:M
for k=1:M
C(i,j) = C(i,j) + A(i,k)*B(k,j);
end
end
end
end
// C
#define M 500
void mmult(double A[M][M],double B[M][M], double C[M][M]){
//double C[M][M];
int i,j,k;
for(i=0; i<M; i++)
for(j=0; j<M; j++){
C[i][j] = 0;
for(k=0; k<M; k++)
C[i][j] += A[i][k]*B[k][j];
}
}
# julia
function mmult(A,B)
(M,N) = size(A);
C = zeros(M,M);
for i=1:M
for j=1:M
for k=1:M
C[i,j] += A[i,k]*B[k,j];
end
end
end
C;
end
Here, I think that Matlab and Julia and a bit clearer, and Julia wins
though the wonders of having "+=". The timing results on 500x500
matrices are:
time in matlab (matmult): 1.229571
time in c (matmult): 0.168296
time in julia (matmult): 0.505222
Third Benchmark: numerical quadrature
Here, we attempt to calculate the integral \int_{x=5}^{15} sin(x) dx by
numerical quadrature, using a simple midpoint rule with computations at
10^7 points.
% matlab
function val=numquad(lb,ub,npoints)
val = 0.0;
for x=lb:(ub-lb)/npoints:ub
val = val + sin(x)/npoints;
end
end
// C
double numquad(double lb,double ub,int npoints){
double val = 0.0;
int i;
for(i=0; i<=npoints; i++){
double x = lb + (ub-lb)*i/npoints;
val += sin(x)/npoints;
}
return(val);
}
# julia
function numquad(lb,ub,npoints)
val = 0.0
for x=lb:(ub-lb)/npoints:ub
val += sin(x)/npoints
end
val
end
The timings are:
time in matlab (numquad): 0.446151
time in c (numquad): 0.201071
time in julia (numquad): 0.273863
Fourth Benchmark: Belief Propagation
Finally, I decided to try a little algorithm similar to what I actually
tend to implement for my research. Roughly speaking, Belief Propagation
<http://en.wikipedia.org/wiki/Belief_propagation> is a repeated sequence
of matrix multiplications, followed by normalization.
% matlab
function x=beliefprop(A,x,N)
for i=1:N
x = A*x;
x = x/sum(x);
end
end
// C
void beliefprop(double A[25][25], double x[25], int N){
int i, n, j;
double x2[25];
for(n=0; n<N; n++){
for(i=0; i<25; i++){
x2[i]=0;
for(j=0; j<25; j++)
x2[i] += A[i][j]*x[j];
}
for(i=0; i<25; i++)
x[i]=x2[i];
double mysum = 0;
for(i=0; i<25; i++)
mysum += x[i];
for(i=0; i<25; i++)
x[i] /= mysum;
}
return;
}
% julia
function beliefprop(A,x,N)
for i=1:N
x = A*exp(x);
x /= sum(x);
end
x
end
Here, I think we can agree that Matlab and Julia are clearer. (Please
don't make fun of me for hardcoding the 25 dimensions in C.) Using a
matrix package for C would probably improve clarity, but perhaps also
slow things down. The results are:
time in matlab (beliefprop): 0.627478
time in c (beliefprop): 0.073564
time in julia (beliefprop): 0.489665
Fifth Benchmark: BP in log-space
In practice, Belief Propagation is often implemented in log-space (to
help avoid numerical under/over-flow.). To simulate an algorithm like
this, I tried changing to propagation to take an exponent before
multiplication, and a logarithm before storage.
% matlab
function x=beliefprop2(A,x,N)
for i=1:N
x = log(A*exp(x));
x = x - log(sum(exp(x)));
end
end
// C
void beliefprop2(double A[25][25], double x[25], int N){
int i, n, j;
double x2[25];
for(n=0; n<N; n++){
for(i=0; i<25; i++){
x2[i]=0;
for(j=0; j<25; j++)
x2[i] += A[i][j]*exp(x[j]);
}
for(i=0; i<25; i++)
x[i]=log(x2[i]);
double mysum = 0;
for(i=0; i<25; i++)
mysum += exp(x[i]);
double mynorm = log(mysum);
for(i=0; i<25; i++)
x[i] -= mynorm;
}
return;
}
# julia
function beliefprop2(A,x,N)
for i=1:N
x = log(A*exp(x));
x -= log(sum(exp(x)));
end
x
end
Life is too short to write C code like that when not necessary. But how
about the speed, you ask?
time in matlab (beliefprop2): 0.662761
time in c (beliefprop2): 0.646153
time in julia (beliefprop2): 0.615113
Sixth Benchmark: Markov Chain Monte Carlo
Here, I implement a simple Metropolis algorithm. For no particular
reason, I use the two-dimensional distribution:
p(x) \propto \exp(\sin(5 x_1) - x_1^2 - x_2^2)
<http://justindomke.files.wordpress.com/2012/09/dist.jpg>
% matlab
function mcmc(x,N)
f = @(x) exp(sin(x(1)*5) - x(1)^2 - x(2)^2);
p = f(x);
for n=1:N
x2 = x + .01*randn(size(x));
p2 = f(x2);
if rand < p2/p
x = x2;
p = p2;
end
end
end
// C
double f(double *x){
return exp(sin(x[0]*5) + x[1]*x[1]);
}
#define pi 3.141592653589793
void mcmc(double *x,int N){
double p = f(x);
int n;
double x2[2];
for(n=0; n<N; n++){
// run Box_Muller to get 2 normal random variables
double U1 = ((double)rand())/RAND_MAX;
double U2 = ((double)rand())/RAND_MAX;
double R1 = sqrt(-2*log(U1))*cos(2*pi*U2);
double R2 = sqrt(-2*log(U1))*sin(2*pi*U2);
x2[0] = x[0] + .01*R1;
x2[1] = x[1] + .01*R2;
double p2 = f(x2);
if(((double)rand())/RAND_MAX< p2/p){
x[0] = x2[0];
x[1] = x2[1];
p = p2;
}
}
}
% julia
function mcmc(x,N)
f(x) = exp(sin(x[1]*5) - x[1]^2 - x[2]^2);
p = f(x);
for n=1:N
x2 = x + .01*randn(size(x));
p2 = f(x2);
if rand() < p2/p
x = x2;
p = p2;
end
end
end
Again, I think that C is far less clear than Matlab or Julia. The
timings are:
time in matlab (mcmc): 7.747716
time in c (mcmc): 0.646153
time in julia (mcmc): 0.613223
Table
Matlab C Julia
fib 14.344 0.005 0.281
matmult 1.229 0.168 0.505
numquad 0.446 0.201 0.273
bp 0.627 0.073 0.489
bp2 0.662 0.646 0.615
mcmc 7.747 0.646 0.613
Conclusions
I'm sure all these programs can be sped up. In particular, I'd bet that
an expert could optimize the C code to beat Julia on |bp2| and |mcmc|.
These are a test of "how fast can Justin Domke make these programs", not
the intrinsic capabilities of the languages. That said, Julia allows for
optional type declarations. I did experiment with these but found
absolutely no speed improvement. (Which is a good or a bad thing,
depending on how you look at life.)
Another surprise to me was how often Matlab's JIT managed a speed within
a reasonable factor of C. (Except when it didn't...)
The main thing that at Matlab programmer will miss in Julia is
undoubtedly plotting. The Julia designers seem to understand the
importance of this ("non-negotiable"). If Julia equalled Matlab's
plotting facilities, Matlab would be in real trouble!
Overall, I think that the killer features of freedom, kinda-sorta-C-like
speed, and ease of use make Julia more likely as a Matlab-killer than
other projects such as R, Sage, Octave, Scipy, etc. (Not to say that
those projects have not succeeded in other ways!) Though Julia's
designers also seem to be targeting current R users, my guess is that
they will have more success with Matlab folks in the short term, since
most Matlab functionality (other than plotting) already exists, while
reproducing R's statistical libraries will be quite difficult. I also
think that Julia would be very attractive to current users of languages
like Lush <http://lush.sourceforge.net/>. Just to never write another
.mex file, I'll very seriously consider Julia for new projects. Other
benefits such as macros, better parallelism support are just bonuses. As
Julia continues to develop, it will become yet more attractive.
There was an interesting discussion on Lambda the Ultimate
<http://lambda-the-ultimate.org/node/4452> about Julia back when it was
announced
*justindomke <http://justindomke.wordpress.com/author/justindomke/>* |
September 17, 2012 at 2:37 am | Tags: c++
<http://justindomke.wordpress.com/?tag=c>, efficiency
<http://justindomke.wordpress.com/?tag=efficiency>, julia
<http://justindomke.wordpress.com/?tag=julia>, matlab
<http://justindomke.wordpress.com/?tag=matlab> | Categories:
Uncategorized <http://justindomke.wordpress.com/?cat=1> | URL:
http://wp.me/pgOXE-j8
Comment
<http://justindomke.wordpress.com/2012/09/17/julia-matlab-and-c/#respond>
See all comments
<http://justindomke.wordpress.com/2012/09/17/julia-matlab-and-c/#comments>
Unsubscribe or change your email settings at Manage Subscriptions
<https://subscribe.wordpress.com/?key=4220e3bc2f2ce5bd9d8f568e895ce9d4&email=sibert%40hawaii.edu>.
*Trouble clicking?* Copy and paste this URL into your browser:
http://justindomke.wordpress.com/2012/09/17/julia-matlab-and-c/
Thanks for flying with WordPress.com <http://wordpress.com>
More information about the Developers
mailing list