permutations of matrix to triangular form?

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

permutations of matrix to triangular form?

David Bateman-3
Thinking yet again about implementing the triangular/banded solvers,
I re-examined what matlab did as described at

http://www.mathworks.com/access/helpdesk/help/techdoc/ref/mldivide.html

Where it seems to test the type of the matrix (and presumably caches
it) before solving with the "\" and "/" operators. Presumably
something similar is done in a few other places like eig, etc.  Doing
it matlabs way has the advantage that the banded and triangular types
are not explicitly needed.

Last night I was investigating how to do something similar with the
sparse class, and basically just wrote a small piece of code to recognize
the type of a matrix (I attach it here if you are interested, though it
needs to be used with the sparse class I've been writing). The interest
of this was to take Andy's back substitution code from his sparse inverse
and use it directly in the "\" and "/" operations, and at the same time
treat banded with LAPACK and diagonal matrices.

In any case, it seems that I can recognize all of the matrix types I
programmed for. But the permuted triangular form is not so clear-cut.
The permutation vector is not explicitly passed and so must be
reconstructed. As this is a test to check which solver will be used,
the means of creating this permutation, or at least checking that it
is possible must be low cost. Frankly, I can't see what matlab is doing.

What I think might be happening is that they are performing an LU
factorization of the matrix and is either L or U of the factorization
is diagonal, then only a single back substitution is needed.  Though
in this case I don't see why matlab makes a distinction between case 6
and case 7 in the docs for mldivide, as there won't be all that much
difference in the computation time.

Does anyone know of a simple, low-cost means of attempting to convert
find where a permutation of a matrix to triangular form exists? Does
anyone have any more understanding of what matlab is doing in this
case?

Regards
David

--
David Bateman                                [hidden email]
Motorola CRM                                 +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax)
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as:

