unexpected behavior when multiplying large matrices

[John W. Eaton]

On  1-Dec-2002, 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 time-step.  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 octave-forge 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/octave-2.1.40/liboctinterp.so (0x40014000)
    liboctave.so => /usr/lib/octave-2.1.40/liboctave.so (0x40431000)
    libcruft.so => /usr/lib/octave-2.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.2-2.so.3 => /usr/lib/libstdc++-libc6.2-2.so.3 (0x40f2b000)
    libc.so.6 => /lib/libc.so.6 (0x40f74000)
    /lib/ld-linux.so.2 => /lib/ld-linux.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
-------------------------------------------------------------