Hoist loop-invariant scaling reciprocal out of the LDP inner loops

[cl445]

Motivation

In the QP-to-LDP transform, daqp_update_M and daqp_update_v divide by work->scaling[...] inside the inner loop, once per element. The divisor is the same across that loop; it depends only on the column, not the inner index. Floating-point division is several times slower than multiplication on most targets, and the gap is widest on the in-order ARM cores used in Raspberry Pis (Cortex-A53/A72), where these dense O(n^2) loops dominate the transform.

Change

Compute the reciprocal once per column and multiply:

for(; j<n; ++j){// Take into account scaling in Rinv
    c_float inv_scaling = 1/work->scaling[n-j-1];
    for(i=0;i<j;++i)
        work->M[disp2-i] += (work->Rinv[--disp]*inv_scaling)*A[disp2-j];
    work->M[disp2-j]=(work->Rinv[--disp]*inv_scaling)*A[disp2-j];
}

and the analogous loop in daqp_update_v. This removes the O(n^2) divisions from the scaling branch and keeps the same multiply order. These branches run on warm updates that reuse the H factorization (the MPC-style path), not on a fresh solve.

Testing

Correctness: full Python test suite passes (21 tests), including test_model_update_bounds, test_model_update_cost_vector, and the warm-start tests that drive these branches.

Speed: standalone C microbenchmark calling daqp_update_M + daqp_update_v in a loop (no Python or Cython overhead), baseline vs this branch, identical compiler flags, median of 3 runs on Apple Silicon (arm64):

n	baseline	this branch	faster
40	12.8 us/iter	11.9 us/iter	~7.5%
80	101.6 us/iter	90.3 us/iter	~11%
160	~889 us/iter	691 us/iter	~19%

The relative gain grows with n: the removed divisions are O(n^2) while the fixed per-call cost is lower order. This is an out-of-order core where FP division is comparatively cheap, so an in-order Cortex-A (Raspberry Pi 3 / Zero 2, where FP64 division is not pipelined) should gain more. The benchmark measures the transform in isolation, not an end-to-end solve.

Downstream CI

PolyDAQP.jl fails on the x86 CI runner, while DAQP's own C, Python, Julia, and MATLAB suites pass, as do the other three downstream packages (ASCertain, LinearMPC, ParametricDAQP). The failing assertion is in the elimination test, which builds a random 6D polytope, projects out three variables, and checks the bounding box matches:

all(abs.(ub[1:3] - ubn) .< δ)   with δ = 1e-8

I reproduced this locally on arm64, running the elimination test over 2000 random seeds with both the unmodified master build and this branch:

• master and this branch give bit-identical bounding boxes on all 2000 seeds (max |Δ| = 0). On this platform the rewrite changes nothing downstream, because the perturbation is below what tips an active-set decision.
• The test still fails on 8 of the 2000 seeds (0.4%) on unmodified master, with the bounding box off by up to 2.5e-4, far above δ = 1e-8.

So the failure is a pre-existing over-approximation in PolyDAQP.eliminate, not a regression from this change: on borderline random polytopes the projection comes out too large by up to ~2.5e-4, which the test's δ = 1e-8 does not tolerate. It is not a solver-tolerance knob either (tightening primal_tol in the redundancy step does not reduce it and makes some instances worse), so it looks like the iterated Fourier-Chernikov step itself. The CI runner is x86, where the transform rounds differently than arm64, so the one seeded CI instance happens to tip there; the same class of instance already fails on master for about 0.4% of seeds.

I filed the analysis and a reproduction as darnstrom/PolyDAQP.jl#3. This PR does not affect that behavior (bit-identical downstream on arm64), so the red check should clear once the issue is addressed, by improving eliminate or relaxing the test's δ.

[darnstrom]

Thank you! This optimization make sense.

With this type of optimization, I think that one can break out yet another factor from the inner loop. So instead of inv_scaling one might use col_scaling as this:

for(; j<n; ++j){// Take into account scaling in Rinv
    c_float col_scaling = A[disp2-j]/work->scaling[n-j-1];
    for(i=0;i<j;++i)
        work->M[disp2-i] += work->Rinv[--disp]*col_scaling;
    work->M[disp2-j]= work->Rinv[--disp]*col_scaling;
}

Do you think this make sense @cl445? I added it in 5556799

[cl445]

Yes, that makes sense to me. It’s strictly better than what I had.
It’s actually closer to the original numerically.