I'm solving a constrained optimization problem, where I have to find the vector $\mathbf{x}$ of dimensionality N x 1. The input vectors are $\mathbf{a}$ and $\mathbf{h}$ with dimensionality N x 1. The error function which I have to minimize is:
$$E(\mathbf{a},\mathbf{x})=\sum_{i=1}^N(a_i-x_i)^2=(\mathbf{a}-\mathbf{x})^T(\mathbf{a}-\mathbf{x}) = \mathbf{x}^T\mathbf{x} - 2\mathbf{a}^T\mathbf{x} + \mathbf{a}^T\mathbf{a}$$
The term $\mathbf{a}^T\mathbf{a}$ can be ignored, as it is just a constant term, and doesn't influence the choice of $\mathbf{x}$, so $E(\mathbf{a},\mathbf{x})=\mathbf{x}^T\mathbf{x} - 2\mathbf{a}^T\mathbf{x}$.
I've to minimize this error function subject to $N$ constraints:
$$x_i \ge \varepsilon_i \; \; \; \; \; \; \forall i \in [1, N] $$
where $\varepsilon_i$ is defined as:
$$\varepsilon_i =\begin{cases} k &for \; i = 1\\ x_{i-1} + h_{i - 1} + \eta &for \; i \gt 1 \end{cases}$$
for some constants $k$ and $\eta$.
So, its a quadratic programming problem. I've to use a QP solver to solve this problem in JavaScript, so I've found numeric.js which provides a function solveQP (The Quadratic Programming function numeric.solveQP() is based on Alberto Santini's quadprog, which is itself a port of the R package quadprog.)
As R package quadprog solves a QP problem of the form: $min(\frac{1}{2} x^T D x - d^T x)$ s.t. inequality constraints $A^Tx \geq b_0$
So, to convert my problem into this problem, I get, $D = I$ (or $D = 2I$) and $d = 2\mathbf{a}$. And my constraints are like:
$$x_1 \ge k$$ $$x_2 \ge x_1 + h_1 + \eta => x_2 - x_1 \ge h_1 + \eta$$ $$x_3 - x_2 \ge h_2 + \eta$$ $$\vdots$$
And so on. So I get: $A = \begin{pmatrix}1 & 0 & 0 \ldots \\-1 & 1 & 0 \ldots\\0 & -1 & 1 \ldots\\\vdots \end{pmatrix} $, and $b_0 = \begin{pmatrix}k\\h_1+\eta\\h_2+\eta\\\vdots\end{pmatrix}$
But when I code this, and call the solveQP method, the solution doesn't even satisfies the constraints.
E.g. for $N=3$, $\eta=2$, $k=5$, $h=\begin{pmatrix}2\\5\\5\end{pmatrix}$, and $a=\begin{pmatrix}8\\15\\16\end{pmatrix}$, it gives the solution as: $\begin{pmatrix}25.5\\28.25\\24.25\end{pmatrix}$, which doesn't satisfy the constraints because $x_2=28.25$ must be greater than or equal to $x_1+h_1+\eta = 25.5+2+2=29.5$ and $x_3=24.25$ must be greater than or equal to $x_2+h_2+\eta = 28.25+5+2=35.25$ (considering the current value of $x_2$). This solution is, if I use $D=I$. And if I use $D=2I$ (which I think, should not make any difference), the solution comes out to be: $\begin{pmatrix}13.999999999999996\\14.499999999999998\\10.499999999999998\end{pmatrix}=\begin{pmatrix}14\\14.5\\10.5\end{pmatrix}$, which again, doesn't even satisfy the constraints.
So, does anyone know a better QP solver for JavaScript? Or am I making any mistake that I'm not able to figure out?
These are the functions from numeric.js that my solution needs:
// functions from numeric.js that my solution needs
function base0to1(A) {
if(typeof A !== "object") { return A; }
var ret = [], i,n=A.length;
for(i=0;i<n;i++) ret[i+1] = base0to1(A[i]);
return ret;
}
function base1to0(A) {
if(typeof A !== "object") { return A; }
var ret = [], i,n=A.length;
for(i=1;i<n;i++) ret[i-1] = base1to0(A[i]);
return ret;
}
function dpofa(a, lda, n, info) {
var i, j, jm1, k, t, s;
for (j = 1; j <= n; j = j + 1) {
info[1] = j;
s = 0;
jm1 = j - 1;
if (jm1 < 1) {
s = a[j][j] - s;
if (s <= 0) {
break;
}
a[j][j] = Math.sqrt(s);
} else {
for (k = 1; k <= jm1; k = k + 1) {
//~ t = a[k][j] - ddot(k - 1, a[1][k], 1, a[1][j], 1);
t = a[k][j];
for (i = 1; i < k; i = i + 1) {
t = t - (a[i][j] * a[i][k]);
}
t = t / a[k][k];
a[k][j] = t;
s = s + t * t;
}
s = a[j][j] - s;
if (s <= 0) {
break;
}
a[j][j] = Math.sqrt(s);
}
info[1] = 0;
}
}
function dposl(a, lda, n, b) {
var i, k, kb, t;
for (k = 1; k <= n; k = k + 1) {
//~ t = ddot(k - 1, a[1][k], 1, b[1], 1);
t = 0;
for (i = 1; i < k; i = i + 1) {
t = t + (a[i][k] * b[i]);
}
b[k] = (b[k] - t) / a[k][k];
}
for (kb = 1; kb <= n; kb = kb + 1) {
k = n + 1 - kb;
b[k] = b[k] / a[k][k];
t = -b[k];
//~ daxpy(k - 1, t, a[1][k], 1, b[1], 1);
for (i = 1; i < k; i = i + 1) {
b[i] = b[i] + (t * a[i][k]);
}
}
}
function dpori(a, lda, n) {
var i, j, k, kp1, t;
for (k = 1; k <= n; k = k + 1) {
a[k][k] = 1 / a[k][k];
t = -a[k][k];
//~ dscal(k - 1, t, a[1][k], 1);
for (i = 1; i < k; i = i + 1) {
a[i][k] = t * a[i][k];
}
kp1 = k + 1;
if (n < kp1) {
break;
}
for (j = kp1; j <= n; j = j + 1) {
t = a[k][j];
a[k][j] = 0;
//~ daxpy(k, t, a[1][k], 1, a[1][j], 1);
for (i = 1; i <= k; i = i + 1) {
a[i][j] = a[i][j] + (t * a[i][k]);
}
}
}
}
function qpgen2(dmat, dvec, fddmat, n, sol, crval, amat,
bvec, fdamat, q, meq, iact, nact, iter, work, ierr) {
var i, j, l, l1, info, it1, iwzv, iwrv, iwrm, iwsv, iwuv, nvl, r, iwnbv,
temp, sum, t1, tt, gc, gs, nu,
t1inf, t2min,
vsmall, tmpa, tmpb,
go;
r = Math.min(n, q);
l = 2 * n + (r * (r + 5)) / 2 + 2 * q + 1;
vsmall = 1.0e-60;
do {
vsmall = vsmall + vsmall;
tmpa = 1 + 0.1 * vsmall;
tmpb = 1 + 0.2 * vsmall;
} while (tmpa <= 1 || tmpb <= 1);
for (i = 1; i <= n; i = i + 1) {
work[i] = dvec[i];
}
for (i = n + 1; i <= l; i = i + 1) {
work[i] = 0;
}
for (i = 1; i <= q; i = i + 1) {
iact[i] = 0;
}
info = [];
if (ierr[1] === 0) {
dpofa(dmat, fddmat, n, info);
if (info[1] !== 0) {
ierr[1] = 2;
return;
}
dposl(dmat, fddmat, n, dvec);
dpori(dmat, fddmat, n);
} else {
for (j = 1; j <= n; j = j + 1) {
sol[j] = 0;
for (i = 1; i <= j; i = i + 1) {
sol[j] = sol[j] + dmat[i][j] * dvec[i];
}
}
for (j = 1; j <= n; j = j + 1) {
dvec[j] = 0;
for (i = j; i <= n; i = i + 1) {
dvec[j] = dvec[j] + dmat[j][i] * sol[i];
}
}
}
crval[1] = 0;
for (j = 1; j <= n; j = j + 1) {
sol[j] = dvec[j];
crval[1] = crval[1] + work[j] * sol[j];
work[j] = 0;
for (i = j + 1; i <= n; i = i + 1) {
dmat[i][j] = 0;
}
}
crval[1] = -crval[1] / 2;
ierr[1] = 0;
iwzv = n;
iwrv = iwzv + n;
iwuv = iwrv + r;
iwrm = iwuv + r + 1;
iwsv = iwrm + (r * (r + 1)) / 2;
iwnbv = iwsv + q;
for (i = 1; i <= q; i = i + 1) {
sum = 0;
for (j = 1; j <= n; j = j + 1) {
sum = sum + amat[j][i] * amat[j][i];
}
work[iwnbv + i] = Math.sqrt(sum);
}
nact = 0;
iter[1] = 0;
iter[2] = 0;
function fn_goto_50() {
iter[1] = iter[1] + 1;
l = iwsv;
for (i = 1; i <= q; i = i + 1) {
l = l + 1;
sum = -bvec[i];
for (j = 1; j <= n; j = j + 1) {
sum = sum + amat[j][i] * sol[j];
}
if (Math.abs(sum) < vsmall) {
sum = 0;
}
if (i > meq) {
work[l] = sum;
} else {
work[l] = -Math.abs(sum);
if (sum > 0) {
for (j = 1; j <= n; j = j + 1) {
amat[j][i] = -amat[j][i];
}
bvec[i] = -bvec[i];
}
}
}
for (i = 1; i <= nact; i = i + 1) {
work[iwsv + iact[i]] = 0;
}
nvl = 0;
temp = 0;
for (i = 1; i <= q; i = i + 1) {
if (work[iwsv + i] < temp * work[iwnbv + i]) {
nvl = i;
temp = work[iwsv + i] / work[iwnbv + i];
}
}
if (nvl === 0) {
return 999;
}
return 0;
}
function fn_goto_55() {
for (i = 1; i <= n; i = i + 1) {
sum = 0;
for (j = 1; j <= n; j = j + 1) {
sum = sum + dmat[j][i] * amat[j][nvl];
}
work[i] = sum;
}
l1 = iwzv;
for (i = 1; i <= n; i = i + 1) {
work[l1 + i] = 0;
}
for (j = nact + 1; j <= n; j = j + 1) {
for (i = 1; i <= n; i = i + 1) {
work[l1 + i] = work[l1 + i] + dmat[i][j] * work[j];
}
}
t1inf = true;
for (i = nact; i >= 1; i = i - 1) {
sum = work[i];
l = iwrm + (i * (i + 3)) / 2;
l1 = l - i;
for (j = i + 1; j <= nact; j = j + 1) {
sum = sum - work[l] * work[iwrv + j];
l = l + j;
}
sum = sum / work[l1];
work[iwrv + i] = sum;
if (iact[i] < meq) {
// continue;
break;
}
if (sum < 0) {
// continue;
break;
}
t1inf = false;
it1 = i;
}
if (!t1inf) {
t1 = work[iwuv + it1] / work[iwrv + it1];
for (i = 1; i <= nact; i = i + 1) {
if (iact[i] < meq) {
// continue;
break;
}
if (work[iwrv + i] < 0) {
// continue;
break;
}
temp = work[iwuv + i] / work[iwrv + i];
if (temp < t1) {
t1 = temp;
it1 = i;
}
}
}
sum = 0;
for (i = iwzv + 1; i <= iwzv + n; i = i + 1) {
sum = sum + work[i] * work[i];
}
if (Math.abs(sum) <= vsmall) {
if (t1inf) {
ierr[1] = 1;
// GOTO 999
return 999;
} else {
for (i = 1; i <= nact; i = i + 1) {
work[iwuv + i] = work[iwuv + i] - t1 * work[iwrv + i];
}
work[iwuv + nact + 1] = work[iwuv + nact + 1] + t1;
// GOTO 700
return 700;
}
} else {
sum = 0;
for (i = 1; i <= n; i = i + 1) {
sum = sum + work[iwzv + i] * amat[i][nvl];
}
tt = -work[iwsv + nvl] / sum;
t2min = true;
if (!t1inf) {
if (t1 < tt) {
tt = t1;
t2min = false;
}
}
for (i = 1; i <= n; i = i + 1) {
sol[i] = sol[i] + tt * work[iwzv + i];
if (Math.abs(sol[i]) < vsmall) {
sol[i] = 0;
}
}
crval[1] = crval[1] + tt * sum * (tt / 2 + work[iwuv + nact + 1]);
for (i = 1; i <= nact; i = i + 1) {
work[iwuv + i] = work[iwuv + i] - tt * work[iwrv + i];
}
work[iwuv + nact + 1] = work[iwuv + nact + 1] + tt;
if (t2min) {
nact = nact + 1;
iact[nact] = nvl;
l = iwrm + ((nact - 1) * nact) / 2 + 1;
for (i = 1; i <= nact - 1; i = i + 1) {
work[l] = work[i];
l = l + 1;
}
if (nact === n) {
work[l] = work[n];
} else {
for (i = n; i >= nact + 1; i = i - 1) {
if (work[i] === 0) {
// continue;
break;
}
gc = Math.max(Math.abs(work[i - 1]), Math.abs(work[i]));
gs = Math.min(Math.abs(work[i - 1]), Math.abs(work[i]));
if (work[i - 1] >= 0) {
temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
} else {
temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
}
gc = work[i - 1] / temp;
gs = work[i] / temp;
if (gc === 1) {
// continue;
break;
}
if (gc === 0) {
work[i - 1] = gs * temp;
for (j = 1; j <= n; j = j + 1) {
temp = dmat[j][i - 1];
dmat[j][i - 1] = dmat[j][i];
dmat[j][i] = temp;
}
} else {
work[i - 1] = temp;
nu = gs / (1 + gc);
for (j = 1; j <= n; j = j + 1) {
temp = gc * dmat[j][i - 1] + gs * dmat[j][i];
dmat[j][i] = nu * (dmat[j][i - 1] + temp) - dmat[j][i];
dmat[j][i - 1] = temp;
}
}
}
work[l] = work[nact];
}
} else {
sum = -bvec[nvl];
for (j = 1; j <= n; j = j + 1) {
sum = sum + sol[j] * amat[j][nvl];
}
if (nvl > meq) {
work[iwsv + nvl] = sum;
} else {
work[iwsv + nvl] = -Math.abs(sum);
if (sum > 0) {
for (j = 1; j <= n; j = j + 1) {
amat[j][nvl] = -amat[j][nvl];
}
bvec[nvl] = -bvec[nvl];
}
}
// GOTO 700
return 700;
}
}
return 0;
}
function fn_goto_797() {
l = iwrm + (it1 * (it1 + 1)) / 2 + 1;
l1 = l + it1;
if (work[l1] === 0) {
// GOTO 798
return 798;
}
gc = Math.max(Math.abs(work[l1 - 1]), Math.abs(work[l1]));
gs = Math.min(Math.abs(work[l1 - 1]), Math.abs(work[l1]));
if (work[l1 - 1] >= 0) {
temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
} else {
temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
}
gc = work[l1 - 1] / temp;
gs = work[l1] / temp;
if (gc === 1) {
// GOTO 798
return 798;
}
if (gc === 0) {
for (i = it1 + 1; i <= nact; i = i + 1) {
temp = work[l1 - 1];
work[l1 - 1] = work[l1];
work[l1] = temp;
l1 = l1 + i;
}
for (i = 1; i <= n; i = i + 1) {
temp = dmat[i][it1];
dmat[i][it1] = dmat[i][it1 + 1];
dmat[i][it1 + 1] = temp;
}
} else {
nu = gs / (1 + gc);
for (i = it1 + 1; i <= nact; i = i + 1) {
temp = gc * work[l1 - 1] + gs * work[l1];
work[l1] = nu * (work[l1 - 1] + temp) - work[l1];
work[l1 - 1] = temp;
l1 = l1 + i;
}
for (i = 1; i <= n; i = i + 1) {
temp = gc * dmat[i][it1] + gs * dmat[i][it1 + 1];
dmat[i][it1 + 1] = nu * (dmat[i][it1] + temp) - dmat[i][it1 + 1];
dmat[i][it1] = temp;
}
}
return 0;
}
function fn_goto_798() {
l1 = l - it1;
for (i = 1; i <= it1; i = i + 1) {
work[l1] = work[l];
l = l + 1;
l1 = l1 + 1;
}
work[iwuv + it1] = work[iwuv + it1 + 1];
iact[it1] = iact[it1 + 1];
it1 = it1 + 1;
if (it1 < nact) {
// GOTO 797
return 797;
}
return 0;
}
function fn_goto_799() {
work[iwuv + nact] = work[iwuv + nact + 1];
work[iwuv + nact + 1] = 0;
iact[nact] = 0;
nact = nact - 1;
iter[2] = iter[2] + 1;
return 0;
}
go = 0;
while (true) {
go = fn_goto_50();
if (go === 999) {
return;
}
while (true) {
go = fn_goto_55();
if (go === 0) {
break;
}
if (go === 999) {
return;
}
if (go === 700) {
if (it1 === nact) {
fn_goto_799();
} else {
while (true) {
fn_goto_797();
go = fn_goto_798();
if (go !== 797) {
break;
}
}
fn_goto_799();
}
}
}
}
}
function solveQP(Dmat, dvec, Amat, bvec, meq, factorized) {
Dmat = base0to1(Dmat);
dvec = base0to1(dvec);
Amat = base0to1(Amat);
var i, n, q,
nact, r,
crval = [], iact = [], sol = [], work = [], iter = [],
message;
meq = meq || 0;
factorized = factorized ? base0to1(factorized) : [undefined, 0];
bvec = bvec ? base0to1(bvec) : [];
// In Fortran the array index starts from 1
n = Dmat.length - 1;
q = Amat[1].length - 1;
if (!bvec) {
for (i = 1; i <= q; i = i + 1) {
bvec[i] = 0;
}
}
for (i = 1; i <= q; i = i + 1) {
iact[i] = 0;
}
nact = 0;
r = Math.min(n, q);
for (i = 1; i <= n; i = i + 1) {
sol[i] = 0;
}
crval[1] = 0;
for (i = 1; i <= (2 * n + (r * (r + 5)) / 2 + 2 * q + 1); i = i + 1) {
work[i] = 0;
}
for (i = 1; i <= 2; i = i + 1) {
iter[i] = 0;
}
qpgen2(Dmat, dvec, n, n, sol, crval, Amat,
bvec, n, q, meq, iact, nact, iter, work, factorized);
message = "";
if (factorized[1] === 1) {
message = "constraints are inconsistent, no solution!";
}
if (factorized[1] === 2) {
message = "matrix D in quadratic function is not positive definite!";
}
return {
solution: base1to0(sol),
value: base1to0(crval),
unconstrained_solution: base1to0(dvec),
iterations: base1to0(iter),
iact: base1to0(iact),
message: message
};
}
My code that constructs $D$, $d$, $A$, and $b_0$, and calls the solveQP function:
// main
// solve QP with 1/2 (x^T) D (x) - (d^T)(x) s.t. (A^T)x >= b
var h = [2, 5, 5]; // vector of heights
var a = [8, 15, 16]; // input vector
var N = h.length;
var eta = 2;
var k = 5;
// print h
document.write("h...\n");
for (i = 0; i < N; i++){
document.write(h[i] + ", ");
}
// print a
document.write("\na...\n");
for (i = 0; i < N; i++){
document.write(a[i] + ", ");
}
// print N
document.write("\nN..." + N + "\n");
// print eta
document.write("\neta..." + eta + "\n");
// print k
document.write("\nk..." + k + "\n");
// construct DMat = 2I, and dvec = 2a
var DMat = [];
for (i = 0; i < N; i++){
DMat[i] = [];
for (j = 0; j < N; j++){
if (i != j){
DMat[i][j] = 0;
}
else{
DMat[i][j] = 2;
}
}
}
// print DMat
document.write("\nDMat...\n");
for (i = 0; i < N; i++){
for (j = 0; j < N; j++){
document.write(DMat[i][j] + ", ");
}
document.write("\n");
}
// now dvec
var dvec = [];
for (i = 0; i < N; i++){
dvec[i] = a[i] * 2;
}
// print dvec
document.write("\ndvec...\n");
for (i = 0; i < N; i++){
document.write(dvec[i] + ", ");
}
// construct AMat and bvec for constraints
var AMat = [];
for (i = 0; i < N; i++){
AMat[i] = [];
for (j = 0; j < N; j++){
if (i == j){
AMat[i][j] = 1;
}
else if ((j+1) == i){
AMat[i][j] = -1;
}
else{
AMat[i][j] = 0;
}
}
}
// print AMat
document.write("\nAMat...\n");
for (i = 0; i < N; i++){
for (j = 0; j < N; j++){
document.write(AMat[i][j] + ", ");
}
document.write("\n");
}
// now bvec
var bvec = [];
for (i = 0; i < N; i++){
if (i == 0){
bvec[i] = k;
}
else{
bvec[i] = h[i - 1] + eta;
}
}
// print bvec
document.write("\nbvec...\n");
for (i = 0; i < N; i++){
document.write(bvec[i] + ", ");
}
// solve QP with 1/2 (xT) D (x) - (dT)(x) s.t. Ax >= b
var res = solveQP(DMat, dvec, AMat, bvec, 0, false);
document.write("\nSolution...\n");
for (i = 0; i < res.solution.length; i++){
document.write(res.solution[i] + ", ");
}
Looks like there was some error in
numeric.js, so I copied the code from Alberto Santini's quadprog on whichnumeric.js'ssolveQPwas based on. Below, it is given:This code uses 1-based indexed arrays instead of standard JS 0-based indexed arrays (So you can leave out the index 0, and only store the values starting from 1).
And I was making one mistake too, which was: the parameter
AMatinsolveQPshould be $A^T$, not $A$. So anyone else looking for JavaScript QP solver, this one is working for me, you can use it too.