Actual source code: bqpip.c
1: /*
2: This file implements a Mehrotra predictor-corrector method for
3: bound-constrained quadratic programs.
5: */
7: #include <../src/tao/quadratic/impls/bqpip/bqpipimpl.h>
8: #include <petscksp.h>
10: static PetscErrorCode QPIPComputeResidual(TAO_BQPIP *qp, Tao tao)
11: {
12: PetscReal dtmp = 1.0 - qp->psteplength;
14: /* Compute R3 and R5 */
16: VecScale(qp->R3, dtmp);
17: VecScale(qp->R5, dtmp);
18: qp->pinfeas = dtmp * qp->pinfeas;
20: VecCopy(qp->S, tao->gradient);
21: VecAXPY(tao->gradient, -1.0, qp->Z);
23: MatMult(tao->hessian, tao->solution, qp->RHS);
24: VecScale(qp->RHS, -1.0);
25: VecAXPY(qp->RHS, -1.0, qp->C);
26: VecAXPY(tao->gradient, -1.0, qp->RHS);
28: VecNorm(tao->gradient, NORM_1, &qp->dinfeas);
29: qp->rnorm = (qp->dinfeas + qp->pinfeas) / (qp->m + qp->n);
30: return 0;
31: }
33: static PetscErrorCode QPIPSetInitialPoint(TAO_BQPIP *qp, Tao tao)
34: {
35: PetscReal two = 2.0, p01 = 1;
36: PetscReal gap1, gap2, fff, mu;
38: /* Compute function, Gradient R=Hx+b, and Hessian */
39: MatMult(tao->hessian, tao->solution, tao->gradient);
40: VecCopy(qp->C, qp->Work);
41: VecAXPY(qp->Work, 0.5, tao->gradient);
42: VecAXPY(tao->gradient, 1.0, qp->C);
43: VecDot(tao->solution, qp->Work, &fff);
44: qp->pobj = fff + qp->d;
48: /* Initialize slack vectors */
49: /* T = XU - X; G = X - XL */
50: VecCopy(qp->XU, qp->T);
51: VecAXPY(qp->T, -1.0, tao->solution);
52: VecCopy(tao->solution, qp->G);
53: VecAXPY(qp->G, -1.0, qp->XL);
55: VecSet(qp->GZwork, p01);
56: VecSet(qp->TSwork, p01);
58: VecPointwiseMax(qp->G, qp->G, qp->GZwork);
59: VecPointwiseMax(qp->T, qp->T, qp->TSwork);
61: /* Initialize Dual Variable Vectors */
62: VecCopy(qp->G, qp->Z);
63: VecReciprocal(qp->Z);
65: VecCopy(qp->T, qp->S);
66: VecReciprocal(qp->S);
68: MatMult(tao->hessian, qp->Work, qp->RHS);
69: VecAbs(qp->RHS);
70: VecSet(qp->Work, p01);
71: VecPointwiseMax(qp->RHS, qp->RHS, qp->Work);
73: VecPointwiseDivide(qp->RHS, tao->gradient, qp->RHS);
74: VecNorm(qp->RHS, NORM_1, &gap1);
75: mu = PetscMin(10.0, (gap1 + 10.0) / qp->m);
77: VecScale(qp->S, mu);
78: VecScale(qp->Z, mu);
80: VecSet(qp->TSwork, p01);
81: VecSet(qp->GZwork, p01);
82: VecPointwiseMax(qp->S, qp->S, qp->TSwork);
83: VecPointwiseMax(qp->Z, qp->Z, qp->GZwork);
85: qp->mu = 0;
86: qp->dinfeas = 1.0;
87: qp->pinfeas = 1.0;
88: while ((qp->dinfeas + qp->pinfeas) / (qp->m + qp->n) >= qp->mu) {
89: VecScale(qp->G, two);
90: VecScale(qp->Z, two);
91: VecScale(qp->S, two);
92: VecScale(qp->T, two);
94: QPIPComputeResidual(qp, tao);
96: VecCopy(tao->solution, qp->R3);
97: VecAXPY(qp->R3, -1.0, qp->G);
98: VecAXPY(qp->R3, -1.0, qp->XL);
100: VecCopy(tao->solution, qp->R5);
101: VecAXPY(qp->R5, 1.0, qp->T);
102: VecAXPY(qp->R5, -1.0, qp->XU);
104: VecNorm(qp->R3, NORM_INFINITY, &gap1);
105: VecNorm(qp->R5, NORM_INFINITY, &gap2);
106: qp->pinfeas = PetscMax(gap1, gap2);
108: /* Compute the duality gap */
109: VecDot(qp->G, qp->Z, &gap1);
110: VecDot(qp->T, qp->S, &gap2);
112: qp->gap = gap1 + gap2;
113: qp->dobj = qp->pobj - qp->gap;
114: if (qp->m > 0) {
115: qp->mu = qp->gap / (qp->m);
116: } else {
117: qp->mu = 0.0;
118: }
119: qp->rgap = qp->gap / (PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0);
120: }
121: return 0;
122: }
124: static PetscErrorCode QPIPStepLength(TAO_BQPIP *qp)
125: {
126: PetscReal tstep1, tstep2, tstep3, tstep4, tstep;
128: /* Compute stepsize to the boundary */
129: VecStepMax(qp->G, qp->DG, &tstep1);
130: VecStepMax(qp->T, qp->DT, &tstep2);
131: VecStepMax(qp->S, qp->DS, &tstep3);
132: VecStepMax(qp->Z, qp->DZ, &tstep4);
134: tstep = PetscMin(tstep1, tstep2);
135: qp->psteplength = PetscMin(0.95 * tstep, 1.0);
137: tstep = PetscMin(tstep3, tstep4);
138: qp->dsteplength = PetscMin(0.95 * tstep, 1.0);
140: qp->psteplength = PetscMin(qp->psteplength, qp->dsteplength);
141: qp->dsteplength = qp->psteplength;
142: return 0;
143: }
145: static PetscErrorCode QPIPComputeNormFromCentralPath(TAO_BQPIP *qp, PetscReal *norm)
146: {
147: PetscReal gap[2], mu[2], nmu;
149: VecPointwiseMult(qp->GZwork, qp->G, qp->Z);
150: VecPointwiseMult(qp->TSwork, qp->T, qp->S);
151: VecNorm(qp->TSwork, NORM_1, &mu[0]);
152: VecNorm(qp->GZwork, NORM_1, &mu[1]);
154: nmu = -(mu[0] + mu[1]) / qp->m;
156: VecShift(qp->GZwork, nmu);
157: VecShift(qp->TSwork, nmu);
159: VecNorm(qp->GZwork, NORM_2, &gap[0]);
160: VecNorm(qp->TSwork, NORM_2, &gap[1]);
161: gap[0] *= gap[0];
162: gap[1] *= gap[1];
164: qp->pathnorm = PetscSqrtScalar(gap[0] + gap[1]);
165: *norm = qp->pathnorm;
166: return 0;
167: }
169: static PetscErrorCode QPIPComputeStepDirection(TAO_BQPIP *qp, Tao tao)
170: {
171: /* Calculate DG */
172: VecCopy(tao->stepdirection, qp->DG);
173: VecAXPY(qp->DG, 1.0, qp->R3);
175: /* Calculate DT */
176: VecCopy(tao->stepdirection, qp->DT);
177: VecScale(qp->DT, -1.0);
178: VecAXPY(qp->DT, -1.0, qp->R5);
180: /* Calculate DZ */
181: VecAXPY(qp->DZ, -1.0, qp->Z);
182: VecPointwiseDivide(qp->GZwork, qp->DG, qp->G);
183: VecPointwiseMult(qp->GZwork, qp->GZwork, qp->Z);
184: VecAXPY(qp->DZ, -1.0, qp->GZwork);
186: /* Calculate DS */
187: VecAXPY(qp->DS, -1.0, qp->S);
188: VecPointwiseDivide(qp->TSwork, qp->DT, qp->T);
189: VecPointwiseMult(qp->TSwork, qp->TSwork, qp->S);
190: VecAXPY(qp->DS, -1.0, qp->TSwork);
191: return 0;
192: }
194: static PetscErrorCode TaoSetup_BQPIP(Tao tao)
195: {
196: TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
198: /* Set pointers to Data */
199: VecGetSize(tao->solution, &qp->n);
201: /* Allocate some arrays */
202: if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
203: if (!tao->stepdirection) VecDuplicate(tao->solution, &tao->stepdirection);
205: VecDuplicate(tao->solution, &qp->Work);
206: VecDuplicate(tao->solution, &qp->XU);
207: VecDuplicate(tao->solution, &qp->XL);
208: VecDuplicate(tao->solution, &qp->HDiag);
209: VecDuplicate(tao->solution, &qp->DiagAxpy);
210: VecDuplicate(tao->solution, &qp->RHS);
211: VecDuplicate(tao->solution, &qp->RHS2);
212: VecDuplicate(tao->solution, &qp->C);
214: VecDuplicate(tao->solution, &qp->G);
215: VecDuplicate(tao->solution, &qp->DG);
216: VecDuplicate(tao->solution, &qp->S);
217: VecDuplicate(tao->solution, &qp->Z);
218: VecDuplicate(tao->solution, &qp->DZ);
219: VecDuplicate(tao->solution, &qp->GZwork);
220: VecDuplicate(tao->solution, &qp->R3);
222: VecDuplicate(tao->solution, &qp->T);
223: VecDuplicate(tao->solution, &qp->DT);
224: VecDuplicate(tao->solution, &qp->DS);
225: VecDuplicate(tao->solution, &qp->TSwork);
226: VecDuplicate(tao->solution, &qp->R5);
227: qp->m = 2 * qp->n;
228: return 0;
229: }
231: static PetscErrorCode TaoSolve_BQPIP(Tao tao)
232: {
233: TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
234: PetscInt its;
235: PetscReal d1, d2, ksptol, sigmamu;
236: PetscReal gnorm, dstep, pstep, step = 0;
237: PetscReal gap[4];
238: PetscBool getdiagop;
240: qp->dobj = 0.0;
241: qp->pobj = 1.0;
242: qp->gap = 10.0;
243: qp->rgap = 1.0;
244: qp->mu = 1.0;
245: qp->dinfeas = 1.0;
246: qp->psteplength = 0.0;
247: qp->dsteplength = 0.0;
249: /* TODO
250: - Remove fixed variables and treat them correctly
251: - Use index sets for the infinite versus finite bounds
252: - Update remaining code for fixed and free variables
253: - Fix inexact solves for predictor and corrector
254: */
256: /* Tighten infinite bounds, things break when we don't do this
257: -- see test_bqpip.c
258: */
259: TaoComputeVariableBounds(tao);
260: VecSet(qp->XU, 1.0e20);
261: VecSet(qp->XL, -1.0e20);
262: if (tao->XL) VecPointwiseMax(qp->XL, qp->XL, tao->XL);
263: if (tao->XU) VecPointwiseMin(qp->XU, qp->XU, tao->XU);
264: VecMedian(qp->XL, tao->solution, qp->XU, tao->solution);
266: /* Evaluate gradient and Hessian at zero to get the correct values
267: without contaminating them with numerical artifacts.
268: */
269: VecSet(qp->Work, 0);
270: TaoComputeObjectiveAndGradient(tao, qp->Work, &qp->d, qp->C);
271: TaoComputeHessian(tao, qp->Work, tao->hessian, tao->hessian_pre);
272: MatHasOperation(tao->hessian, MATOP_GET_DIAGONAL, &getdiagop);
273: if (getdiagop) MatGetDiagonal(tao->hessian, qp->HDiag);
275: /* Initialize starting point and residuals */
276: QPIPSetInitialPoint(qp, tao);
277: QPIPComputeResidual(qp, tao);
279: /* Enter main loop */
280: tao->reason = TAO_CONTINUE_ITERATING;
281: while (1) {
282: /* Check Stopping Condition */
283: gnorm = PetscSqrtScalar(qp->gap + qp->dinfeas);
284: TaoLogConvergenceHistory(tao, qp->pobj, gnorm, qp->pinfeas, tao->ksp_its);
285: TaoMonitor(tao, tao->niter, qp->pobj, gnorm, qp->pinfeas, step);
286: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
287: if (tao->reason != TAO_CONTINUE_ITERATING) break;
288: /* Call general purpose update function */
289: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
290: tao->niter++;
291: tao->ksp_its = 0;
293: /*
294: Dual Infeasibility Direction should already be in the right
295: hand side from computing the residuals
296: */
298: QPIPComputeNormFromCentralPath(qp, &d1);
300: VecSet(qp->DZ, 0.0);
301: VecSet(qp->DS, 0.0);
303: /*
304: Compute the Primal Infeasiblitiy RHS and the
305: Diagonal Matrix to be added to H and store in Work
306: */
307: VecPointwiseDivide(qp->DiagAxpy, qp->Z, qp->G);
308: VecPointwiseMult(qp->GZwork, qp->DiagAxpy, qp->R3);
309: VecAXPY(qp->RHS, -1.0, qp->GZwork);
311: VecPointwiseDivide(qp->TSwork, qp->S, qp->T);
312: VecAXPY(qp->DiagAxpy, 1.0, qp->TSwork);
313: VecPointwiseMult(qp->TSwork, qp->TSwork, qp->R5);
314: VecAXPY(qp->RHS, -1.0, qp->TSwork);
316: /* Determine the solving tolerance */
317: ksptol = qp->mu / 10.0;
318: ksptol = PetscMin(ksptol, 0.001);
319: KSPSetTolerances(tao->ksp, ksptol, 1e-30, 1e30, PetscMax(10, qp->n));
321: /* Shift the diagonals of the Hessian matrix */
322: MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES);
323: if (!getdiagop) {
324: VecCopy(qp->DiagAxpy, qp->HDiag);
325: VecScale(qp->HDiag, -1.0);
326: }
327: MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY);
328: MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY);
330: KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);
331: KSPSolve(tao->ksp, qp->RHS, tao->stepdirection);
332: KSPGetIterationNumber(tao->ksp, &its);
333: tao->ksp_its += its;
334: tao->ksp_tot_its += its;
336: /* Restore the true diagonal of the Hessian matrix */
337: if (getdiagop) {
338: MatDiagonalSet(tao->hessian, qp->HDiag, INSERT_VALUES);
339: } else {
340: MatDiagonalSet(tao->hessian, qp->HDiag, ADD_VALUES);
341: }
342: MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY);
343: MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY);
344: QPIPComputeStepDirection(qp, tao);
345: QPIPStepLength(qp);
347: /* Calculate New Residual R1 in Work vector */
348: MatMult(tao->hessian, tao->stepdirection, qp->RHS2);
349: VecAXPY(qp->RHS2, 1.0, qp->DS);
350: VecAXPY(qp->RHS2, -1.0, qp->DZ);
351: VecAYPX(qp->RHS2, qp->dsteplength, tao->gradient);
353: VecNorm(qp->RHS2, NORM_2, &qp->dinfeas);
354: VecDot(qp->DZ, qp->DG, gap);
355: VecDot(qp->DS, qp->DT, gap + 1);
357: qp->rnorm = (qp->dinfeas + qp->psteplength * qp->pinfeas) / (qp->m + qp->n);
358: pstep = qp->psteplength;
359: step = PetscMin(qp->psteplength, qp->dsteplength);
360: sigmamu = (pstep * pstep * (gap[0] + gap[1]) + (1 - pstep) * qp->gap) / qp->m;
362: if (qp->predcorr && step < 0.9) {
363: if (sigmamu < qp->mu) {
364: sigmamu = sigmamu / qp->mu;
365: sigmamu = sigmamu * sigmamu * sigmamu;
366: } else {
367: sigmamu = 1.0;
368: }
369: sigmamu = sigmamu * qp->mu;
371: /* Compute Corrector Step */
372: VecPointwiseMult(qp->DZ, qp->DG, qp->DZ);
373: VecScale(qp->DZ, -1.0);
374: VecShift(qp->DZ, sigmamu);
375: VecPointwiseDivide(qp->DZ, qp->DZ, qp->G);
377: VecPointwiseMult(qp->DS, qp->DS, qp->DT);
378: VecScale(qp->DS, -1.0);
379: VecShift(qp->DS, sigmamu);
380: VecPointwiseDivide(qp->DS, qp->DS, qp->T);
382: VecCopy(qp->DZ, qp->RHS2);
383: VecAXPY(qp->RHS2, -1.0, qp->DS);
384: VecAXPY(qp->RHS2, 1.0, qp->RHS);
386: /* Approximately solve the linear system */
387: MatDiagonalSet(tao->hessian, qp->DiagAxpy, ADD_VALUES);
388: if (!getdiagop) {
389: VecCopy(qp->DiagAxpy, qp->HDiag);
390: VecScale(qp->HDiag, -1.0);
391: }
392: MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY);
393: MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY);
395: /* Solve using the previous tolerances that were set */
396: KSPSolve(tao->ksp, qp->RHS2, tao->stepdirection);
397: KSPGetIterationNumber(tao->ksp, &its);
398: tao->ksp_its += its;
399: tao->ksp_tot_its += its;
401: if (getdiagop) {
402: MatDiagonalSet(tao->hessian, qp->HDiag, INSERT_VALUES);
403: } else {
404: MatDiagonalSet(tao->hessian, qp->HDiag, ADD_VALUES);
405: }
406: MatAssemblyBegin(tao->hessian, MAT_FINAL_ASSEMBLY);
407: MatAssemblyEnd(tao->hessian, MAT_FINAL_ASSEMBLY);
408: QPIPComputeStepDirection(qp, tao);
409: QPIPStepLength(qp);
410: } /* End Corrector step */
412: /* Take the step */
413: dstep = qp->dsteplength;
415: VecAXPY(qp->Z, dstep, qp->DZ);
416: VecAXPY(qp->S, dstep, qp->DS);
417: VecAXPY(tao->solution, dstep, tao->stepdirection);
418: VecAXPY(qp->G, dstep, qp->DG);
419: VecAXPY(qp->T, dstep, qp->DT);
421: /* Compute Residuals */
422: QPIPComputeResidual(qp, tao);
424: /* Evaluate quadratic function */
425: MatMult(tao->hessian, tao->solution, qp->Work);
427: VecDot(tao->solution, qp->Work, &d1);
428: VecDot(tao->solution, qp->C, &d2);
429: VecDot(qp->G, qp->Z, gap);
430: VecDot(qp->T, qp->S, gap + 1);
432: /* Compute the duality gap */
433: qp->pobj = d1 / 2.0 + d2 + qp->d;
434: qp->gap = gap[0] + gap[1];
435: qp->dobj = qp->pobj - qp->gap;
436: if (qp->m > 0) qp->mu = qp->gap / (qp->m);
437: qp->rgap = qp->gap / (PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0);
438: } /* END MAIN LOOP */
439: return 0;
440: }
442: static PetscErrorCode TaoView_BQPIP(Tao tao, PetscViewer viewer)
443: {
444: return 0;
445: }
447: static PetscErrorCode TaoSetFromOptions_BQPIP(Tao tao, PetscOptionItems *PetscOptionsObject)
448: {
449: TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
451: PetscOptionsHeadBegin(PetscOptionsObject, "Interior point method for bound constrained quadratic optimization");
452: PetscOptionsInt("-tao_bqpip_predcorr", "Use a predictor-corrector method", "", qp->predcorr, &qp->predcorr, NULL);
453: PetscOptionsHeadEnd();
454: KSPSetFromOptions(tao->ksp);
455: return 0;
456: }
458: static PetscErrorCode TaoDestroy_BQPIP(Tao tao)
459: {
460: TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
462: if (tao->setupcalled) {
463: VecDestroy(&qp->G);
464: VecDestroy(&qp->DG);
465: VecDestroy(&qp->Z);
466: VecDestroy(&qp->DZ);
467: VecDestroy(&qp->GZwork);
468: VecDestroy(&qp->R3);
469: VecDestroy(&qp->S);
470: VecDestroy(&qp->DS);
471: VecDestroy(&qp->T);
473: VecDestroy(&qp->DT);
474: VecDestroy(&qp->TSwork);
475: VecDestroy(&qp->R5);
476: VecDestroy(&qp->HDiag);
477: VecDestroy(&qp->Work);
478: VecDestroy(&qp->XL);
479: VecDestroy(&qp->XU);
480: VecDestroy(&qp->DiagAxpy);
481: VecDestroy(&qp->RHS);
482: VecDestroy(&qp->RHS2);
483: VecDestroy(&qp->C);
484: }
485: KSPDestroy(&tao->ksp);
486: PetscFree(tao->data);
487: return 0;
488: }
490: static PetscErrorCode TaoComputeDual_BQPIP(Tao tao, Vec DXL, Vec DXU)
491: {
492: TAO_BQPIP *qp = (TAO_BQPIP *)tao->data;
494: VecCopy(qp->Z, DXL);
495: VecCopy(qp->S, DXU);
496: VecScale(DXU, -1.0);
497: return 0;
498: }
500: /*MC
501: TAOBQPIP - interior-point method for quadratic programs with
502: box constraints.
504: Notes:
505: This algorithms solves quadratic problems only, the Hessian will
506: only be computed once.
508: Options Database Keys:
509: . -tao_bqpip_predcorr - use a predictor/corrector method
511: Level: beginner
512: M*/
514: PETSC_EXTERN PetscErrorCode TaoCreate_BQPIP(Tao tao)
515: {
516: TAO_BQPIP *qp;
518: PetscNew(&qp);
520: tao->ops->setup = TaoSetup_BQPIP;
521: tao->ops->solve = TaoSolve_BQPIP;
522: tao->ops->view = TaoView_BQPIP;
523: tao->ops->setfromoptions = TaoSetFromOptions_BQPIP;
524: tao->ops->destroy = TaoDestroy_BQPIP;
525: tao->ops->computedual = TaoComputeDual_BQPIP;
527: /* Override default settings (unless already changed) */
528: if (!tao->max_it_changed) tao->max_it = 100;
529: if (!tao->max_funcs_changed) tao->max_funcs = 500;
530: #if defined(PETSC_USE_REAL_SINGLE)
531: if (!tao->catol_changed) tao->catol = 1e-6;
532: #else
533: if (!tao->catol_changed) tao->catol = 1e-12;
534: #endif
536: /* Initialize pointers and variables */
537: qp->n = 0;
538: qp->m = 0;
540: qp->predcorr = 1;
541: tao->data = (void *)qp;
543: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
544: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
545: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
546: KSPSetType(tao->ksp, KSPCG);
547: KSPSetTolerances(tao->ksp, 1e-14, 1e-30, 1e30, PetscMax(10, qp->n));
548: return 0;
549: }