RE: sparce matrices

[adler]

On Sat, 13 Apr 2002, Mike Miller <address@hidden> wrote:
>It would be great if I could specify the non-zero elements of a big block-
>diagonal matrix, call it X, say, then do things like chol(X)*Y without
>having to allocate a lot of memory.  I'm dealing with situations where I
>might have, say, 1,000 8x8 blocks.  That means there would be 64,000
>non-zero elements in a matrix with 64,000,000 elements.
>
>Can octave deal with this kind of thing efficiently?  I've been dealing
>with it by making 1,000 distinct 8x8 matrices, doing 1,000 Cholesky
>factorizations, multiplying the 1,000 resulting matrices by 1,000 pieces
>of Y and then assembling the products into one big matrix.

There are several solutions available to you in octave-forge.

If you have banded sparse matrices, then the
chol routines in octave-forge/extra/symband/ will do what you
want.

For general sparse matrices, you should use the stuff in
octave-forge/main/sparse. This currently does not support
chol - but then for arbitrary sparse matrices, chol is
typically a bad idea: 1) chol will fill out your matrix
out to the furthest from the diagonal nonzero, 2) chol
does not do the column reordering that is so important
for numerical stability for arbitrary large sparse matrices.

using splu, you can calculate
   [L,U,Pr,Pc] = splu(A)
       returns unit lower triangular L, upper triangular U,
         and permutation matrices Pr,Pc with Pr*A*Pc' = L*U.

If A is symmetric, then you can calculate a diagonal matrix
D such that
    A= P*L*D*L'*P'

and if D is non negative, then
    R= L*sqrt(D);
    A= P*R*R'*P';

For arbitrary large sparse matrices, this permuted chol
will most likely give better performance (speed and numerical).
(Note: I haven't actually tested any of this)

andy

_______________________________________
Andy Adler,                address@hidden





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