[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
unexpected behavior when multiplying large matrices
From: 
John W. Eaton 
Subject: 
unexpected behavior when multiplying large matrices 
Date: 
Sun, 1 Dec 2002 20:41:47 0600 
On 1Dec2002, E. Joshua Rigler <address@hidden> wrote:
 I have a rather complicated operation I want to perform with some large
 matrices. This is part of an adaptive optimization algorithm that I
 seed with the identity matrix, then update with each timestep. That
 really isn't important though.

 My question is why is there such a huge difference in the processing
 time required if, for example:

 octave:29> x = eye(1000);
 octave:30> y = [1:1000]' * [1:1000];
 octave:31> z = [1:1000]'
 octave:32> tic; a = 2 * (x  ((x * z * z' * x)/((1/2) + z' * x * z) ) ); toc
 ans = 0.59119
 octave:33> tic; b = 2 * (y  ((y * z * z' * y)/((1/2) + z' * y * z) ) ); toc
 ans = 42.064
 octave:34> whos
 ...
 *** local user variables:

 prot type rows cols name
 ==== ==== ==== ==== ====
 rwd matrix 1000 1000 a
 rwd matrix 1000 1000 b
 rwd matrix 1000 1000 x
 rwd matrix 1000 1000 y
 rwd matrix 1000 1 z

 octave:35>

 It almost seems like there must be some sort of 'sparse' optimization
 involved that knows that x is a diagonal matrix, and avoids the element
 by element multiplication when it isn't necessary. I'm using 2.1.35
 with all _but_ the most recent octaveforge package installed. Thanks.
Octave doesn't do any special optimization for identity matrices but
the Fortran version of level 3 blas routine dgemm does include some
checks for ones and zeros, and avoids multiplication and addition
operations where it can. Perhaps this sort of optimization accounts
for the timing difference that you see? I don't know whether this is
also true of the accelerated blas packages. It does not seem to be:
octave2.1:9> x = eye(1000);
octave2.1:10> y = [1:1000]' * [1:1000];
octave2.1:11> z = [1:1000]';
octave2.1:12> tic; a = 2 * (x  ((x * z * z' * x)/((1/2) + z' * x * z))); toc
ans = 2.1231
octave2.1:13> tic; b = 2 * (y  ((y * z * z' * y)/((1/2) + z' * y * z))); toc
ans = 2.0958
The above results were obtained on a 1.4GHz Athlon system running
Debian and Octave is linked with ATLAS:
$ ldd `type p octave`
liboctinterp.so => /usr/lib/octave2.1.40/liboctinterp.so (0x40014000)
liboctave.so => /usr/lib/octave2.1.40/liboctave.so (0x40431000)
libcruft.so => /usr/lib/octave2.1.40/libcruft.so (0x40640000)
liblapack.so.2 => /usr/lib/atlas/liblapack.so.2 (0x406a0000)
libblas.so.2 => /usr/lib/atlas/libblas.so.2 (0x40bcd000)
libfftw.so.2 => /usr/lib/libfftw.so.2 (0x40e6b000)
libreadline.so.4 => /lib/libreadline.so.4 (0x40e9f000)
libncurses.so.5 => /lib/libncurses.so.5 (0x40ec9000)
libdl.so.2 => /lib/libdl.so.2 (0x40f07000)
libm.so.6 => /lib/libm.so.6 (0x40f0a000)
libstdc++libc6.22.so.3 => /usr/lib/libstdc++libc6.22.so.3 (0x40f2b000)
libc.so.6 => /lib/libc.so.6 (0x40f74000)
/lib/ldlinux.so.2 => /lib/ldlinux.so.2 (0x40000000)
jwe

Octave is freely available under the terms of the GNU GPL.
Octave's home on the web: http://www.octave.org
How to fund new projects: http://www.octave.org/funding.html
Subscription information: http://www.octave.org/archive.html
