How to use casadi.dulmageMendelsohn for debugging purposes?

[Adrian Bürger]

Hello everyone,

first about the background: I have problems inverting a matrix using casadi.solve, since the matrix appears to be singular according to casadi.sprank. But from what I know, this matrix shouldn't be singular in principle, so I think I have a bug somewhere in my code which I wasn't yet able to resolve.

I got to know that for determination of the matrix rank with casadi.sprank, you use casadi.dulmageMendelsohn which is described to be an implementation of the Dulmage-Mendelsohn matrix decomposition algorithm, as I got to know e. g. from https://github.com/casadi/casadi/issues/601. So in principle, one should be able to - in addition to the matrix rank - get also information about the problematic (i. e. over- and under-determined) matrix parts from this function, which is something that I think would be of great use for debugging purposes (not only for this application, probably).

Within the issue mentioned above, the documentation for the Matlab dmperm function is referenced which shows i. a. the structure of the decomposed matrix and the several blocks contained. But from what I see, the depiction in the Matlab documentation shows a upper block diagonal matrix, while the CasADi implementation generates a lower block triangular matrix. So my questions would be:

- Is there yet some documentation about the structuring of the results returned by casadi.dulmageMendelsohn, apart for the example at https://github.com/casadi/casadi/blob/master/docs/api/examples/matrix/dulmageMendelsohn.py, that shows more clearly about the meaning of the returned block indices? Or can it directly be related to the according Matlab documentation (also for the slicing part, which I suppose to probably be different)?

- If I use this method for examination of matrices constructed from MX symbolics (which is what I would do), the output of slicing these matrices would probably be rather cryptic. What I thought of was using the casadi.Sparsity.spy function, which would probably yield already good information about the structure. In addition this this method, are there other possibilities or tricks to improve the clearness of the returned values of MX slicings?

Thank all of you really much for your help! And also, if there are other (and better) ways to address this problems, any hints are of course highly appreciated!

Best, Adrian

[Joel Andersson]

Hello Adrian,

The "casadi.dulmageMendelsohn" function, which has been renamed "casadi.btf" in CasADi 3.0, just returns the block-triangular form. I.e. a reordering of the rows and columns that would allow solving a sequence of smaller linear systems instead of one large linear system. CasADi's implementation is based on the one in CSparse, and if you want to understand what it does, you should have a look at Tim Davis' book "Direct Methods for Sparse Linear Systems".

About the slicings for MX, it's not really descriptive now, right. At some point, the slicing nodes in MX will need to be refactored to scale better with dimensions, but it has not been on the top of my list of priorities. Anyway, such a refactoring would also improve the readability.

Joel

[Adrian Bürger]

Hello Joel,

thank you for your answer and explanations, I will now have a look at the literature you suggested me!

Best, Adrian