forum.alglib.net

Currently, I try to solve the following linear system with the QP-BLEIC solver

AtA * x = AtB, subject to:
0<= x_i <= 1.0 and
sum(x_i) = 1.0;

In order to prove my implementation, I have a test case, that creates 100s of equations systems with AtA of size 2x2 to 4x4. In the final case AtA will never be bigger than 8x8. My initialization values are close to the final solution. To validate the result I compute x = AtA^(-1) AtB with Eigen. However, the results returned by Alglib are far from the initial and unconstrained solution. Following the examples provided with Alglib I set:

all n entries in bndl to 0.0
all n entries in bndu to 1.0
all n entries in scaling to 1.0
all n+1 entries in linear constraint (c) to 1.0
type of linear constraint is ct ="[0]"
with n being the size of vector x

I have tried default termination values as well as low error bounds (1e-10) and high maximum iterations (1000). In most cases, the rep variable return: inner/outeriterations=1 or 2; terminationtype=4 (success). Have I missed something?

Below one example:

Matrix AtA:
0.422638 0.430536 0.441594
0.430536 0.470939 0.437771
0.441594 0.437771 0.545798

Vector AtB:
0.460166 0.485326 0.499367

Starting Values:
0.782251 0.213387 0.00436184

AtA*x_(start):
0.424406 0.439189 0.441233


Solution AlgLib:
1 0 -6.27744e+66

Solution Eigen (unconstrained):
0.126134 0.627357 0.30969

Hi!

1. Are you sure that you set ALL the necessary constraints and parameters (including problem size) correctly?

The trailing ...E66 in "1 0 -6.27744e+66" is very suspicious - almost all ALGLIB solvers (including QP-BLEIC) respect box constraints. So returning value which is so much out of bounds may indicate that you (say) accidentally created 2-dimensional problem instead of 3-dimensional one.

2. I just tried your data with QP solver, and it returned an answer with all constraints being satisfied.

However, I have to point that your problem is formulated as "solve linear system subject to constraints". But QP problem formulation is: "minimize x'Ax+b'x subject to constraints". Simply feeding your A and b to QP solver will return some answer, but it will be meaningless answer - a correct solution to completely different problem.

In particular, ALGLIB returned [1, 0, 0], which satisfies all constraints and seem to minimize x'Ax+b'x subject to constraints - but it is far from what you want it to do: to solve linear equations system.

So, you should reformulate your problem in QP-like way if you want to use QP solver.

Hello Sergey,

indeed the error is on my side. I expected that the size of the solver is derived from the input vectors/matrices and forgot to adjust the size in minqpcreate(size, state). This solves the issue of getting values above the box constraints.

Thank you also for pointing me that my problem formulation is incorrect! I will have a more detailed look at QP optimization. If you have a good recommendation for people who haven't a higher degree in math, feel free to share. :)

Hi!

1. You can treat this task as minimization of 0.5*|AtA*x-AtB|^2, which is equal to 0.5*(AtA*x-AtB)^T * (AtA*x-AtB) = 0.5*x'(AtA'*AtA)x - (AtB*AtA)x + const. In this form solution of QP problem will minimize error in the right part subject to constraints.

2. From the other side, this "AtA / AtB" notation suggests that your linear equations system can be result of normal equations approach to some least squares problem A*x~B. If it is the case, then it means that you "squared" the least squares system in order to solve it, but now you have to square it one more time. As result, condition number can become prohibitively high.

The better approach in such cases is to minimize directly 0.5*|Ax-B|^2 = 0.5*x'(A'A)x - (B'A)x + const.

Hmm... it looks exactly as what you tried to do!!! :) Well, this idea may work - but only provided that your AtA/AtB are in fact A'*A and A'*B for some A and B. It does not have to work for "general" AtA and AtB.