Issues when `tril` is chosen for a new ldlsolver

[mipals]

Hi Paul,

I am trying to add PanuaPardiso as a possible solver for the KKT system. Preferably I would like to use tril as in Julia (PanuaPardiso wants the upper part in a CSR format, so if I get tril in CSC I can simply swap the row/cols in order to get triu in CSR format). However, I've encountered two issues so far:

1. Adding the following to a new PanuaPardisoDirectLDLSolver result in an initial KKT system for which the diagonal is not completely filled (for a triu format the diagonal is filled - as is required by PanuaPardiso).

fn required_matrix_shape() -> MatrixTriangle {
      MatrixTriangle::Tril
  }

2. Later the code also panics (unsurprisingly) when it hits the following line https://github.com/oxfordcontrol/Clarabel.rs/blob/d1ca25da23bf3a925bcde3315a2f4d10361c02de/src/algebra/csc/core.rs#L302

I've tried to look into the code, but seems a little complicated, so my hope is that you might know an easy fix.

Best, Mikkel

[goulart-paul]

It sounds like a bug and not something you did. This would be unsurprising because I don’t think we have any solvers in the Rust implementation that want inputs int lower triangular form, so probably not stress tested. It works in Julia though, so probably something is just a bit out of sync.

I will take a look, but can’t promise to do so until next week. Are you doing this entirely within rust or trying to do so through the C++ interface? Probably it is a rust issue, but would be helpful to know as a starting point.

P

On 13 Jun 2024, at 16:42, Mikkel Paltorp @.***> wrote:

Hi Paul,

I am trying to add PanuaPardiso as a possible solver for the KKT system. Preferably I would like to use tril as in Julia (PanuaPardiso wants the upper part in a CSR format, so if I get tril in CSC I can simply swap the row/cols in order to get triu in CSR format). However, I've encountered two issues so far:

1. Adding the following to a new PanuaPardisoDirectLDLSolver result in an initial KKT system for which the diagonal is not completely filled (for a triu format the diagonal is filled - as is required by PanuaPardiso).

fn required_matrix_shape() -> MatrixTriangle { MatrixTriangle::Tril }

1. Later the code also panics (unsurprisingly) when it hits the following line https://github.com/oxfordcontrol/Clarabel.rs/blob/d1ca25da23bf3a925bcde3315a2f4d10361c02de/src/algebra/csc/core.rs#L302

I've tried to look into the code, but seems a little complicated, so my hope is that you might know an easy fix.

Best, Mikkel

— Reply to this email directly, view it on GitHubhttps://github.com/oxfordcontrol/Clarabel.rs/issues/117, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ABUFVKKLJMICB3QP6AMU6TTZHG4XNAVCNFSM6AAAAABJIVP74WVHI2DSMVQWIX3LMV43ASLTON2WKOZSGM2TCNBXGQ4TGNA. You are receiving this because you are subscribed to this thread.Message ID: @.***>

[mipals]

Hi Paul,

Im doing it entirely in rust. I think I've gotten it to work after investigating it a little further.

The assembly is actually working correctly. The problem ended up being that count_diagonal_entries() is not working for lower-triangular matrices. Since I used this function to check if the diagonal of the KKT system was there i reported the diagonal as missing some entires (which was actually there). The same bug seem to be present in the Julia implementation.

The second issue was fixed by compiling using --release as debug_assert! only applies the assertion when compiling in debug mode. However, I am not sure why it is there to begin with.

Cheers, Mikkel

[mipals]

At some point I would like to solve QPs with a direct solver that takes the tril rather than triu. I believe that the error in count_diagonal_entries() will become an issue in this case. It is an easy fix (counting the first rather than last element in a column in this case). I can make a PR here and for the Julia version if want.

Another thing I suspect can become an issue is that P should be given in the same MatrixTriangle as supported by the direct solver. For an unexperienced user that changes the underlying direct solver without the changing the form of P will run into issues.

[goulart-paul]

I am not clear about what the error condition is here. As I understand things, the solver is expecting to always get the input matrix P in CSC upper triangular form, regardless of the internal formatting of the KKT matrix specified for whichever direct KKT solver is used.

The particular choice of KKT matrix shape (i.e. triu or tril) is determined by the shape argument here, but shape there refers to the shape of the matrix being assembled, not the shape of the P matrix it is being assembled from.

The internal function count_diagonal_entries is only called from one place (i.e. during KKT assembly here) and only operates on P. I don't see what problem generalizing that function to allow for lower triangular inputs would be solving, because lower triangular inputs should never be passed to it regardless of solver type.

Note that KKT assembly code has some internal unit tests already that exercise this process, for both triu and tril KKT matrix formats (see here). If there is some bug in assembly of the KKT matrix to tril form, shouldn't one of those tests fail?

Possibly related: the solver is assuming that any inputs you pass to it are in CSC format with no repeats in and with entries appearing in ascending row order within each column. If that is not true that the assembly describing above could get mangled. I added some utilities for checking / fixing inputs that aren't in this standard form in #140, which should go into the next minor release. Maybe the issue is related to this instead?

P

PS. Apologies for only just around to this issue now.

(related PR is #143)

[mipals]

You are correct. In short the function currently does not work for tril matrices, but it does not matter since P is assumed to be given as triu (I think I got confused here and thought it had to match the underlying KKT solver). However, the reason I actually went back to this was that I anticipated that this would become an issue in the GPU version which currently assumes that a full matrix is given. Having looked at the GPU code again I can see that the function has already been modified to search the full column rather just the first element.

[goulart-paul]

OK. I made it a bit more general anyway, mostly so that the assumptions in the code would be clearer. Implemented in #145. I think #143 the way I have done it there is a bit faster since it just switches on a triu/tril marker, and therefore excludes (for now at least) the possibility of passing a full matrix. That's easily extended to a third (full) case when the time comes.

[mipals]

Looks good. I agree that it makes it more clear and will avoid confusion in the future.

[goulart-paul]

Fix included in v0.10.0.