Re: sparse matrix operations

[Miroslaw Kwasniak]

On Tue, Nov 16, 2004 at 11:32:28AM -0700, E. Joshua Rigler wrote:
> For example, suppose I only want the diagonal of the product of two
> square matrices.  As it stands right now in Octave (with octave-forge),
> C=A*B will result in a full C matrix, and may potentially use up more
> memory than I wish to allow. 

Yes

octave:20> a=[1 3 1;0 0 0;1 1 2],b=a',c=b*a,
a =
  1  3  1
  0  0  0
  1  1  2
b =
  1  0  1
  3  0  1
  1  0  2
c =
   2   4   3
   4  10   5
   3   5   5

> Since I only want the diagonal of C, I may
> as well only perform those operations that place a vector product at
> C(i,i), thereby saving CPU cycles and memory.
> 
> I have a simple .m file that does this, but of course it is very slow in
> the Octave interpreter since it must loop over each of the desired
> non-zero elements of C.  I was hoping somebody out there might have
> something similar, but in a .oct format, or could show me a trick for
> doing this with standard Octave operators, before I go and try to
> reinvent this wheel.

I don't know how smart is sparse implementation in octave but you can try

octave:21> sum(a'.*b)
ans =

   2  10   5

In which c=a'.*b should't use more memory then smallest of one of a,b. It
doesn't save so mach memory as in place operation but save more cpu the
m-file loop.

My first test gives resonable results:

octave:49> sp_diag
N = 1000
mult sp=0 t=0.075325 e=0
sum  sp=0 t=0.046321 e=0
loop sp=0 t=1.667849 e=0
mult sp=1 t=0.116958 e=0
sum  sp=1 t=0.157900 e=0
loop sp=1 t=5.076967 e=0

sum metod is better then mult on full matrix, but on sparse execution time
isn't to bad :)

Mirek


file: sp_diag.m
-----------------------------------
N=1000
clear a0 a b c c0;

a0=[1 0 7;0 0 3;1 0 2];
a=a0;
a=[a a a a a a a a a a]';
a=[a a a a a a a a a a];

L=size(a,1);

for sp=0:1
  if sp; a=sparse(a); end
  b=a';

  tic
  for n=1:N;
    c=diag(a*b);
  end
  printf('mult sp=%d t=%f e=%d\n',sp,toc,0);
  c0=c;

  tic
  for n=1:N;
    c=sum(a'.*b);
  end
  printf('sum  sp=%d t=%f e=%d\n',sp,toc,any(c'~=c0));

  tic
  for n=1:N;
    c=zeros(L,1);
    for l=1:L
      c(l)=sum(a(l,:).*b(:,l)');
    end
  end
  printf('loop sp=%d t=%f e=%d\n',sp,toc,any(c~=c0));
end
-----------------------------------



-------------------------------------------------------------
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
-------------------------------------------------------------