Efficient inversion of sparse symbolic matrices

[Adrian Bürger]

Hello everyone,

the core of my problem can already be taken from the title: how can I efficiently compute inverse matrices of sparse symbolic matrices? I hope that in the following, I can explain the problem in an understandable way.

A summary of the computations I have to do can be seen here: https://gist.github.com/adbuerger/89ef1e4f47661ae7c1a4

Apart from the casadi.mul()-call for computing h1 that can take some time (within the scope of half a minute), the matrix inversion (at least, until now) seems to be the most time consuming element of the computations. If I get it correctly, I now basically have the following choices:

- I can either use the SX or the MX class, and also convert the previous computation (which is some parameter estimation using IPOPT with a NLP set up using the collocation method with Lagrange polynomials) to meet the requirements for this. When I use the MX class, the previous computations for my small test example take up to 130.4 s, while when using the SX class everything can be done in 0.4 (!) s, so using MX-symbolics I imagine not to finish in acceptable time (if ever :-). So, I think I'm bound to use the SX-class.

- I can either use casadi.inv() or casadi.solve() and a suitable identity matrix to compute the inverse, but I was told it's usually better to use casadi.solve() to do this. When I try to invert using casadi.solve(), I quickly run out of memory (because more than ~14 GB of memory are needed), so I can not make a statement on the computational time. Using casadi.inv() again takes very long, I cancelled the computation after several minutes because this would already be too long for my application.

Additionally I have to say that at least for me there seems to be no chance of avoiding the matrix inversion.

As you are more experienced with the internal function principles of CasADi, can you advise me what I can do to speed up these computations? If I can provide you additional information on my problem, please tell me!

Thank you all very much for your support!

Best regards,

Adrian

[Greg Horn]

You want to use casadi.solve(), but with a better sparse linear solver. The default linear solver is a convenience and isn't meant for large systems.

I'm actually not sure if you need to use LinearSolver class, or if you can do something like casadi.solve(A,b,solver=CSparse)

[Joel Andersson]

Hi!

You should try to use CasADi matrices instead of numpy matrices. numpy.zeros, numpy.eye give dense matrices. SX.eye, MX(3,4), etc. are sparse matrices. If you want to use "standard" datatype, you could use e.g. scipy.sparse.csc_matrix, but best is probably to work with CasADi matrices directly. Note that SX(2,3) is an all-zero sparse 2-by-3 matrix, SX.zeros(2,3) is a dense 2-by-3 matrix. This should make a huge difference.

Other things: Calculating a Jacobian and multiplying from left/right often makes little sense. That can often better be calculated as directional derivatives.

Your code looks very messy. The first thing you should do is probably to write it down cleanly on paper and simplify the expressions before throwing them on CasADi.

Best regards,

Joel

[Adrian Bürger]

Hello Joel,

thank you for your hints!

 

You should try to use CasADi matrices instead of numpy matrices. numpy.zeros, numpy.eye give dense matrices. SX.eye, MX(3,4), etc. are sparse matrices. If you want to use "standard" datatype, you could use e.g. scipy.sparse.csc_matrix, but best is probably to work with CasADi matrices directly. Note that SX(2,3) is an all-zero sparse 2-by-3 matrix, SX.zeros(2,3) is a dense 2-by-3 matrix. This should make a huge difference.

When using the CasADi sparse matrices, the computations are much faster, this is great!

 

Other things: Calculating a Jacobian and multiplying from left/right often makes little sense. That can often better be calculated as directional derivatives.

I'm afraid I do not really understand what you mean, I'm sorry... can you show to me this with a little example?

 

Your code looks very messy. The first thing you should do is probably to write it down cleanly on paper and simplify the expressions before throwing them on CasADi.

This formula for computing the covariance matrix to achieve a statement of the goodness of the fit for estimated parameters can be found in several texts of H. G. Bock and related people and it's always the same formulation, I'm not sure if it can be simplified any further in, at least if I am able to do this. This is the complete expression:

Please see also the upcoming response to Greg's answer!

Best regards, Adrian

[Adrian Bürger]

Hello Greg,

thank you for the hint with the sparse solver! But as always, after one question got answered, two others arise:

1. I tried casadi.solve(A, b, "csparse") with A, b SX symbolics, and recieved the following error:

