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;qp->dinfeas=1.0;qp->pinfeas=1.0;
86: while ((qp->dinfeas+qp->pinfeas)/(qp->m+qp->n) >= qp->mu) {
87: VecScale(qp->G,two);
88: VecScale(qp->Z,two);
89: VecScale(qp->S,two);
90: VecScale(qp->T,two);
92: QPIPComputeResidual(qp,tao);
94: VecCopy(tao->solution,qp->R3);
95: VecAXPY(qp->R3,-1.0,qp->G);
96: VecAXPY(qp->R3,-1.0,qp->XL);
98: VecCopy(tao->solution,qp->R5);
99: VecAXPY(qp->R5,1.0,qp->T);
100: VecAXPY(qp->R5,-1.0,qp->XU);
102: VecNorm(qp->R3,NORM_INFINITY,&gap1);
103: VecNorm(qp->R5,NORM_INFINITY,&gap2);
104: qp->pinfeas=PetscMax(gap1,gap2);
106: /* Compute the duality gap */
107: VecDot(qp->G,qp->Z,&gap1);
108: VecDot(qp->T,qp->S,&gap2);
110: qp->gap = gap1+gap2;
111: qp->dobj = qp->pobj - qp->gap;
112: if (qp->m>0) {
113: qp->mu=qp->gap/(qp->m);
114: } else {
115: qp->mu=0.0;
116: }
117: qp->rgap=qp->gap/(PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0);
118: }
119: return 0;
120: }
122: static PetscErrorCode QPIPStepLength(TAO_BQPIP *qp)
123: {
124: PetscReal tstep1,tstep2,tstep3,tstep4,tstep;
126: /* Compute stepsize to the boundary */
127: VecStepMax(qp->G,qp->DG,&tstep1);
128: VecStepMax(qp->T,qp->DT,&tstep2);
129: VecStepMax(qp->S,qp->DS,&tstep3);
130: VecStepMax(qp->Z,qp->DZ,&tstep4);
132: tstep = PetscMin(tstep1,tstep2);
133: qp->psteplength = PetscMin(0.95*tstep,1.0);
135: tstep = PetscMin(tstep3,tstep4);
136: qp->dsteplength = PetscMin(0.95*tstep,1.0);
138: qp->psteplength = PetscMin(qp->psteplength,qp->dsteplength);
139: qp->dsteplength = qp->psteplength;
140: return 0;
141: }
143: static PetscErrorCode QPIPComputeNormFromCentralPath(TAO_BQPIP *qp,PetscReal *norm)
144: {
145: PetscReal gap[2],mu[2],nmu;
147: VecPointwiseMult(qp->GZwork,qp->G,qp->Z);
148: VecPointwiseMult(qp->TSwork,qp->T,qp->S);
149: VecNorm(qp->TSwork,NORM_1,&mu[0]);
150: VecNorm(qp->GZwork,NORM_1,&mu[1]);
152: nmu=-(mu[0]+mu[1])/qp->m;
154: VecShift(qp->GZwork,nmu);
155: VecShift(qp->TSwork,nmu);
157: VecNorm(qp->GZwork,NORM_2,&gap[0]);
158: VecNorm(qp->TSwork,NORM_2,&gap[1]);
159: gap[0]*=gap[0];
160: gap[1]*=gap[1];
162: qp->pathnorm=PetscSqrtScalar(gap[0]+gap[1]);
163: *norm=qp->pathnorm;
164: return 0;
165: }
167: static PetscErrorCode QPIPComputeStepDirection(TAO_BQPIP *qp,Tao tao)
168: {
169: /* Calculate DG */
170: VecCopy(tao->stepdirection,qp->DG);
171: VecAXPY(qp->DG,1.0,qp->R3);
173: /* Calculate DT */
174: VecCopy(tao->stepdirection,qp->DT);
175: VecScale(qp->DT,-1.0);
176: VecAXPY(qp->DT,-1.0,qp->R5);
178: /* Calculate DZ */
179: VecAXPY(qp->DZ,-1.0,qp->Z);
180: VecPointwiseDivide(qp->GZwork,qp->DG,qp->G);
181: VecPointwiseMult(qp->GZwork,qp->GZwork,qp->Z);
182: VecAXPY(qp->DZ,-1.0,qp->GZwork);
184: /* Calculate DS */
185: VecAXPY(qp->DS,-1.0,qp->S);
186: VecPointwiseDivide(qp->TSwork,qp->DT,qp->T);
187: VecPointwiseMult(qp->TSwork,qp->TSwork,qp->S);
188: VecAXPY(qp->DS,-1.0,qp->TSwork);
189: return 0;
190: }
192: static PetscErrorCode TaoSetup_BQPIP(Tao tao)
193: {
194: TAO_BQPIP *qp =(TAO_BQPIP*)tao->data;
196: /* Set pointers to Data */
197: VecGetSize(tao->solution,&qp->n);
199: /* Allocate some arrays */
200: if (!tao->gradient) {
201: VecDuplicate(tao->solution,&tao->gradient);
202: }
203: if (!tao->stepdirection) {
204: VecDuplicate(tao->solution,&tao->stepdirection);
205: }
206: if (!tao->XL) {
207: VecDuplicate(tao->solution,&tao->XL);
208: VecSet(tao->XL,PETSC_NINFINITY);
209: }
210: if (!tao->XU) {
211: VecDuplicate(tao->solution,&tao->XU);
212: VecSet(tao->XU,PETSC_INFINITY);
213: }
215: VecDuplicate(tao->solution,&qp->Work);
216: VecDuplicate(tao->solution,&qp->XU);
217: VecDuplicate(tao->solution,&qp->XL);
218: VecDuplicate(tao->solution,&qp->HDiag);
219: VecDuplicate(tao->solution,&qp->DiagAxpy);
220: VecDuplicate(tao->solution,&qp->RHS);
221: VecDuplicate(tao->solution,&qp->RHS2);
222: VecDuplicate(tao->solution,&qp->C);
224: VecDuplicate(tao->solution,&qp->G);
225: VecDuplicate(tao->solution,&qp->DG);
226: VecDuplicate(tao->solution,&qp->S);
227: VecDuplicate(tao->solution,&qp->Z);
228: VecDuplicate(tao->solution,&qp->DZ);
229: VecDuplicate(tao->solution,&qp->GZwork);
230: VecDuplicate(tao->solution,&qp->R3);
232: VecDuplicate(tao->solution,&qp->T);
233: VecDuplicate(tao->solution,&qp->DT);
234: VecDuplicate(tao->solution,&qp->DS);
235: VecDuplicate(tao->solution,&qp->TSwork);
236: VecDuplicate(tao->solution,&qp->R5);
237: qp->m=2*qp->n;
238: return 0;
239: }
241: static PetscErrorCode TaoSolve_BQPIP(Tao tao)
242: {
243: TAO_BQPIP *qp = (TAO_BQPIP*)tao->data;
244: PetscInt its;
245: PetscReal d1,d2,ksptol,sigmamu;
246: PetscReal gnorm,dstep,pstep,step=0;
247: PetscReal gap[4];
248: PetscBool getdiagop;
250: qp->dobj = 0.0;
251: qp->pobj = 1.0;
252: qp->gap = 10.0;
253: qp->rgap = 1.0;
254: qp->mu = 1.0;
255: qp->dinfeas = 1.0;
256: qp->psteplength = 0.0;
257: qp->dsteplength = 0.0;
259: /* TODO
260: - Remove fixed variables and treat them correctly
261: - Use index sets for the infinite versus finite bounds
262: - Update remaining code for fixed and free variables
263: - Fix inexact solves for predictor and corrector
264: */
266: /* Tighten infinite bounds, things break when we don't do this
267: -- see test_bqpip.c
268: */
269: TaoComputeVariableBounds(tao);
270: VecSet(qp->XU,1.0e20);
271: VecSet(qp->XL,-1.0e20);
272: VecPointwiseMax(qp->XL,qp->XL,tao->XL);
273: VecPointwiseMin(qp->XU,qp->XU,tao->XU);
274: VecMedian(qp->XL,tao->solution,qp->XU,tao->solution);
276: /* Evaluate gradient and Hessian at zero to get the correct values
277: without contaminating them with numerical artifacts.
278: */
279: VecSet(qp->Work,0);
280: TaoComputeObjectiveAndGradient(tao,qp->Work,&qp->d,qp->C);
281: TaoComputeHessian(tao,qp->Work,tao->hessian,tao->hessian_pre);
282: MatHasOperation(tao->hessian,MATOP_GET_DIAGONAL,&getdiagop);
283: if (getdiagop) {
284: MatGetDiagonal(tao->hessian,qp->HDiag);
285: }
287: /* Initialize starting point and residuals */
288: QPIPSetInitialPoint(qp,tao);
289: QPIPComputeResidual(qp,tao);
291: /* Enter main loop */
292: tao->reason = TAO_CONTINUE_ITERATING;
293: while (1) {
295: /* Check Stopping Condition */
296: gnorm = PetscSqrtScalar(qp->gap + qp->dinfeas);
297: TaoLogConvergenceHistory(tao,qp->pobj,gnorm,qp->pinfeas,tao->ksp_its);
298: TaoMonitor(tao,tao->niter,qp->pobj,gnorm,qp->pinfeas,step);
299: (*tao->ops->convergencetest)(tao,tao->cnvP);
300: if (tao->reason != TAO_CONTINUE_ITERATING) break;
301: /* Call general purpose update function */
302: if (tao->ops->update) {
303: (*tao->ops->update)(tao, tao->niter, tao->user_update);
304: }
305: tao->niter++;
306: tao->ksp_its = 0;
308: /*
309: Dual Infeasibility Direction should already be in the right
310: hand side from computing the residuals
311: */
313: QPIPComputeNormFromCentralPath(qp,&d1);
315: VecSet(qp->DZ,0.0);
316: VecSet(qp->DS,0.0);
318: /*
319: Compute the Primal Infeasiblitiy RHS and the
320: Diagonal Matrix to be added to H and store in Work
321: */
322: VecPointwiseDivide(qp->DiagAxpy,qp->Z,qp->G);
323: VecPointwiseMult(qp->GZwork,qp->DiagAxpy,qp->R3);
324: VecAXPY(qp->RHS,-1.0,qp->GZwork);
326: VecPointwiseDivide(qp->TSwork,qp->S,qp->T);
327: VecAXPY(qp->DiagAxpy,1.0,qp->TSwork);
328: VecPointwiseMult(qp->TSwork,qp->TSwork,qp->R5);
329: VecAXPY(qp->RHS,-1.0,qp->TSwork);
331: /* Determine the solving tolerance */
332: ksptol = qp->mu/10.0;
333: ksptol = PetscMin(ksptol,0.001);
334: KSPSetTolerances(tao->ksp,ksptol,1e-30,1e30,PetscMax(10,qp->n));
336: /* Shift the diagonals of the Hessian matrix */
337: MatDiagonalSet(tao->hessian,qp->DiagAxpy,ADD_VALUES);
338: if (!getdiagop) {
339: VecCopy(qp->DiagAxpy,qp->HDiag);
340: VecScale(qp->HDiag,-1.0);
341: }
342: MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
343: MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);
345: KSPSetOperators(tao->ksp,tao->hessian,tao->hessian_pre);
346: KSPSolve(tao->ksp,qp->RHS,tao->stepdirection);
347: KSPGetIterationNumber(tao->ksp,&its);
348: tao->ksp_its += its;
349: tao->ksp_tot_its += its;
351: /* Restore the true diagonal of the Hessian matrix */
352: if (getdiagop) {
353: MatDiagonalSet(tao->hessian,qp->HDiag,INSERT_VALUES);
354: } else {
355: MatDiagonalSet(tao->hessian,qp->HDiag,ADD_VALUES);
356: }
357: MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
358: MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);
359: QPIPComputeStepDirection(qp,tao);
360: QPIPStepLength(qp);
362: /* Calculate New Residual R1 in Work vector */
363: MatMult(tao->hessian,tao->stepdirection,qp->RHS2);
364: VecAXPY(qp->RHS2,1.0,qp->DS);
365: VecAXPY(qp->RHS2,-1.0,qp->DZ);
366: VecAYPX(qp->RHS2,qp->dsteplength,tao->gradient);
368: VecNorm(qp->RHS2,NORM_2,&qp->dinfeas);
369: VecDot(qp->DZ,qp->DG,gap);
370: VecDot(qp->DS,qp->DT,gap+1);
372: qp->rnorm = (qp->dinfeas+qp->psteplength*qp->pinfeas)/(qp->m+qp->n);
373: pstep = qp->psteplength;
374: step = PetscMin(qp->psteplength,qp->dsteplength);
375: sigmamu = (pstep*pstep*(gap[0]+gap[1]) + (1 - pstep)*qp->gap)/qp->m;
377: if (qp->predcorr && step < 0.9) {
378: if (sigmamu < qp->mu) {
379: sigmamu = sigmamu/qp->mu;
380: sigmamu = sigmamu*sigmamu*sigmamu;
381: } else {
382: sigmamu = 1.0;
383: }
384: sigmamu = sigmamu*qp->mu;
386: /* Compute Corrector Step */
387: VecPointwiseMult(qp->DZ,qp->DG,qp->DZ);
388: VecScale(qp->DZ,-1.0);
389: VecShift(qp->DZ,sigmamu);
390: VecPointwiseDivide(qp->DZ,qp->DZ,qp->G);
392: VecPointwiseMult(qp->DS,qp->DS,qp->DT);
393: VecScale(qp->DS,-1.0);
394: VecShift(qp->DS,sigmamu);
395: VecPointwiseDivide(qp->DS,qp->DS,qp->T);
397: VecCopy(qp->DZ,qp->RHS2);
398: VecAXPY(qp->RHS2,-1.0,qp->DS);
399: VecAXPY(qp->RHS2,1.0,qp->RHS);
401: /* Approximately solve the linear system */
402: MatDiagonalSet(tao->hessian,qp->DiagAxpy,ADD_VALUES);
403: if (!getdiagop) {
404: VecCopy(qp->DiagAxpy,qp->HDiag);
405: VecScale(qp->HDiag,-1.0);
406: }
407: MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
408: MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);
410: /* Solve using the previous tolerances that were set */
411: KSPSolve(tao->ksp,qp->RHS2,tao->stepdirection);
412: KSPGetIterationNumber(tao->ksp,&its);
413: tao->ksp_its += its;
414: tao->ksp_tot_its += its;
416: if (getdiagop) {
417: MatDiagonalSet(tao->hessian,qp->HDiag,INSERT_VALUES);
418: } else {
419: MatDiagonalSet(tao->hessian,qp->HDiag,ADD_VALUES);
420: }
421: MatAssemblyBegin(tao->hessian,MAT_FINAL_ASSEMBLY);
422: MatAssemblyEnd(tao->hessian,MAT_FINAL_ASSEMBLY);
423: QPIPComputeStepDirection(qp,tao);
424: QPIPStepLength(qp);
425: } /* End Corrector step */
427: /* Take the step */
428: dstep = qp->dsteplength;
430: VecAXPY(qp->Z,dstep,qp->DZ);
431: VecAXPY(qp->S,dstep,qp->DS);
432: VecAXPY(tao->solution,dstep,tao->stepdirection);
433: VecAXPY(qp->G,dstep,qp->DG);
434: VecAXPY(qp->T,dstep,qp->DT);
436: /* Compute Residuals */
437: QPIPComputeResidual(qp,tao);
439: /* Evaluate quadratic function */
440: MatMult(tao->hessian,tao->solution,qp->Work);
442: VecDot(tao->solution,qp->Work,&d1);
443: VecDot(tao->solution,qp->C,&d2);
444: VecDot(qp->G,qp->Z,gap);
445: VecDot(qp->T,qp->S,gap+1);
447: /* Compute the duality gap */
448: qp->pobj = d1/2.0 + d2+qp->d;
449: qp->gap = gap[0]+gap[1];
450: qp->dobj = qp->pobj - qp->gap;
451: if (qp->m > 0) {
452: qp->mu = qp->gap/(qp->m);
453: }
454: qp->rgap = qp->gap/(PetscAbsReal(qp->dobj) + PetscAbsReal(qp->pobj) + 1.0);
455: } /* END MAIN LOOP */
456: return 0;
457: }
459: static PetscErrorCode TaoView_BQPIP(Tao tao,PetscViewer viewer)
460: {
461: return 0;
462: }
464: static PetscErrorCode TaoSetFromOptions_BQPIP(PetscOptionItems *PetscOptionsObject,Tao tao)
465: {
466: TAO_BQPIP *qp = (TAO_BQPIP*)tao->data;
468: PetscOptionsHead(PetscOptionsObject,"Interior point method for bound constrained quadratic optimization");
469: PetscOptionsInt("-tao_bqpip_predcorr","Use a predictor-corrector method","",qp->predcorr,&qp->predcorr,NULL);
470: PetscOptionsTail();
471: KSPSetFromOptions(tao->ksp);
472: return 0;
473: }
475: static PetscErrorCode TaoDestroy_BQPIP(Tao tao)
476: {
477: TAO_BQPIP *qp = (TAO_BQPIP*)tao->data;
479: if (tao->setupcalled) {
480: VecDestroy(&qp->G);
481: VecDestroy(&qp->DG);
482: VecDestroy(&qp->Z);
483: VecDestroy(&qp->DZ);
484: VecDestroy(&qp->GZwork);
485: VecDestroy(&qp->R3);
486: VecDestroy(&qp->S);
487: VecDestroy(&qp->DS);
488: VecDestroy(&qp->T);
490: VecDestroy(&qp->DT);
491: VecDestroy(&qp->TSwork);
492: VecDestroy(&qp->R5);
493: VecDestroy(&qp->HDiag);
494: VecDestroy(&qp->Work);
495: VecDestroy(&qp->XL);
496: VecDestroy(&qp->XU);
497: VecDestroy(&qp->DiagAxpy);
498: VecDestroy(&qp->RHS);
499: VecDestroy(&qp->RHS2);
500: VecDestroy(&qp->C);
501: }
502: PetscFree(tao->data);
503: return 0;
504: }
506: static PetscErrorCode TaoComputeDual_BQPIP(Tao tao,Vec DXL,Vec DXU)
507: {
508: TAO_BQPIP *qp = (TAO_BQPIP*)tao->data;
510: VecCopy(qp->Z,DXL);
511: VecCopy(qp->S,DXU);
512: VecScale(DXU,-1.0);
513: return 0;
514: }
516: /*MC
517: TAOBQPIP - interior-point method for quadratic programs with
518: box constraints.
520: Notes:
521: This algorithms solves quadratic problems only, the Hessian will
522: only be computed once.
524: Options Database Keys:
525: . -tao_bqpip_predcorr - use a predictor/corrector method
527: Level: beginner
528: M*/
530: PETSC_EXTERN PetscErrorCode TaoCreate_BQPIP(Tao tao)
531: {
532: TAO_BQPIP *qp;
534: PetscNewLog(tao,&qp);
536: tao->ops->setup = TaoSetup_BQPIP;
537: tao->ops->solve = TaoSolve_BQPIP;
538: tao->ops->view = TaoView_BQPIP;
539: tao->ops->setfromoptions = TaoSetFromOptions_BQPIP;
540: tao->ops->destroy = TaoDestroy_BQPIP;
541: tao->ops->computedual = TaoComputeDual_BQPIP;
543: /* Override default settings (unless already changed) */
544: if (!tao->max_it_changed) tao->max_it=100;
545: if (!tao->max_funcs_changed) tao->max_funcs = 500;
546: #if defined(PETSC_USE_REAL_SINGLE)
547: if (!tao->catol_changed) tao->catol=1e-6;
548: #else
549: if (!tao->catol_changed) tao->catol=1e-12;
550: #endif
552: /* Initialize pointers and variables */
553: qp->n = 0;
554: qp->m = 0;
556: qp->predcorr = 1;
557: tao->data = (void*)qp;
559: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
560: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
561: KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);
562: KSPSetType(tao->ksp,KSPCG);
563: KSPSetTolerances(tao->ksp,1e-14,1e-30,1e30,PetscMax(10,qp->n));
564: return 0;
565: }