albertosantini / quadprog

Module for solving quadratic programming problems with constraints

Hi there 👋

She sat behind the black pieces and said carefully in Russian, "Would you like to play chess?" 

-- Walter Tevis, The Queen's Gambit, 1983

quadprog's Issues

quadprog.js solveQP issue

Hello, Alberto Santini!

I’ve tried to use (and thank you for) quadprog.js from your repository and found that there are significant differences in your solution and original R solution (fortran-based library quadprod).
You have an example on that page that shows the input and output of solveQP. I’ve tried this example with the following data and got strange result that differed from R solution.

var Dmat = [], dvec = [], Amat = [], bvec = [], res;

Dmat[1] = [];
Dmat[2] = [];
Dmat[3] = [];
Dmat[4] = [];
Dmat[5] = [];
Dmat[6] = [];
Dmat[7] = [];
Dmat[1][1] = 0.00045232933192250496;
Dmat[2][1] = 0.0005097526372655511;
Dmat[3][1] = 0.0005725870424150262;
Dmat[4][1] = 0.0005048974812299229;
Dmat[5][1] = -1.116682787785219e-05;
Dmat[6][1] = 0.00019563544837137415;
Dmat[7][1] = 0.00033093688577306117;

Dmat[1][2] = 0.0005097526372655511;
Dmat[2][2] = 0.0006951277835926323;
Dmat[3][2] = 0.0006417805981316768;
Dmat[4][2] = 0.0006696702746960787;
Dmat[5][2] = -1.9542315203644345e-06;
Dmat[6][2] = 0.00024214025160680232;
Dmat[7][2] = 0.00038568575529115387;

Dmat[1][3] = 0.0005725870424150262;
Dmat[2][3] = 0.0006417805981316768;
Dmat[3][3] = 0.0008404721337760725;
Dmat[4][3] = 0.0006537024901320917;
Dmat[5][3] = -1.237174668388066e-05;
Dmat[6][3] = 0.00025408683849470514;
Dmat[7][3] = 0.00047162617359422656;

Dmat[1][4] = 0.0005048974812299229;
Dmat[2][4] = 0.0006696702746960787;
Dmat[3][4] = 0.0006537024901320917;
Dmat[4][4] = 0.0008531442886388518;
Dmat[5][4] = 1.6526537671997837e-05;
Dmat[6][4] = 0.0002641850198282845;
Dmat[7][4] = 0.0003581148091890727;

Dmat[1][5] = -1.116682787785219e-05;
Dmat[2][5] = -1.9542315203644345e-06;
Dmat[3][5] = -1.237174668388066e-05;
Dmat[4][5] = 1.6526537671997837e-05;
Dmat[5][5] = 6.024397455732469e-05;
Dmat[6][5] = 2.8301634815247453e-05;
Dmat[7][5] = -2.639682028670009e-05;

Dmat[1][6] = 0.00019563544837137415;
Dmat[2][6] = 0.00024214025160680232;
Dmat[3][6] = 0.00025408683849470514;
Dmat[4][6] = 0.0002641850198282845;
Dmat[5][6] = 2.8301634815247453e-05;
Dmat[6][6] = 0.00014176037608338397;
Dmat[7][6] = 0.00014458666076982795;

Dmat[1][7] = 0.00033093688577306117;
Dmat[2][7] = 0.00038568575529115387;
Dmat[3][7] = 0.00047162617359422656;
Dmat[4][7] = 0.0003581148091890727;
Dmat[5][7] = -2.639682028670009e-05;
Dmat[6][7] = 0.00014458666076982795;
Dmat[7][7] = 0.0022785517723467757;

Amat[1] = [];
Amat[2] = [];
Amat[3] = [];
Amat[4] = [];
Amat[5] = [];
Amat[6] = [];
Amat[7] = [];

Amat[1][1] = 1;
Amat[2][1] = 1;
Amat[3][1] = 1;
Amat[4][1] = 1;
Amat[5][1] = 1;
Amat[6][1] = 1;
Amat[7][1] = 1;

Amat[1][2] = 1;
Amat[2][2] = 0;
Amat[3][2] = 0;
Amat[4][2] = 0;
Amat[5][2] = 0;
Amat[6][2] = 0;
Amat[7][2] = 0;