[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary


Reply | Threaded
Open this post in threaded view
|

Re: permutations of matrix to triangular form?

Andy Adler
On Wed, 8 Dec 2004, David Bateman wrote:

> Thinking yet again about implementing the triangular/banded solvers,
> I re-examined what matlab did as described at
>
> http://www.mathworks.com/access/helpdesk/help/techdoc/ref/mldivide.html
>
> Where it seems to test the type of the matrix (and presumably caches
> it) before solving with the "\" and "/" operators. Presumably
> something similar is done in a few other places like eig, etc.  Doing
> it matlabs way has the advantage that the banded and triangular types
> are not explicitly needed.
>
> Last night I was investigating how to do something similar with the
> sparse class, and basically just wrote a small piece of code to recognize
> the type of a matrix (I attach it here if you are interested, though it
> needs to be used with the sparse class I've been writing). The interest
> of this was to take Andy's back substitution code from his sparse inverse
> and use it directly in the "\" and "/" operations, and at the same time
> treat banded with LAPACK and diagonal matrices.

David,

This doesn't answer your question, but I think that this approach
(of having a small number of storage types (full,sparse) with 'hints'
about the solver) offers a nice set of possibilities.

Currently, your approach autodetects the matrix type, but these
properties could be made available to the user, such as

   a= sparse(diag(1:3));
   get(a, 'storage-type')
   a= set(a, 'storage-type', 'diagonal')

or some similar syntax. This would allow

1. Users who already know something about the matrix to write more
   efficient code

2. The possibilitiy of including more detailed information to the
   solver. (ie. iterative, precision thresholds )

   a= set(a, 'solver', 'iterative');

   Currently, doing this in Matlab with sparse code makes it look very
   unlike the simple solver cases. For example, in some finite element
   code I'm working on (http://cvs.sf.net/viewcvs.py/eidors3d/eidors3d/algorithms/np_2003/forward_solver.m)
   instead of
        V= E\I;
   one does
        %Permute the rows and columns to make the factors sparser
        E = E(pp,pp);
        In = I(pp,:);
        rr(pp)=1:max(size(pp));
        U = cholinc(E,tol/10);
        q_c =  U' \ In;
        Vn = U \ q_c;
        %De-permute the result for Cholesky
        V = Vn(rr,:);

I know you've suggested a similar idea before. I'm just pointing out that
it works with the framework you're working on now.

Another issue: There is the possibility of user's giving the wrong info
about the matrix.  My opinion would be that that is their problem.

--
Andy Adler <adler AT site.uOttawa.ca> 1(613)562-5800x6218


Reply | Threaded
Open this post in threaded view
|

Re: permutations of matrix to triangular form?

David Bateman-3
According to Andy Adler <[hidden email]> (on 12/08/04):
> David,
>
> This doesn't answer your question, but I think that this approach
> (of having a small number of storage types (full,sparse) with 'hints'
> about the solver) offers a nice set of possibilities.

Not my idea yet, as all I've done is implement detecting of the matrix
type as per the matlab mldivide spec and in the same precedence. Including
that in the sparse solvers themselves is the next step

> Currently, your approach autodetects the matrix type, but these
> properties could be made available to the user, such as
>
>    a= sparse(diag(1:3));
>    get(a, 'storage-type')
>    a= set(a, 'storage-type', 'diagonal')

The matlab way is that diagonal and banded matrices are only available
in the sparse type, as why would you store them any other way. Finding
a diagonal, banded or unpermuted triangular matrix is a releatively
low cost operation, and so allowing the user to specify the type just
introduces the possiblity that they'll get it wrong without a huge
speed gain. I'd hesitate to do that.

> 2. The possibilitiy of including more detailed information to the
>    solver. (ie. iterative, precision thresholds )
>
>    a= set(a, 'solver', 'iterative');

This is what John was against in a previous thread as the setting the
solver type can be disassociated from the operation itself and increases
the potential for confusion.

>    Currently, doing this in Matlab with sparse code makes it look very
>    unlike the simple solver cases. For example, in some finite element
>    code I'm working on (http://cvs.sf.net/viewcvs.py/eidors3d/eidors3d/algorithms/np_2003/forward_solver.m)
>    instead of
>         V= E\I;
>    one does
>         %Permute the rows and columns to make the factors sparser
>         E = E(pp,pp);
>         In = I(pp,:);
>         rr(pp)=1:max(size(pp));
>         U = cholinc(E,tol/10);
>         q_c =  U' \ In;
>         Vn = U \ q_c;
>         %De-permute the result for Cholesky
>         V = Vn(rr,:);

I agree this is ugly, but it is safer. Perhaps the way around this is
to introduce mldivide/mrdivide functions (we'll have to do it anyway for
matlab style "@" classes), but with an additional argument that can force
a particular solver to be used. In this manner, we remain compatiable,
don't disassociate the command for the solver type from the operator
itself and make the above simplier. Which is simplier


    E = set(E, 'solver', 'cholesky');
    V= E\I;

or

    V = mldivide (E, I, 'cholesky')

This of course means that "mldivide" would have to be explicitly used
in such cases.

> I know you've suggested a similar idea before. I'm just pointing out that
> it works with the framework you're working on now.
>
> Another issue: There is the possibility of user's giving the wrong info
> about the matrix.  My opinion would be that that is their problem.

Its too easy in your scheme to have the setting of the solver far away
from its use, with potential to introduce horrible bugs that are difficult
to find. One debugging aid in that case is the spparms('spmoni') flag
that can be used to detect the solver type used. Then the problem becomes
one of educating users. I what a simple life, so I feel that its probably
still bad to have such a scheme.

Cheers
David

--
David Bateman                                [hidden email]
Motorola CRM                                 +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax)
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as:

[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary


Reply | Threaded
Open this post in threaded view
|

Re: permutations of matrix to triangular form?

Andy Adler
On Wed, 8 Dec 2004, David Bateman wrote:

> >    Currently, doing this in Matlab with sparse code makes it look very
> >    unlike the simple solver cases. For example, in some finite element
> >    code I'm working on (http://cvs.sf.net/viewcvs.py/eidors3d/eidors3d/algorithms/np_2003/forward_solver.m)
> >    instead of
> >         V= E\I;
> >    one does
> >         %Permute the rows and columns to make the factors sparser
> >         E = E(pp,pp);
> >         In = I(pp,:);
> >         rr(pp)=1:max(size(pp));
> >         U = cholinc(E,tol/10);
> >         q_c =  U' \ In;
> >         Vn = U \ q_c;
> >         %De-permute the result for Cholesky
> >         V = Vn(rr,:);
>
> I agree this is ugly, but it is safer. Perhaps the way around this is
> to introduce mldivide/mrdivide functions (we'll have to do it anyway for
> matlab style "@" classes), but with an additional argument that can force
> a particular solver to be used. In this manner, we remain compatiable,
> don't disassociate the command for the solver type from the operator
> itself and make the above simplier. Which is simplier
>
>
>     E = set(E, 'solver', 'cholesky');
>     V= E\I;

This is basically what I was suggesting.

While I know that this isn't the point of what you're currently
thinking of, it seems to me that the ability to add this kind
of information to matrices is powerful.

For example, one thing I would love to be able to do, is have a matrix
object which would allow arbitrary pieces of cached computations to be
tacked onto it. My current code tends to keep this information in
global variables - not the cleanest solution.

I realize I haven't carefully looked at your new sparse code yet.
I hope to do that this week.

Thanks,
--
Andy Adler <adler AT site.uOttawa.ca> 1(613)562-5800x6218



Reply | Threaded
Open this post in threaded view
|

Re: permutations of matrix to triangular form?

David Bateman-3
In reply to this post by David Bateman-3

> >
> >     E = set(E, 'solver', 'cholesky');
> >     V= E\I;
>
> This is basically what I was suggesting.

Yes, but my question was, isn't

    V = mldivide (E, I, 'cholesky')

simpler, as you can't seperate the set of the solver type from the solve
itself....

> While I know that this isn't the point of what you're currently
> thinking of, it seems to me that the ability to add this kind
> of information to matrices is powerful.
>
> For example, one thing I would love to be able to do, is have a matrix
> object which would allow arbitrary pieces of cached computations to be
> tacked onto it. My current code tends to keep this information in
> global variables - not the cleanest solution.

It seems to me that a better way of attacking this would be to have
classes and define your type as a sparse, with a string to select
the solver type. You can then overload the operators \ and / with
some selection code in it. With that capability, you'd have what
you want..

> I realize I haven't carefully looked at your new sparse code yet.
> I hope to do that this week.

I'm not going to touch this code in the next few weeks as I'm on
vacation. So it won't change in the next few weeks... However, I
checked-in last night most of the polymorphic solver code, including
diagonal, permuted diagonal, upper/lower triangular and banded
solvers. Its missing permuted upper/lower triangular, and cholesky
solvers still, as well as the QR solvers that were always missing.
I haven't benchmarked the code yet, but I've included test cases that
fully exercises this code. This should be capable of replacing Paul's
tri-diagonal and the SymBand stuff in octave-forge in the long run.
In any case, it would be interesting to have some benchmarking of the
special casing of the solver to see how it compares... The
spparms("bandden") option can be used to change between the banded
and the UMFPACK solvers if you want to look at this..

Regards
David

--
David Bateman                                [hidden email]
Motorola CRM                                 +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax)
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as:

[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary