Sparsity masks for sparity exploting solvers

Feature Requests
1

Hi everyone,

I would like to discuss a feature that could unlock the full potential of the sparse QP solvers in acados, OSQP and Clarabel, in particular when using them inside SQP.

Current situation

OSQP and Clarabel work on sparse matrices in CSC format and exploit sparsity in all their linear algebra. When building these matrices, acados already exploits the stage structure of the OCP, asonly the stage blocks are stored, nothing between them. Within each block, however, all entries are written, including the zeros. In typical MPC problems many of these entries are structurally zero. The cost matrices are often diagonal, and dynamics and constraint matrices frequently contain many zeros. The solvers then factorize and multiply with these zeros as if they were nonzeros.

Proposed feature

acados could allow defining sparsity patterns for the matrices A, B, Q, S, R, C, D (dense by default). In contrast to the one-sided constraint masks, which can change between solves, the patterns would be fixed at solver creation (e.g. as part of the dims), since they determine the memory and the symbolic factorization of the sparse solvers. The interface code of the sparse solvers would then build the CSC matrices only from the entries in the pattern.

  • In the QP interface, the user would set the patterns together with the dimensions.
  • In the NLP interface, the patterns could come from the sparsity of the CasADi functions.
  • The modules between NLP and QP solver (integrator, Hessian approximation, regularization, condensing) would carry the patterns through at solver creation. Some steps are less straightforward (e.g. fill-in in A/B when integrating continuous dynamics, exact Hessians through the integrator, eigenvalue-based regularization like MIRROR/PROJECT), but in the worst case a module can always report a dense block.

Initial results

As we are currently comparing QP solvers, our student @jmiddelberg prototyped OSQP and Clarabel interfaces that detect the pattern from the data . On our humanoid walking MPC ( a QP with N = 9, 13 states, 12 inputs, 32 general constraints per stage), it reduces the median solve time by a factor of about 3:
The data-driven detection works, but each rebuild costs time and the warm start. A pattern defined at creation avoids that.

Questions

  • What do you think of such a feature?
  • Where would such a pattern live at best - in the dims?
  • Would you accept a stepwise contribution, starting with the QP interface and the OSQP/Clarabel interfaces, and extending to the NLP interface later?

Curious to hear your opinions!

Best,
Franek

2

Hi Franek,

I think this makes sense, if we you want to use these QP solvers efficiently!
Of course it makes sense to start with the QP solvers.

I think the pattern could live in the `in` structure of the module. This way, we are more flexible in creating the sparsity pattern later as the memory does only depend on the dims. The goal is to improve runtime and the potential memory savings are not so important. When using the QP solver within an NLP solver, we anyway have a dense copy of the OCP-structured problem.

Best,
Jonathan

3

Hi Jonathan,

okay that makes sense to me. The ocp_qp_in struct is currently part of HPIPM. Does it make sense to add the sparsity pattern to into HPIPM?

Best,
Franek

4

Hi Franek,

ocp_qp_in is defined in acados but only as an alias to the HPIPM struct:

typedef struct d_ocp_qp ocp_qp_in;

I am not sure yet where the sparsity structure should live.

If we only want to change the ocp_qp_* definition in acados, we could do something like

typedef struct
{
    struct d_ocp_qp_dim hpipm_dims;
    int with_sparsity;
} ocp_qp_dims;

typedef struct
{
    struct d_ocp_qp hpipm_ocp;
    struct ocp_qp_sparsity *sparsity;
} ocp_qp_in;

However, this will require a lot of changes in acados, which could be avoided when adding the sparsity information directly in HPIPM.

@giaf should comment if he would “allow” adding the sparsity information in HPIPM QP data.

Best,

Jonathan

5

Hi all,

in general I would expect sparse QP solvers to not be competitive with structure exploiting ones for most problems. My feeling is that well optimized dense linear algebra would be maybe an order of magnitude faster than sparse one: so I would expect to see sparse QP solvers to perform better only on really sparse problems.

On top of that, sparsity is often lost through processing. For example, even if the ODE is very sparse, the sensitivities returned by the integration methods tend to be much less sparse. Another example is partial condensing, where the block condensing would further destroy structure (so you should take care of properly processing or just invalidating the sparsity structure with performing partial condensing within acados). And in problems with long horizon and little controls, partial condensing can help structured solvers such as HPIPM a lot.

That said, there can still be cases where sparse QP solvers can be the best option, if the sparsity structure within the data matrices is exploited. In HPIPM there is only a field used to flag for a diagonal Hessian, and I don’t remember if this is processed systematically through e.g. (partial) condensing, I need to double check that.

As of now I would not add explicitly sparsity information to HPIPM, because this would not be used in HPIPM itself, and it would not be practical to maintain a structure that is only used elsewhere.

But I can add a void pointer to the dim and qp structures, such that you can use it however you like to pass stuff around together with the QP structures, without having to change acados internals.

In acados we have already been doing that for quite some time: HPIPM solution structs have a member void *misc that is unused in HPIPM, and it is used in acados to pass around the QP info structure. acados takes exclusive care of this structure, comprising its memory allocation.

Would this make the job for you?

6

Hi all,

from my perspective an additional void pointer in qp_in_* would work, but @FreyJo you are the expert.
Regarding partial condensing I would imagine running the x0-elemination and partial condensing on the sparsity pattern (0/1 matrices), would obtain the sparisty pattern of the condensed problem.

For context: our aim is exactly to compare sparsity- against structure-exploiting QP solvers, especially also under consideration of condensing. Initial results on bipedal walking QPs have shown that OSQP’s solve time drops by a factor of 2.6 when the zeros inside the stage blocks are left out.
I am happy to share some numbers once I have them.

Best,
Franek

7

Hi all,

I agree using void *misc in qp_in_* and dims_* would be a good solution.

In principle, the sparsity pattern propagation should work. One pitfall could be that qp_solver_cond_ric_alg requires Hessian blocks Q with respect to x to be strictly positive definite, which a matrix full of ones does not satisfy.

Best,

Jonathan

8

Hi all,

The sparsity pattern out of the (partial) condensing algorithm routine is the same regardless of what condensing algorithm is used. So you can just use a condensing algorithm allowing for indefinite hessian.

That said, I would expect that in general the (partially) condensed Hessian and dynamics quickly get nearly full, except for one stage at the beginning of each block.

The (partially) condensed constraint matrix has some structure, being block triangular (or better, at the moment it consists of two block triangular matrices glued together, but it can be made into a single block triangular matrix by changing the order of general constraints in condensing in HPIPM, I should do that at some point). So here there are structural nonzeros to be exploited, but within each structural block sparsity would probably get lost quickly too.

9

Thanks @giaf for adding the void pointers to the HPIPM structs! As soon as I find some time I will try to implement the sparsity pattern handling into acados.