NotImplementedError: Wrong number or type of arguments for overloaded function 'solve'.

  Possible Python usages are:

    solve(array(double) ,array(double) ,str,Dictionary)

    solve(array(double) ,array(double) ,str)

    solve(array(int) ,array(int) )

    solve(array(double) ,array(double) )

    solve(SX,SX)

    solve(MX,MX,str,Dictionary)

    solve(MX,MX,str)

    solve(MX,MX)

You have: solve(SX, SX, str)

So for some reason, you can only supply an additional string when using doubles or MX symbolics, but not when using SX symbolics. Is this intended?

2. I then tried to use the LinearSolver-class as you suggested, and as far as I can tell it's working, but I was able to discover something else: when I set up the collocation method as depicted within the CasADi example pack in vdp_collocation2.py, I will initialize variables for the states that will not be used within the NLP, but handed to the solver. This happens because when initializing

V = struct_symMX([

      (

       entry("X",repeat=[nk+1,d+1],struct=states),

       entry("U",repeat=[nk],shape=nu)

      )

    ])

also for the last time point, the variables for the states will be repeated as for the other points in time, and therefor the last d variables will never be used. Within to the minimization, these variables will then just be set to zero, and not cause any problem here (I think). If I now (with an example model that contains 4 states) try to do as follows:

h3 = ca.vertcat((h1, h2))

lsol = ca.LinearSolver("csparse", h3.sparsity())

I receive an error message telling me

RuntimeError: The assertion "!sparsity.isSingular()" on line 48 of file "/home/casadibot/slaves/linuxbot/release-update-minimal/source/casadi/core/function/linear_solver_internal.cpp" failed. 

LinearSolverInternal::init: singularity - the matrix is structurally rank-deficient. sprank(J)=2566 (in stead of 2578)

which is a discrepancy of  12 lines, which matches "4 states repeated 3 times to often". For another example I tried that contains only two states, the discrepancy is 6, which fits to my considerations. So it seems that I now have to do one of the following things:

a) either I have to prevent the states from being repeated to often, so I would have do something like "append only one more state-vector at the end of the X-entry" or

b) I have to calculate the derivatives w. r. t. all but these variables by somehow "filtering out" the last 3 repetitions values from the struct.

For both things, I do not know how to achieve them :-) but I think that a) would be the much cleaner solution. Can this be achieved in some elegant way?

Thank you for your answer!

Best regards,

Adrian

[Joel Andersson]

Hi Adrian and sorry for slow answer!

Some remarks:

• I don't understand why you use "nlpIn/nlOut" to construct the expressions. From what I can tell, they won't be passed to an NLP solver.
  
• About the linear solve: You create an inverse, and then multiply the inverse with your right-hand-side, i.e. "mul(solve(A, eye), rhs)". It is more efficient and more stable to solve for the right-hand-side, i.e. "solve(A, rhs)"
  
• You use a matrix multiplication with "mul(horzcat(eye(n), zeros(n,m)), A)" to get the first n rows of a matrix. It should be more efficient to use slice: A[:n, :]

Have a look at https://gist.github.com/jaeandersson/eb0c5203806b25ed3ae9. Maybe that helps you further.

Best regards,

Joel

[Joel Andersson]

Hi!

Concerning this second message, solve(A, b, "csparse") is the correct syntax, but it only works for MX not for SX. SX doesn't support embedded function calls (which a call to a linear solver is). The "solve" algorithm used in SX is different, it uses a QR factorization with some reordering to avoid zeros along the diagonal, but without partial pivoting.

Best regards,

Joel

[Joris Gillis]

Hello Adrian,

For scalability of your method, it is better to focus on the MX-oriented solution.

It is not possibly to specify a special-case repeat as you describe.
One can however do as follows:

```

V = struct_symMX([

      (

       entry("X",repeat=[nk,d+1],struct=states),

                        entry("XF",struct=states),

       entry("U",repeat=[nk],shape=nu)

      )

    ])

```

Or alternatively:

```

V = struct_symMX([

      (

       entry("X",repeat=nk+1,struct=states),

                        entry("X_helper",repeat=[nk,d],struct=states),

       entry("U",repeat=[nk],shape=nu)

      )

    ])

```

Best regards,
  Joris Gillis

[Adrian Bürger]

Hello Joel, Joris and Greg,

sorry for the late reply, but with your combined help, I was now able to solve my tasks!

Thank you again very much for your help! I'll mark this one as solved.

Best regards,

Adrian