[Haskell-cafe] matrix computations based on the GSL

[Alberto Ruiz]

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)