# Infeasible Primal-Dual Quadratic Programming with Box Constraints Jan de Leeuw Version 08, May 09, 2017 Note: This is a working paper which will be expanded/updated frequently. All suggestions for improvement are welcome. The directory [gifi.stat.ucla.edu/boxqp](http://gifi.stat.ucla.edu/boxqp) has a pdf version, the complete Rmd file with all code chunks, the R and C code, and the bib file. #Introduction The problem we address in this paper is minimization of $$\label{E:f} f(x)=\frac12 x'Cx+d'x$$ over $a\leq x\leq b$ (elementwise). Here $C$ is positive definite, so $f$ is strictly convex. Also we assume without loss of generality that \$ab_i\text{ or }(x_i^{(k)}=b_i\text{ and }\mu_i^{(k)}\geq 0)\}\\ S^{(k)}&=\{i\mid a_i #include inline int VINDEX(const int); inline int MINDEX(const int, const int, const int); inline int VINDEX(const int i) { return i - 1; } inline int MINDEX(const int i, const int j, const int n) { return (i - 1) + (j - 1) * n; } void boxcqp(const double *, const double *, double *, double *, double *, const double *, const double *, double *, double *, const int *, int *, int *, double *, const int *); void in(const double *, const double *, const double *, const double *, const double *, const int *, int *, int *); void sb(double *x, double *lambda, double *mu, const double *lower, const double *upper, const int *ind, const int *n); void ss(const double *, const double *, double *, double *, double *, const double *, const double *, double *, double *, const int *, const int *, const int *); void co(const double *, const double *, const double *, const double *, const double *, const int *, const int *, bool *); void ff(const double *, const double *, const double *, const int *, double *); void dposv(const int *, const int *, double *, double *, int *); void mprint(const double *, const int *, const int *); void check(const double *, const double *, const double *, const double *, const double *, const double *, const double *, const int *, const int *); void dposv(const int *n, const int *nrhs, double *a, double *b, int *status) { lapack_int info = 0, nn = (lapack_int)*n, nr = (lapack_int)*nrhs; info = LAPACKE_dposv(LAPACK_COL_MAJOR, 'U', nn, nr, a, nn, b, nn); *status = info; return; } void boxcqp(const double *a, const double *b, double *x, double *lambda, double *mu, const double *lower, const double *upper, double *amat, double *bvec, const int *n, int *ind, int *itel, double *f, const int *itmax) { int nzero = 0, ktel = 1, isit = 0; bool *opt = (bool *)&isit; while (true) { (void)in(x, lambda, mu, lower, upper, n, ind, &nzero); (void)sb(x, lambda, mu, lower, upper, n, ind); (void)ss(a, b, x, lambda, mu, lower, upper, amat, bvec, n, ind, &nzero); (void)co(x, lambda, mu, lower, upper, n, ind, opt); if ((ktel == *itmax) || (*opt == true)) { break; } ktel++; } *itel = ktel; (void)ff(a, b, x, n, f); return; } void in(const double *x, const double *lambda, const double *mu, const double *lower, const double *upper, const int *n, int *ind, int *nzero) { int nvar = *n, nz = 0; for (int i = 1; i <= nvar; i++) { int ii = VINDEX(i); if (x[ii] < lower[ii]) { ind[ii] = -1; } if (x[ii] > upper[ii]) { ind[ii] = 1; } if ((x[ii] > lower[ii]) && (x[ii] < upper[ii])) { ind[ii] = 0; nz++; } if (x[ii] == lower[ii]) { if (lambda[ii] >= 0) { ind[ii] = -1; } else { ind[ii] = 0; nz++; } } if (x[ii] == upper[ii]) { if (mu[ii] >= 0) { ind[ii] = 1; } else { ind[ii] = 0; nz++; } } } *nzero = nz; return; } void sb(double *x, double *lambda, double *mu, const double *lower, const double *upper, const int *n, const int *ind) { int nvar = *n; for (int i = 1; i <= nvar; i++) { int ii = VINDEX(i); if (ind[ii] == -1) { x[ii] = lower[ii]; mu[ii] = 0.0; } if (ind[ii] == 1) { x[ii] = upper[ii]; lambda[ii] = 0.0; } if (ind[ii] == 0) { lambda[ii] = 0.0; mu[ii] = 0.0; } } return; } void ss(const double *a, const double *b, double *x, double *lambda, double *mu, const double *lower, const double *upper, double *amat, double *bvec, const int *n, const int *ind, const int *nzero) { int nvar = *n, nz = *nzero, status = 0, k = 1, nr = 1; for (int i = 1; i <= nvar; i++) { int ii = VINDEX(i); if (ind[ii] == 0) { int kk = VINDEX(k); bvec[kk] = -b[ii]; int l = 1; for (int j = 1; j <= nvar; j++) { int jj = VINDEX(j); double aa = a[MINDEX(i, j, nvar)]; if (ind[jj] == -1) { bvec[kk] -= aa * lower[jj]; } if (ind[jj] == 1) { bvec[kk] -= aa * upper[jj]; } if (ind[jj] == 0) { amat[MINDEX(k, l, nz)] = aa; l++; } } k++; } } (void)dposv(nzero, &nr, amat, bvec, &status); k = 0; for (int i = 1; i <= nvar; i++) { int ii = VINDEX(i); if (ind[ii] == 0) { x[ii] = bvec[k]; k++; } } for (int i = 1; i <= nvar; i++) { int ii = VINDEX(i); if (ind[ii] == -1) { lambda[ii] = b[ii]; for (int j = 1; j <= nvar; j++) { int jj = VINDEX(j); double aa = a[MINDEX(i, j, nvar)]; lambda[ii] += aa * x[jj]; } } if (ind[ii] == 1) { mu[ii] = -b[ii]; for (int j = 1; j <= nvar; j++) { int jj = VINDEX(j); double aa = a[MINDEX(i, j, nvar)]; mu[ii] -= aa * x[jj]; } } } return; } void co(const double *x, const double *lambda, const double *mu, const double *lower, const double *upper, const int *n, const int *ind, bool *opt) { int nvar = *n; bool isit = true; for (int i = 1; i <= nvar; i++) { int ii = VINDEX(i); if (ind[ii] == 0) { if (x[ii] < lower[ii]) { isit = false; } if (x[ii] > upper[ii]) { isit = false; } } if (ind[ii] == 1) { if (mu[ii] < 0) { isit = false; } } if (ind[ii] == -1) { if (lambda[ii] < 0) { isit = false; } } } *opt = isit; return; } void ff(const double *a, const double *b, const double *x, const int *n, double *f) { int nvar = *n; double sa = 0.0, sb = 0.0; for (int i = 1; i <= nvar; i++) { int ii = VINDEX(i); sb += b[ii] * x[ii]; for (int j = 1; j <= nvar; j++) { int jj = VINDEX(j); double aa = a[MINDEX(i, j, nvar)]; sa += aa * x[ii] * x[jj]; } } *f = sa / 2.0 + sb; } ` #References