Sparse Toeplitz matrices

Marco Caliari marco.caliari at univr.it
Fri Apr 10 10:33:13 CDT 2009


Dear maintainers,

I usually use finite differences (sparse) matrices and I find quite 
convenient to use the function toeplitz.m to generate them. However, I 
discovered that

toeplitz(sparse([-2,1,zeros(1,N-2)]))

can be quite slow if N is big.
I wrote the attached function which should be used only to generate sparse 
Toeplitz matrices. Against Octave 3.1.54 I can see the following speedups:

octave-3.1.54:1> N = 5000;
octave-3.1.54:2> tic,A = toeplitz(sparse([-2,1,zeros(1,N-2)]));,toc
Elapsed time is 4.78363 seconds.
octave-3.1.54:3> tic,B = mytoeplitz(sparse([-2,1,zeros(1,N-2)]));,toc
Elapsed time is 0.095167 seconds.
octave-3.1.54:4> A-B
ans =

Compressed Column Sparse (rows = 5000, cols = 5000, nnz = 0 [0%])
octave-3.1.54:5> 
tic,A=toeplitz(sparse([0,-1,zeros(1,N-2)]),sparse([0,1,zeros(1,N-2)]));,toc
Elapsed time is 3.8 seconds.
octave-3.1.54:6> 
tic,B=mytoeplitz(sparse([0,-1,zeros(1,N-2)]),sparse([0,1,zeros(1,N-2)]));,toc
Elapsed time is 0.021 seconds.
octave-3.1.54:7> A-B
ans =

Compressed Column Sparse (rows = 5000, cols = 5000, nnz = 0 [0%])

If interested, I can try to make a patch against toeplitz.m.

Best regards,

Marco
-------------- next part --------------
function A = mytoeplitz(c,r)
if (nargin == 1)
  r = c;
  c = conj (c);
  c(1) = conj (c(1));
endif
if (r (1) != c (1))
  warning ("toeplitz: column wins diagonal conflict");
  r(1) = 0;
endif
n = length(c);
c = c(:).';
r = r(:).';
cidx = find(c);
ridx = find(r);
A = spdiags(repmat(full(c(cidx)),n,1),1-cidx,n,n)+...
spdiags(repmat(full(r(ridx)),n,1),ridx-1,n,n);


More information about the Octave-maintainers mailing list