On Thursday 16 March 2006 18:13, Frederik Eaton wrote:
Also, in my experiments (with matrix inversion) it seems, subjectively, that Octave is about 5 or so times faster for operations on large matrices. Presumably you've tested this as well, do you have any comparison results?
Frederik, you are right, in my machine Octave is at least 3x faster than gsl (v1.5). Too much, specially in a simple matrix multiplication, and I don't know why. See below the times measured in the C side for the PCA example. I will look into this... Alberto gsl 1.5 GHC> main Loading package haskell98-1.0 ... linking ... done. GSL Wrapper submatrix: 0 s GSL Wrapper constant: 0 s GSL Wrapper multiplyR (gsl_blas_dgemm): 0 s GSL Wrapper vector_scale: 0 s GSL Wrapper vector_scale: 1 s GSL Wrapper gsl_vector_add: 0 s GSL Wrapper trans: 0 s GSL Wrapper multiplyR (gsl_blas_dgemm): 36 s <---- (784x5000) x (5000 x 784) GSL Wrapper vector_scale: 0 s GSL Wrapper trans: 0 s GSL Wrapper vector_scale: 0 s GSL Wrapper gsl_vector_add: 0 s GSL Wrapper toScalar: 0 s GSL Wrapper eigensystem: 11 s <---------- (784x784) [337829.65625394206,253504.7867502207, etc. (97.95 secs, 13096315340 bytes) (includes file loading, to be optimized) GNU Octave, version 2.1.57 (i686-pc-linux-gnu). octave:6> testmnist x - repmat(mean(x),rows(x),1) 0.46144 (xc'*xc)/rows(x) 11.623 <-------------- eig 4.3873 <-------------- 3.3776e+05 2.5345e+05 2.0975e+05 etc. ----------octave------------------ load mnist.txt x = mnist(:,1:784); d = mnist(:,785); t0=time(); xc = x - repmat(mean(x),rows(x),1); disp("x - repmat(mean(x),rows(x),1)"); disp(time()-t0) t0=time(); mc = (xc'*xc)/rows(x); disp("(xc'*xc)/rows(x)"); disp(time()-t0) t0=time(); [v,l]=eig(mc); disp("eig"); disp(time()-t0) disp(flipud(diag(l))(1:10)); ------------- GSL + Haskell --------- module Main where import GSL import GSL.Util mean x = sumCols x ./. fromIntegral (rows x) center x = x |-| fromRows (replicate (rows x) (mean x)) cov x = (trans xc <> xc) ./. n where n = fromIntegral (rows x -1) xc = center x m ./. r = m <> (1/r ::Double) test :: M -> [Double] test x = take 10 (toList lambdas) where mc = cov x (lambdas, _) = eig mc main = do m <- fromFile "mnist.txt" let x = subMatrix 0 (rows m-1) 0 (cols m -2) m print (test x)