ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Mon, 28 Oct 2013 18:52:08 +0100Speed comparison with numerical arrays: sage, R, matlabhttps://ask.sagemath.org/question/10671/speed-comparison-with-numerical-arrays-sage-r-matlab/Dear all
I am a sage newbie coming from matlab, R, and Maple, and although I really like many features, I find that my sage code for handling numerical arrays is much more complex and much slower than the corresponding code in R and Matlab.
Here's an example:
# Simuating a sample path of standard Brownian motion
I want to simulate a random path of standard Brownian motion. The input is a vector of time points, and I want returned a (random) trajectory, sampled at these time points.
To do this I need: 1) Fast random number generation 2) Fast array operations: Diff, multiplication, and cumsum. I found that in sage I need to go through NumPy, so I wrote:
sage: import numpy
sage: # General routine to sample a Brownian motion on a time grid
sage: def rB(tvec):
sage: dt = numpy.diff(numpy.array(tvec))
sage: sdt = sqrt(dt)
sage: RNG = RealDistribution('gaussian',1)
sage: dB = sdt * numpy.array([RNG.get_random_element() for j in range(dt.__len__())],dtype=float)
sage: B1 = numpy.array([normalvariate(0,sqrt(tvec[0]))],dtype=float)
sage: B = numpy.cumsum(numpy.concatenate((B1,dB)))
sage: return(B)
I then time it, testing with 1 million time points:
sage: tvec = srange(0,1e3,1e-3)
sage: time B = rB(tvec)
I get a time of 3.5 seconds on my machine, a standard Laptop running Ubuntu, running sage from the shell. Using prun, I see that half of the time, roughly, is spent generating random numbers and the other half mainly constructing numerical arrays.
For comparison, the R code would be (here I assume that tvec[1]=0)
> tvec <- seq(0,1000,1e-3)
> system.time(B <- c(0,cumsum(rnorm(length(tvec)-1,sd=sqrt(diff(tvec))))))
which executes in 0.24 seconds. The matlab code would be
>> tvec = linspace(0,1e3,1e6+1);
>> tic;B = [0,cumsum(randn(1,length(tvec)-1) .* sqrt(diff(tvec)))];toc
which runs in 0.06 seconds.
You see that my sage code is not only clumsier, it is also between 10 and 100 times slower, so I guess I am doing something wrong.
I would really appreciate if you could point me in the direction of cleaner and faster code.
Cheers,
Uffe
Uffe_H_ThygesenMon, 28 Oct 2013 18:52:08 +0100https://ask.sagemath.org/question/10671/