Amat[1][3] = 0;
Amat[2][3] = 1;
Amat[3][3] = 0;
Amat[4][3] = 0;
Amat[5][3] = 0;
Amat[6][3] = 0;
Amat[7][3] = 0;

Amat[1][4] = 0;
Amat[2][4] = 0;
Amat[3][4] = 1;
Amat[4][4] = 0;
Amat[5][4] = 0;
Amat[6][4] = 0;
Amat[7][4] = 0;

Amat[1][5] = 0;
Amat[2][5] = 0;
Amat[3][5] = 0;
Amat[4][5] = 1;
Amat[5][5] = 0;
Amat[6][5] = 0;
Amat[7][5] = 0;

Amat[1][6] = 0;
Amat[2][6] = 0;
Amat[3][6] = 0;
Amat[4][6] = 0;
Amat[5][6] = 1;
Amat[6][6] = 0;
Amat[7][6] = 0;

Amat[1][7] = 0;
Amat[2][7] = 0;
Amat[3][7] = 0;
Amat[4][7] = 0;
Amat[5][7] = 0;
Amat[6][7] = 1;
Amat[7][7] = 0;

Amat[1][8] = 0;
Amat[2][8] = 0;
Amat[3][8] = 0;
Amat[4][8] = 0;
Amat[5][8] = 0;
Amat[6][8] = 0;
Amat[7][8] = 1;

Amat[1][9] = -1;
Amat[2][9] = 0;
Amat[3][9] = 0;
Amat[4][9] = 0;
Amat[5][9] = 0;
Amat[6][9] = 0;
Amat[7][9] = 0;

Amat[1][10] = 0;
Amat[2][10] = -1;
Amat[3][10] = 0;
Amat[4][10] = 0;
Amat[5][10] = 0;
Amat[6][10] = 0;
Amat[7][10] = 0;

Amat[1][11] = 0;
Amat[2][11] = 0;
Amat[3][11] = -1;
Amat[4][11] = 0;
Amat[5][11] = 0;
Amat[6][11] = 0;
Amat[7][11] = 0;

Amat[1][12] = 0;
Amat[2][12] = 0;
Amat[3][12] = 0;
Amat[4][12] = -1;
Amat[5][12] = 0;
Amat[6][12] = 0;
Amat[7][12] = 0;

Amat[1][13] = 0;
Amat[2][13] = 0;
Amat[3][13] = 0;
Amat[4][13] = 0;
Amat[5][13] = -1;
Amat[6][13] = 0;
Amat[7][13] = 0;

Amat[1][14] = 0;
Amat[2][14] = 0;
Amat[3][14] = 0;
Amat[4][14] = 0;
Amat[5][14] = 0;
Amat[6][14] = -1;
Amat[7][14] = 0;

Amat[1][15] = 0;
Amat[2][15] = 0;
Amat[3][15] = 0;
Amat[4][15] = 0;
Amat[5][15] = 0;
Amat[6][15] = 0;
Amat[7][15] = -1;

bvec[1] = 1;
bvec[2] = 0;
bvec[3] = 0;
bvec[4] = 0;
bvec[5] = 0;
bvec[6] = 0;
bvec[7] = 0;
bvec[8] = 0;
bvec[9] = -0.5;
bvec[10] = -0.5;
bvec[11] = -0.5;
bvec[12] = -0.5;
bvec[13] = -0.5;
bvec[14] = -0.5;
bvec[15] = -0.5;

var increment_by = 0.001, iteration = 21,
    step = increment_by * iteration;
dvec[1] = 0.003196462492322401 * step;
dvec[2] = 0.002010657315490846 * step;
dvec[3] = 0.0033973917499075636 * step;
dvec[4] = 0.001382156372148105 * step;
dvec[5] = 0.0013320376697428297 * step;
dvec[6] = 0.0019232309058722838 * step;
dvec[7] = 0.0021818304169943905 * step;


res = solveQP(Dmat, dvec, Amat, bvec);
//res.solution => [0.4772818853979638,0,0,0,0.5,0,0.022718114602036023]
//original R solution => [0.09037104540326746,7.304277169042432e-17,-5.3326210492157636e-17,-1.9732858652101513e-18,0.5,0.4005283006854242,0.0091006539113083]

R solution looks more smooth and reasonable. I am confused with the results and have a question:

• Do you know (or can you explain) why the solutions made by the algorithm in quadprog.js differ so badly from R solution? (At least the difference is in order - tenths vs hundredth)

