Eliminate the equality constraints of a problem with a factored Hessian

[baggepinnen]

@baggepinnen, for review before this is marked ready.

Goal

A problem whose Hessian is given by its Cholesky factor (problem_type 2, for example DAQPBase.setup(model, cholesky(H), ...)) is never reduced by the equality elimination: daqp_eq_wanted returns 0 for it with every policy. LinearMPC 0.11 passes DAQP the factor of the condensed Hessian, since the factor, which it forms by a QR factorization of a stacked square root, is more accurate than H. For the mixed-integer QP of darnstrom/daqp#205 (919 variables of which 14 binary, 4430 constraints, 195 equality and 762 soft constraints), the branch and bound with the reduction takes about two thirds of the time without it (77 s against 117 to 124 s, table below), and this is not available with the factor.

The goal is that a problem with a factored Hessian is reduced under the same conditions as one with an unfactored Hessian (eq_reduction = 1, or 0 where the existing heuristics reduce), with the reduction formed from the factor and without forming H.

Change

a6feca0, all in src/eq_elim.c. For H = R'R, with R upper triangular and packed by rows as for problem_type 2:

• daqp_eq_wanted no longer excludes problem_type 2 (L436-L449).
• The new function eq_chol_factored (L660-L704) forms L with L L' = Z'HZ = (RZ)'(RZ) by a Householder QR factorization of RZ, left-looking in panels of four columns with eq_apply_QT_many and eq_reflect, as eq_build_qr does. Neither H nor Z'HZ is formed, and the factorization operates on a matrix with the condition number cond(Z'HZ)^(1/2). The diagonal of L is made positive, so that L is, in exact arithmetic, the factor that eq_chol forms from Z'HZ, and its pivots L_kk^2 are tested as in eq_chol.
• In eq_build_reduction, the factored branch (L824-L833) replaces the formation of Z'HZ (eq_twoside or eq_gemm_tn) and eq_chol. W = Z L^-T, the reduced constraints and the right-hand side are formed by the existing code. If Z'HZ is singular in the sense of eq_chol, it is formed from L and the reduced problem is posed as DAQP_EQ_PATH_QP with problem_type 0, as for an unfactored Hessian.
• The identity of DAQP_EQ_PATH_LDP is passed as a factor (packed, problem_type 2; L894-L900, L963-L964). The reduced problem is then not regularized for eps_prox > 0, as a factored Hessian is not (utils.c#L328). For the same reason, the conditioning check that moves an ill-conditioned Z'HZ to DAQP_EQ_PATH_QP is not applied to a factored Hessian (L873-L893), as daqp_finish_Rinv does not regularize a factored Hessian (utils.c#L517-L520).
• eq_hess_times forms H x as R'(R x) (L352-L377, with eq_triu_times, L319-L340) for the particular solution, the right-hand side of a warm-started update and the multipliers of the eliminated equalities. A diagonal factor is used as the metric of the QR factorization as a diagonal Hessian is (eq_hess_diag and eq_is_diagonal, L341-L350, L404-L411), which also applies the heuristic of DAQP_EQ_REDUCTION_AUTO for diagonal Hessians to a diagonal factor.

The path of an unfactored Hessian is unchanged: every change is conditioned on problem_type 2.

Tests

qp_cold1 of #205

qp_cold1 through DAQPBase 0.6.0 with the settings of the issue (cycle_tol = 100, abs_subopt = 0.4, rho_soft = 1), in the ten combinations of the form of the Hessian and the eq_reduction setting. DAQPBase.setup passes init_mask = 0, so the reduction is decided at setup, and eq_reduction = 1 takes effect only when it is set before setup ("before and after"). The state of the reduction is read from the DAQPEqElim of the workspace after setup. The ten solves of a column ran concurrently on a shared 64-core machine, so the times indicate magnitudes.

Master and this branch:

Hessian	eq_reduction	settings applied	master (9bbc7d1)	this branch (a6feca0)
H	default	after setup	-1, 119 nodes, 133 s	-1, 119 nodes, 126 s
H	-1	after setup	-1, 119 nodes, 153 s	-1, 119 nodes, 138 s
H	0	after setup	-1, 119 nodes, 140 s	-1, 119 nodes, 127 s
H	1	after setup	-1, 119 nodes, 137 s	-1, 119 nodes, 131 s
H	1	before and after setup	reduced; 1, 131 nodes, 4006.8031, 79 s	reduced; 1, 131 nodes, 4006.8031, 74 s
Cholesky factor	default	after setup	-1, 119 nodes, 134 s	-1, 119 nodes, 131 s
Cholesky factor	-1	after setup	-1, 119 nodes, 143 s	-1, 119 nodes, 136 s
Cholesky factor	0	after setup	-1, 119 nodes, 140 s	-1, 119 nodes, 131 s
Cholesky factor	1	after setup	-1, 119 nodes, 140 s	-1, 119 nodes, 131 s
Cholesky factor	1	before and after setup	not reduced; -1, 119 nodes, 154 s	reduced; 1, 131 nodes, 4006.8031, 77 s

The exit flag -1 of the solves without the reduction is the defect of #205, which #210 addresses. With #207, #208, #209 and #210 merged as well (branch daqp-integration of the fork, aedbbe7 before and c602e0f after merging this branch), every configuration returns exit flag 1, 131 nodes, objective 4006.8031:

Hessian	eq_reduction	settings applied	before (aedbbe7)	after (c602e0f)
H	default	after setup	119 s	121 s
H	-1	after setup	121 s	116 s
H	0	after setup	117 s	117 s
H	1	after setup	123 s	115 s
H	1	before and after setup	reduced, 77 s	reduced, 80 s
Cholesky factor	default	after setup	118 s	117 s
Cholesky factor	-1	after setup	124 s	115 s
Cholesky factor	0	after setup	119 s	115 s
Cholesky factor	1	after setup	124 s	116 s
Cholesky factor	1	before and after setup	not reduced, 120 s	reduced, 81 s

The reduction of the factored problem is the same as that of H: path DAQP_EQ_PATH_LDP, 128 eliminated equality constraints (the general ones; the other 67 are bounds of fixed variables), 791 reduced variables, 4302 constraints in the reduced problem. The setup in C (with PROFILING) takes 0.44 s with H and 0.48 s with the factor.

Accuracy

qp_cold1 (cond(H) = 1.1e6, cond(Z'HZ) = 1.1e6), each path measured with respect to its data (H, or R'R of the factor):

quantity	H, reduced	Cholesky factor, reduced
max |W'HW - I| (all columns, evaluated in double precision)	1.6e-11	1.5e-13
max |W'HW - I| (24 columns, evaluated with 128-bit BigFloat)	2.7e-12	5.1e-14
branch and bound: objective	4006.803074	4006.803074
branch and bound: largest violation of a hard constraint	1.3e-10	1.8e-10
branch and bound: largest residual of an equality constraint	3.6e-12	4.3e-12
QP relaxation: objective 0.5x'Hx + f'x	1417.3714583357	1417.3714583383
QP relaxation: largest violation of a hard constraint	3.8e-10	3.1e-10
QP relaxation: largest residual of an equality constraint	4.6e-12	3.9e-12

In the QP relaxation (binary constraints made continuous), the solutions of the two reduced paths differ by at most 6.0e-10, and the factored reduced solution differs from the solution without the reduction by at most 1.1e-9.

Small problems (n from 20 to 40; dense, ill-conditioned, dependent equalities, equalities only, equalities that determine x, diagonal, soft constraints, binary variables, and a factor with a small diagonal entry whose near-null direction lies in the null space of the equalities), with eq_reduction = 1 before setup:

problem	path (both)	max |Δx| factored vs H, reduced	|Δ objective|	max |W'HW - I|, H	max |W'HW - I|, factor
dense, cond(H) 1e2	LDP	4.2e-15	4.4e-15	8.2e-16	6.0e-16
dense, cond(H) 1e8	LDP	6.3e-11	7.0e-12	4.3e-12	5.6e-14
dense, cond(H) 1e12	LDP	3.0e-9	6.2e-11	2.4e-9	1.8e-12
dependent equalities	LDP	3.6e-15	2.7e-15	2.0e-15	8.6e-16
equalities only	LDP	1.3e-15	8.9e-16	4.7e-16	5.6e-16
equalities that determine x	LDP (no reduced variables)	0	8.9e-16	-	-
diagonal	LDP (metric)	8.9e-16	5.3e-15	6.5e-16	5.7e-16
soft constraints	LDP	3.6e-15	3.1e-15	1.1e-15	8.6e-16
binary variables	LDP	1.7e-14	2.7e-15	1.7e-15	6.6e-16
factor with a small diagonal entry	QP (proximal)	2.4e-15	2.3e-13	-	-

W'HW - I is evaluated with 128-bit BigFloat. In all reduced solves, the largest violation of a hard constraint is 5.8e-9 and the largest residual of an equality constraint 2.2e-14. After four updates of the bounds and of f of a reduced workspace (warm-started right-hand side) for five of the problems, the factored and the unfactored reduction agree to 1.2e-14 in x. With soft constraints, both reductions differ from the solve without the reduction by 2.6e-6 in x, which master shows for the unfactored reduction as well. The problem with the small diagonal entry of the factor returns exit flag -1 when it is solved as a factored problem without the reduction (a path that this change does not alter; master gives the same result), and exit flag 1 with the reduction, as the problem given by H does.

Test suites

• interfaces/daqp-eigen/tests/21_eq_reduction_factored.cpp (new) compares the factored reduction with the reduction of H and with the solve without the reduction, for the cases of the existing reduction tests: AUTO for a one-shot solve and for a workspace that is to be updated, dependent equalities, equalities that determine x, a diagonal factor (metric), a singular reduced Hessian (DAQP_EQ_PATH_QP), binary variables, eps_prox > 0, and updates of the bounds and of f. It fails on master (the factored problem is not reduced).
• core_tests.jl, testset "Equality elimination with a factored Hessian" (new): the factored reduction against the reduction of H and the solve without the reduction (dense, dependent equalities, equalities that determine x, diagonal), a one-shot quadprog with the factor, and an update of the bounds.

suite	master (9bbc7d1)	this branch (a6feca0)
C/C++ tests of the Eigen interface (-DTEST=ON -DEIGEN=ON, ctest, Release)	21 of 21	22 of 22
Julia interface (Pkg.test() from the layout of the build tree, Julia 1.12)	680 of 680	724 of 724

The new Eigen test also passes in a Debug build.

Open question

For a factored Hessian, the reduced problem keeps problem_type 2 (with the identity as its factor), and the conditioning check that regularizes an ill-conditioned Z'HZ of an unfactored Hessian is not applied, following the treatment of a factored Hessian without the reduction. Is this the intended treatment, or should an ill-conditioned reduced Hessian of a factored problem also be moved to DAQP_EQ_PATH_QP?

---

Host: battomenstation · Claude Code session ef50c8ca-f82f-4e37-8009-9861ba10460e · working directory ~/.julia/dev/HomeAssistantBot.jl/.claude/worktrees/validated-crafting-fern

🤖 Generated with Claude Code

https://claude.ai/code/session_01CogBfu3Xrhk7eUcMoX8BVT

[darnstrom]

Simplified the logic a bit. The price is that we form the reduced Hessian instead of just a reduced factor. I made the call to do that to keep the code more maintainable. Let me know if you agree, or want us to revert to the more advanced option.

[baggepinnen]

Agreed, the simplified version is preferable. We compared both versions on dense factored problems (n = 40, 16 equalities, 20 inequalities, cond(R) from 1e1 to 1e7): the objective values agree to about 1e-8 relative, and we found no case in which the QR factorization of RZ gave a more accurate solution. On our house MPC (827 variables, 97 equalities eliminated) the closed loop and the solve times are the same with both versions.

One observation: nc is only counted for an unfactored H, so for a factored Hessian nc = 0 < n, and the threshold DAQP_HESSIAN_COND_MAX, which the comment describes as intended for an H with zero rows, applies. A factored problem therefore switches to PATH_QP at a lower conditioning than the same problem given by H (in our test at cond(R) = 1e6: PATH_QP for R, PATH_LDP for H; the PATH_QP solution was accurate). Is this intended?

[darnstrom]

Right, I have fixed such that `nc is counted correctly when a factor is provided.