IterativeRefinementSolver performance

[apete]

@Programmer-Magnus

In the IterativeRefinementSolver the primal and dual scaling factors are always treated as independant and different from each other.

Solving quadratic programs to high precision using scaled
iterative refinement:

Note that all calculations except the expensive solves of the QPs are done in rational
precision. Algorithm 1 uses only one scaling factor Δk for primal and dual infeasibility
to avoid the scaling of the quadratic term of the objective. Keeping this matrix and
the constraint matrix A fixed gives QP solvers the possibility to reuse the internal
factorization of the basis system between refinements, as the transformation does not
change the basis. Hence one can perform hotstarts with the underlying QP solver
which is crucial for the practical performance of the algorithm. This comes at the cost
of only scaling either primal or dual infeasibilities as required, especially if they differ a lot, possibly slowing convergence.

I think treating them as two different or one common scaling factor should be configurable, and one common should be the default.

Provided you fully exploit the fact that they are common this should be a huge performance gain.

[Programmer-Magnus]

If this will make calls, after the first, to the ConvexSolver faster, it should be tested. I just need to read the article again. I guess that the change is very small.

I have another experimental idea for this routine. I believe the Quadruple precision type might be replaced by new matrix-vector methods. I hope it will give the same precision by using Kahn summation (or version of it). That might lead to much less memory being allocated. This will require more coding at a later time.
The rational is that Quadruple precision is only used to compute the residual at higher precision. If that can be achieved in another way, It can be replaced. It probably even work to do this in normal double precision, but with less improvement. That should also be tested.
https://en.wikipedia.org/wiki/Kahan_summation_algorithm

[apete]

The key performance gain, I assume, will come from not creating a complete new sub-problem at each master iteration. Doesn't matter if you use Quadruple or something else.

Having one common scaling factor is required to avoid recalculating the Q-inverse. (Deriving that scaling factor; my first attempt would be to only look at the primal infeasibilities. Or why not just scale it with something like 1_000 every time.)

The IterativeRefinementSolver should use one instance of ConstrainedSolver (some subclass of that) that it manipulates. The main work here is building what's missing to be able to manipulate and restart the ConstrainedSolver . All we need to do is replace the c, be and bi matrices and deactivate inequality constraints that are now infeasible (something like that) and then re-solve from there.

[apete]

I sketched out an implementation in a branch https://github.com/optimatika/ojAlgo/tree/refinement

Just got curious to try this and didn't want to mess with your code.

The tests are passing, but I have no idea about performance. Perhaps you can try this on some of your problems. (Haven't really done any performance optimisations yet.)

Found and fixed a few surrounding issues, like loss of precision when extracting the solutions, that should be moved to the main branch regardless of the rest.

It is now possible to solve problems to higher than double-precision. That wasn't possible before, because there was always some intermediate calculation done in double-precision.

[Programmer-Magnus]

Great. I will have a look at the refinement branch.

[Programmer-Magnus]

Looks more optimized than present code since checking for zeroes.

1. For me testQP1 works better when replacing,

        double errorI = iterBI.aggregateAll(Aggregator.LARGEST).doubleValue();
   

with
double errorI = iterBI.negate().aggregateAll(Aggregator.MAXIMUM).doubleValue();

It computes the max violation in inequalities. If iterBI is positive its fine, but it should not be negative.

Random comments:
2. I would find it useful to find/save an error norm somewhere to not have to compute it from the result afterwards.
I can accept large errors in the solution if they are rare. Today I check the error when ojAlgo say it failed, and if error is small enough I still accept it. But if possible I want a more precise result.

3. The iterations with noTries in end of IterativeRefinementSolver tries to take a big step in the scaling, going back if solution fails. If a big step mostly succeed this saves iterations.

IterativeRefinementSolver uses relative errors to try to handle large values in Q and C.

[apete]

Yes, the inequalities error should be one-sided. (fixed that)

Would it be enough if you could trust the state to indicate the quality of the solution? With larger errors set the state to APPROXIMATE.

I got the impression IterativeRefinementSolver2 only needed a few iterations to complete, but I'm sure the strategy can be improved. Just don't want to complicate too much before it's actually needed.

[apete]

@Programmer-Magnus With your recent changes to IterativeRefinementSolver and IterativeRefinementSolverDouble or my attempt at IterativeRefinementSolver2, are there any performance gains? I mean actually faster solve times on your real life problems.

[apete]

@Programmer-Magnus Curious about this.

[Programmer-Magnus]

I just wanted to say that I am working a bit on this now. But it is not as easy as it sounds. The first tiny benchmark showed no difference in speed.