I have no idea why is it happening and it seems like a bug.

P.S. In your another repository you use quadprog.js to optimize portfolios. Why do you use both R script and nativeGetOptimialPortfolio functions? Do you feel that native solution doesn’t work always properly?

Number.EPSILON support

Node.js currently doesn't support Number.EPSILON. Mozilla reports that it isn't available in several browsers and its value is about 2.2204460492503130808472633361816E-16.
I ran the previous version that included the vsmall calculation and found its value to be 1.4272476927059598e-15. Since javascript is numerically stable across platforms it is likely that this value will work in multiple environments. I obtained the same value using Safari browser on JSFiddle. I tested Number.EPSILON and 1.43e-15 as initial values to the vsmall calculation with surprising results. Assuming the vsmall calculation is related to the calculations used in the algorithm, then Number.EPSILON may be too small. When I used the Number.EPSILON constant in my test problems I obtained the same results.

Can you change line 99 from:
vsmall = Number.EPSILON;
to:
vsmall = Number.EPSILON || 2.2204460492503130808472633361816E-16;

Since vsmall is only used to test values near zero, it seems like Number.EPSILON is usable. In reviewing the code, I noticed that one comparison uses '<=' while the others use '<'.

Solver founds no solution though it should

Hello,
I'm using your QP. I have a set of solvable inputs. With the solver I get around 20% failure when output.Solution contains NaN. How can I tune (some threshold) the accuracy of the solver, so that I can get solution in 100% of solvable inputs?

My inputs are in format of
minimize 1/2 * x^T D x + d^T x
where A x >= b

I'm calling:
var output = solveQP(Dmat, dvec, Amat, bvec, meq, factorized);

One example of solvable input that returns NaN:
JSON.stringify(Dmat)
"[[],[0,2,0,0,0,0,0],[0,0,2,0,0,0,0],[0,0,0,2,0,0,0],[0,0,0,0,2,0,0],[0,0,0,0,0,2,0],[0,0,0,0,0,0,2]]"
JSON.stringify(dvec)
"[0,-10,-10,-30,-20,-30,-20]"
JSON.stringify(Amat)
"[[],[0,1,-1,0,0,0,0,-1,0,0,1,0,0,0,0,0],[0,0,0,1,-1,0,0,-1,0,0,0,1,0,0,0,0],[0,0,0,1,-1,0,0,0,-1,0,0,0,1,0,0,0],[0,0,0,1,-1,0,0,0,0,-1,0,0,0,1,0,0],[0,0,0,0,0,1,-1,0,-1,0,0,0,0,0,1,0],[0,0,0,0,0,1,-1,0,0,-1,0,0,0,0,0,1]]"
JSON.stringify(bvec)
"[0,3.999,-4.001,4.999,-5.001,0.999,-1.001,-7.001,-1.001,-2.001,-0.001,-0.001,-0.001,-0.001,-0.001,-0.001]"
JSON.stringify(meq)
"0"
JSON.stringify(factorized)
"[0,0]"
This returns:
factorized = [0, 0]
output = {"solution":[,NaN,NaN,NaN,NaN,NaN,NaN],"Lagrangian":[,NaN,0,NaN,0,NaN,0,0,0,0,0,0,0,0,NaN,0],"value":[,NaN],"unconstrained_solution":[0,-4.999999999999999,-4.999999999999999,-14.999999999999996,-9.999999999999998,-14.999999999999996,-9.999999999999998],"iterations":[,6,1],"iact":[,3,5,1,14,0,0,0,0,0,0,0,0,0,0,0],"message":""}

If there is more than one possible solution with equal minimum value, it's enough to get any solution.

The constraints are of three types:

1. x2 + x3 + x4 = 5,
   this returned NaNs for many inputs, so I substituted it with two inequalities:
   x2 + x3 + x4 >= 4.999 and
   x2 + x3 + x4 <= 5.001
   This substitution helped a lot.
2. sum of some variables <= some positive integer (eg. 7 as in bvec[7])
   this I substituted with:
   sum of some variables <= some positive integer + epsilon
   eg.
   sum of some variables <= 7.001
3. each variable >= 0
   that I substituted with:
   each variable >= -epsilon
   eg.
   each variable >= -0.001

The best I would like to get only integers in solution but I'm noticed you wrote the solver that returns reals.