Actual source code: tron.c

  1: /*$Id$*/

  3: #include "tron.h"       /*I "tao_solver.h" I*/

  5: /* TRON Routines */
  6: static int TaoGradProjections(TAO_SOLVER,TAO_TRON *);
  7: static int TronCheckOptimalFace(TaoVec*, TaoVec*, TaoVec*, TaoVec*, TaoVec*, 
  8:                                 TaoIndexSet*, TaoIndexSet*, TaoTruth *optimal);

 10: /*------------------------------------------------------------*/
 13: static int TaoSetDown_TRON(TAO_SOLVER tao, void*solver)
 14: {
 15:   TAO_TRON *tron = (TAO_TRON *)solver;
 16:   int      info;
 17:   /* Free allocated memory in TRON structure */
 18:   TaoFunctionBegin;
 19:   
 20:   info = TaoVecDestroy(tron->X_New);CHKERRQ(info);
 21:   info = TaoVecDestroy(tron->G_New);CHKERRQ(info);
 22:   info = TaoVecDestroy(tron->DX);CHKERRQ(info);tron->DX=0;
 23:   info = TaoVecDestroy(tron->Work);CHKERRQ(info);
 24:   info = TaoVecDestroy(tron->DXFree);CHKERRQ(info);
 25:   info = TaoVecDestroy(tron->R);CHKERRQ(info);
 26:   info = TaoVecDestroy(tron->G);CHKERRQ(info);
 27:   info = TaoVecDestroy(tron->PG);CHKERRQ(info);
 28:   info = TaoVecDestroy(tron->XL);CHKERRQ(info);
 29:   info = TaoVecDestroy(tron->XU);CHKERRQ(info);
 30:   
 31:   info = TaoIndexSetDestroy(tron->Free_Local);CHKERRQ(info);
 32:   info = TaoIndexSetDestroy(tron->TT);CHKERRQ(info);
 33:   info = TaoMatDestroy(tron->Hsub);CHKERRQ(info);

 35:   info = TaoSetLagrangianGradientVector(tao,0);CHKERRQ(info);
 36:   info = TaoSetStepDirectionVector(tao,0);CHKERRQ(info);
 37:   info = TaoSetVariableBounds(tao,0,0);CHKERRQ(info);

 39:   TaoFunctionReturn(0);
 40: }

 42: /*------------------------------------------------------------*/
 45: static int TaoSetOptions_TRON(TAO_SOLVER tao, void*solver)
 46: {
 47:   TAO_TRON  *tron = (TAO_TRON *)solver;
 48:   int        info;
 49:   TaoInt     ival;
 50:   TaoTruth flg;

 52:   TaoFunctionBegin;

 54:   info = TaoOptionsHead("Newton Trust Region Method for bound constrained optimization");CHKERRQ(info);

 56:   info = TaoOptionInt("-tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg);
 57:   CHKERRQ(info);

 59:   info = TaoOptionInt("-redistribute","Redistribute Free variables (> 1 processors, only)","TaoPetscISType",1,&ival,&flg); CHKERRQ(info);

 61:   info = TaoOptionName("-submatrixfree","Mask full matrix instead of extract submatrices","TaoPetscISType",&flg); CHKERRQ(info);

 63:   info = TaoOptionsTail();CHKERRQ(info);
 64:   info = TaoLineSearchSetFromOptions(tao);CHKERRQ(info);

 66:   TaoFunctionReturn(0);
 67: }

 69: /*------------------------------------------------------------*/
 72: static int TaoView_TRON(TAO_SOLVER tao,void*solver)
 73: {
 74:   TAO_TRON  *tron = (TAO_TRON *)solver;
 75:   int        info;

 77:   TaoFunctionBegin;
 78:   /*
 79:   info = TaoPrintf1(tao," Variables, Total: %d,",tron->n);
 80:   info = TaoPrintf3(tao,"Free: %d,  Binding: %d \n",
 81:                     tron->n_free, tron->n - tron->n_free,
 82:                     tron->n_bind);CHKERRQ(info);
 83:   info = TaoPrintf1(tao,"            Equal lower bound: %d,",
 84:                     tron->n_lower);CHKERRQ(info);
 85:   info = TaoPrintf1(tao,"  Equal upper bound: %d \n",
 86:                     tron->n_upper);CHKERRQ(info);
 87:   */
 88:   info = TaoPrintInt(tao," Total PG its: %d,",tron->total_gp_its);CHKERRQ(info);
 89:   info = TaoPrintDouble(tao," PG tolerance: %4.3f \n",tron->pg_ftol);CHKERRQ(info);
 90:   info = TaoLineSearchView(tao);CHKERRQ(info);
 91:   info = TaoPrintStatement(tao,"  Linear Solver minimizes quadratic over Trust Region: \n");CHKERRQ(info);

 93:   TaoFunctionReturn(0);
 94: }


 97: /* ---------------------------------------------------------- */
100: static int TaoSetUp_TRON(TAO_SOLVER tao, void*solver){

102:   int info;
103:   TAO_TRON *tron = (TAO_TRON *)solver;
104:   TaoVec* X;
105:   TaoMat *HH;
106:   TaoIndexSet *TIS;

108:   TaoFunctionBegin;
109:   info = TaoGetSolution(tao,&tron->X);CHKERRQ(info); X=tron->X;
110:   info = TaoGetHessian(tao,&tron->H);CHKERRQ(info);  HH=tron->H;

112:   /* Allocate some arrays */
113:   info = X->Clone(&tron->DX); CHKERRQ(info);
114:   info = X->Clone(&tron->X_New); CHKERRQ(info);
115:   info = X->Clone(&tron->G_New); CHKERRQ(info);
116:   info = X->Clone(&tron->Work); CHKERRQ(info);
117:   info = X->Clone(&tron->DXFree); CHKERRQ(info);
118:   info = X->Clone(&tron->R); CHKERRQ(info);
119:   info = X->Clone(&tron->G); CHKERRQ(info);
120:   info = X->Clone(&tron->PG); CHKERRQ(info);
121:   info = X->Clone(&tron->XL); CHKERRQ(info);
122:   info = X->Clone(&tron->XU); CHKERRQ(info);

124:   info = TaoSetLagrangianGradientVector(tao,tron->PG);CHKERRQ(info);
125:   info = TaoSetStepDirectionVector(tao,tron->DX);CHKERRQ(info);
126:   info = TaoSetVariableBounds(tao,tron->XL,tron->XU);CHKERRQ(info);

128:   info = X->GetDimension(&tron->n); CHKERRQ(info);
129:   
130:   info = X->CreateIndexSet(&tron->Free_Local); CHKERRQ(info);
131:   info = tron->Free_Local->Duplicate(&tron->TT); CHKERRQ(info);

133:   TIS=tron->Free_Local;
134:   info = tron->H->CreateReducedMatrix(TIS,TIS,&tron->Hsub); CHKERRQ(info);
135:   info = TaoCreateLinearSolver(tao,HH,220,0); CHKERRQ(info);

137:   info = TaoCheckFGH(tao);CHKERRQ(info);

139:   TaoFunctionReturn(0);
140: }



146: static int TaoSolve_TRON(TAO_SOLVER tao, void*solver){

148:   TAO_TRON *tron = (TAO_TRON *)solver;
149:   int info;
150:   TaoInt lsflag,iter=0;
151:   TaoTerminateReason reason;
152:   TaoTruth optimal_face=TAO_FALSE,success;
153:   double prered,actred,delta,f,f_new,f_full,rhok,gnorm,gdx,xdiff,stepsize;
154:   TaoVec *XU, *XL;
155:   TaoVec *X,  *G;
156:   TaoVec *PG=tron->PG;
157:   TaoVec *R=tron->R, *DXFree=tron->DXFree;
158:   TaoVec *X_New=tron->X_New, *G_New=tron->G_New;
159:   TaoVec *DX=tron->DX, *Work=tron->Work;
160:   TaoMat *H, *Hsub=tron->Hsub;
161:   TaoIndexSet *Free_Local = tron->Free_Local, *TIS=tron->TT;

163:   TaoFunctionBegin;

165:   // Get initial trust region radius
166:   info = TaoGetInitialTrustRegionRadius(tao, &tron->delta); CHKERRQ(info);
167:   if (tron->delta <= 0) {
168:     SETERRQ(1, "Initial trust region radius must be positive");
169:   }

171:   // Get vectors we will need
172:   info = TaoGetSolution(tao, &X); CHKERRQ(info);
173:   info = TaoGetGradient(tao, &G); CHKERRQ(info);
174:   info = TaoGetHessian(tao, &H); CHKERRQ(info);
175:   info = TaoGetVariableBounds(tao, &XL, &XU); CHKERRQ(info);

177:   // Check that upper bound greater than lower bound
178:   info = TaoEvaluateVariableBounds(tao, XL, XU); CHKERRQ(info);

180:   tron->pgstepsize=1.0;

182:   /*   Project the current point onto the feasible set */
183:   info = X->Median(XL,X,XU); CHKERRQ(info);
184:   
185:   info = TaoComputeMeritFunctionGradient(tao,X,&tron->f,G);CHKERRQ(info);
186:   info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
187:   
188:   /* Project the gradient and calculate the norm */
189:   //  info = G_New->CopyFrom(G);CHKERRQ(info);
190:   info = PG->BoundGradientProjection(G,XL,X,XU);CHKERRQ(info);
191:   info = PG->Norm2(&tron->gnorm); CHKERRQ(info);

193:   if (tron->delta <= 0.0){
194:     tron->delta=TaoMax(tron->gnorm*tron->gnorm,1.0);
195:     //    tron->delta = TAO_INFINITY;
196:   }

198:   tron->stepsize=tron->delta;

200:   info = TaoMonitor(tao,iter++,tron->f,tron->gnorm,0.0,tron->delta,&reason);
201:   CHKERRQ(info);

203:   while (reason==TAO_CONTINUE_ITERATING){
204:     
205:     info = TaoGradProjections(tao,tron); CHKERRQ(info);

207:     info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
208:     info = Free_Local->GetSize(&tron->n_free); CHKERRQ(info);
209:     f=tron->f; delta=tron->delta; gnorm=tron->gnorm; 

211:     if (tron->n_free > 0){
212:       
213:       info = TaoComputeHessian(tao,X,H);CHKERRQ(info);

215:       /* Create a reduced linear system */
216:       info = R->SetReducedVec(G,Free_Local);CHKERRQ(info);
217:       info = R->Negate(); CHKERRQ(info);
218:       info = DXFree->SetReducedVec(DX,Free_Local);CHKERRQ(info);
219:       info = DXFree->SetToZero(); CHKERRQ(info);

221:       info = Hsub->SetReducedMatrix(H,Free_Local,Free_Local);CHKERRQ(info);

223:       info = TaoPreLinearSolve(tao,Hsub);CHKERRQ(info);

225:       while (1) {

227:          /* Approximately solve the reduced linear system */

229:         info = TaoLinearSolveTrustRegion(tao, Hsub, R, DXFree, delta, &success); CHKERRQ(info);

231:         info=DX->SetToZero(); CHKERRQ(info);
232:         info=DX->ReducedXPY(DXFree,Free_Local);CHKERRQ(info);

234:         info = G->Dot(DX,&gdx); CHKERRQ(info);
235:         info = PetscInfo1(tao,"Expected decrease in function value: %14.12e\n",gdx); CHKERRQ(info);

237:         stepsize=1.0; f_new=f;
238:         info = X_New->CopyFrom(X); CHKERRQ(info);
239:         info = G_New->CopyFrom(G); CHKERRQ(info);
240:         
241:         info = TaoLineSearchApply(tao,X_New,G_New,DX,Work,
242:                                   &f_new,&f_full,&stepsize,&lsflag);
243:         CHKERRQ(info);
244:         info = H->Multiply(DX,Work); CHKERRQ(info);
245:         info = Work->Aypx(0.5,G); CHKERRQ(info);
246:         info = Work->Dot(DX,&prered); CHKERRQ(info);
247:         actred = f_new - f;
248:         
249:         if (actred<0) rhok=TaoAbsScalar(-actred/prered);
250:         else rhok=0.0;

252:         /* Compare actual improvement to the quadratic model */
253:         if (rhok > tron->eta1) { /* Accept the point */

255:           info = DX->Waxpby(1.0,X_New,-1.0, X); CHKERRQ(info);
256:           info = DX->Norm2(&xdiff); CHKERRQ(info);
257:           xdiff*=stepsize;

259:           /* Adjust trust region size */
260:           if (rhok < tron->eta2 ){
261:             delta = TaoMin(xdiff,delta)*tron->sigma1;
262:           } else if (rhok > tron->eta4 ){
263:             delta= TaoMin(xdiff,delta)*tron->sigma3;
264:           } else if (rhok > tron->eta3 ){
265:             delta=TaoMin(xdiff,delta)*tron->sigma2;
266:           }

268:           info =  PG->BoundGradientProjection(G_New,XL,X_New,XU);
269:           CHKERRQ(info);
270:           info = PG->Norm2(&gnorm);  CHKERRQ(info);
271:           info = TronCheckOptimalFace(X_New,XL,XU,G_New,PG, Free_Local, TIS,
272:                                       &optimal_face); CHKERRQ(info);          
273:           if (stepsize < 1 || optimal_face==TAO_FALSE || reason!=TAO_CONTINUE_ITERATING ){
274:             f=f_new;
275:             info = X->CopyFrom(X_New); CHKERRQ(info);
276:             info = G->CopyFrom(G_New); CHKERRQ(info);
277:             break;
278:           }
279:           if (delta<=1e-30){
280:             break;
281:           }
282:         } 
283:         else if (delta <= 1e-30) {
284:           break;
285:         }
286:         else {
287:           delta /= 4.0;
288:         }
289:       } /* end linear solve loop */
290:       
291:     } else {
292:       
293:       actred=0;
294:       info =  Work->BoundGradientProjection(G,XL,X,XU);
295:       CHKERRQ(info);
296:       info = Work->Norm2(&gnorm);  CHKERRQ(info);
297:       /* if there were no free variables, no cg method */

299:     }

301:     tron->f=f;tron->gnorm=gnorm; tron->actred=actred; tron->delta=delta;
302:     info = TaoMonitor(tao,iter,f,gnorm,0.0,delta,&reason); CHKERRQ(info);
303:     if (reason!=TAO_CONTINUE_ITERATING) break;
304:     iter++;
305:     
306:   }  /* END MAIN LOOP  */

308:   TaoFunctionReturn(0);
309: }


314: static int TaoGradProjections(TAO_SOLVER tao,TAO_TRON *tron)
315: {
316:   int info;
317:   TaoInt lsflag=0,i;
318:   TaoTruth sameface=TAO_FALSE;
319:   double actred=-1.0,actred_max=0.0;
320:   double f_new, f_full;
321:   TaoVec *DX=tron->DX,*XL=tron->XL,*XU=tron->XU,*Work=tron->Work;
322:   TaoVec *X=tron->X,*G=tron->G;
323:   TaoIndexSet *TT1=tron->Free_Local, *TT2=tron->TT, *TT3;
324:   /*
325:      The gradient and function value passed into and out of this
326:      routine should be current and correct.
327:      
328:      The free, active, and binding variables should be already identified
329:   */
330:   
331:   TaoFunctionBegin;
332:   
333:   info = TaoGetSolution(tao,&X);CHKERRQ(info);
334:   info = TaoGetGradient(tao,&G);CHKERRQ(info);
335:   info = TaoGetVariableBounds(tao,&XL,&XU);CHKERRQ(info);

337:   info = TT1->WhichBetween(XL,X,XU); CHKERRQ(info);

339:   for (i=0;i<tron->maxgpits;i++){

341:     if ( -actred <= (tron->pg_ftol)*actred_max) break;
342:   
343:     tron->gp_iterates++; tron->total_gp_its++;      
344:     f_new=tron->f;

346:     info = DX->ScaleCopyFrom(-1.0,G); CHKERRQ(info);

348:     info = TaoLineSearchApply(tao,X,G,DX,Work,
349:                               &f_new,&f_full,&tron->pgstepsize,&lsflag);
350:     CHKERRQ(info);

352:     /* Update the iterate */
353:     actred = f_new - tron->f;
354:     actred_max = TaoMax(actred_max,-(f_new - tron->f));
355:     tron->f = f_new;

357:     info = TT2->WhichBetween(XL,X,XU); CHKERRQ(info);
358:     info = TT2->IsSame(TT1,&sameface);  CHKERRQ(info);
359:     if (sameface==TAO_TRUE) {
360:       break;
361:     } else {
362:       //      info = TT1->WhichBetween(XL,X,XU); CHKERRQ(info);
363:       TT3=TT2;
364:       TT2=TT1;
365:       TT1=TT3;
366:     }

368:   }
369:   
370:   TaoFunctionReturn(0);
371: }


376: static int TronCheckOptimalFace(TaoVec *X,TaoVec *XL,TaoVec*XU,TaoVec *PG,TaoVec*W,
377:                                 TaoIndexSet*Free_Local, TaoIndexSet*TT,
378:                                 TaoTruth *optimal)
379: {
380:   int info;
381:   TaoInt n_free;
382:   double rr;
383:   TaoTruth same;

385:   TaoFunctionBegin;
386:   *optimal = TAO_FALSE;

388:   /* Also check to see if the active set is the same */

390:   info = TT->WhichBetween(XL,X,XU); CHKERRQ(info);
391:   info = Free_Local->IsSame(TT,&same); CHKERRQ(info);
392:   info = Free_Local->GetSize(&n_free); CHKERRQ(info);
393:   if (same == TAO_FALSE){
394:     info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
395:     *optimal = TAO_FALSE;
396:     TaoFunctionReturn(0);
397:   } else {
398:     *optimal = TAO_TRUE;
399:   }

401:   info = W->CopyFrom(PG); CHKERRQ(info);
402:   info = W->Negate(); CHKERRQ(info);

404:   info = W->BoundGradientProjection(W,XL,X,XU); CHKERRQ(info);
405:   info = W->Axpy(1.0,PG); CHKERRQ(info);

407:   info = W->Norm2(&rr); CHKERRQ(info);
408:   if (rr>0) *optimal = TAO_FALSE;

410:   *optimal = TAO_FALSE;
411:   /*
412:     info = tron->TT->whichNonNegative(W); CHKERRQ(info);
413:     info = tron->TT->GetSize(&n); CHKERRQ(info);
414:     if (n==0) *optimal = TAO_TRUE;
415:   */
416:   TaoFunctionReturn(0);
417: }



423: int TaoDefaultMonitor_TRON(TAO_SOLVER tao,void *dummy)
424: {
425:   int info;
426:   TaoInt its,nfree,nbind;
427:   double fct,gnorm;
428:   TAO_TRON *tron;

430:   TaoFunctionBegin;
431:   info = TaoGetSolutionStatus(tao,&its,&fct,&gnorm,0,0,0);CHKERRQ(info);
432:   info = TaoGetSolverContext(tao,"tao_tron",(void**)&tron); CHKERRQ(info);
433:   if (tron){
434:     nfree=tron->n_free;
435:     nbind=tron->n_bind;
436:     info=TaoPrintInt(tao,"iter = %d,",its); CHKERRQ(info);
437:     info=TaoPrintDouble(tao," Function value: %g,",fct); CHKERRQ(info);
438:     info=TaoPrintDouble(tao,"  Residual: %g \n",gnorm);CHKERRQ(info);
439:     
440:     info=TaoPrintInt(tao," free vars = %d,",nfree); CHKERRQ(info);
441:     info=TaoPrintInt(tao," binding vars = %d\n",nbind); CHKERRQ(info);
442:   }
443:   TaoFunctionReturn(0);
444: }

446: int TaoSetMaxGPIts(TAO_SOLVER tao, int its){
447:   int info;
448:   TAO_TRON  *tron;

450:   TaoFunctionBegin;

452:   info = TaoGetSolverContext(tao,"tao_tron",(void**)&tron); CHKERRQ(info);
453:   if (tron){
454:     tron->maxgpits     = its;
455:   }
456:   TaoFunctionReturn(0);
457: }


462: static int TaoGetDualVariables_TRON(TAO_SOLVER tao, TaoVec* DXL, TaoVec* DXU, void *solver){

464:   TAO_TRON *tron = (TAO_TRON *) solver;
465:   TaoVec  *G=tron->G,*GP=tron->Work;
466:   TaoVec  *X,*XL,*XU;
467:   int       info;

469:   TaoFunctionBegin;
470:   info = TaoGetSolution(tao,&X); CHKERRQ(info);
471:   info = TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
472:   info = GP->BoundGradientProjection(G,XL,X,XU); CHKERRQ(info);

474:   info = DXL->Waxpby(-1.0,G,1.0,GP); CHKERRQ(info);
475:   info = DXU->SetToZero(); CHKERRQ(info);
476:   info = DXL->PointwiseMaximum(DXL,DXU); CHKERRQ(info);

478:   info = DXU->Waxpby(-1.0,GP,1.0,G); CHKERRQ(info);
479:   info = GP->SetToZero(); CHKERRQ(info);
480:   info = DXU->PointwiseMinimum(GP,DXU); CHKERRQ(info);

482:   TaoFunctionReturn(0);
483: }

485: /*------------------------------------------------------------*/
489: int TaoCreate_TRON(TAO_SOLVER tao)
490: {
491:   TAO_TRON *tron;
492:   int      info;

494:   TaoFunctionBegin;

496:   info = TaoNew(TAO_TRON,&tron); CHKERRQ(info);
497:   info = PetscLogObjectMemory(tao,sizeof(TAO_TRON)); CHKERRQ(info);

499:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_TRON,(void*)tron); CHKERRQ(info);
500:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_TRON,TaoSetDown_TRON); CHKERRQ(info);
501:   info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_TRON); CHKERRQ(info);
502:   info=TaoSetTaoViewRoutine(tao,TaoView_TRON); CHKERRQ(info);
503:   info=TaoSetTaoDualVariablesRoutine(tao,TaoGetDualVariables_TRON); CHKERRQ(info);

505:   info = TaoSetMaximumIterates(tao, 50); CHKERRQ(info);
506:   info = TaoSetTolerances(tao, 1e-10, 1e-10, 0, 0); CHKERRQ(info);

508:   info = TaoSetTrustRegionRadius(tao, 1.0); CHKERRQ(info);
509:   info = TaoSetTrustRegionTolerance(tao, 1.0e-12); CHKERRQ(info);

511:   /* Initialize pointers and variables */
512:   tron->n            = 0;
513:   tron->delta        = -1.0;
514:   tron->maxgpits     = 3;
515:   tron->pg_ftol      = 0.001;

517:   tron->eta1         = 1.0e-4;
518:   tron->eta2         = 0.25;
519:   tron->eta3         = 0.50;
520:   tron->eta4         = 0.90;

522:   tron->sigma1       = 0.5;
523:   tron->sigma2       = 2.0;
524:   tron->sigma3       = 4.0;

526:   tron->gp_iterates  = 0; /* Cumulative number */
527:   tron->cgits        = 0; /* Current iteration */
528:   tron->total_gp_its = 0;
529:   tron->cg_iterates  = 0;
530:   tron->total_cgits  = 0;
531:  
532:   tron->n_bind       = 0;
533:   tron->n_free       = 0;
534:   tron->n_upper      = 0;
535:   tron->n_lower      = 0;

537:   tron->DX=0;
538:   tron->DXFree=0;
539:   tron->R=0;
540:   tron->X_New=0;
541:   tron->G_New=0;
542:   tron->Work=0;
543:   tron->Free_Local=0;
544:   tron->TT=0;
545:   tron->Hsub=0;

547:   info = TaoCreateMoreThuenteBoundLineSearch(tao, 0 , 0); CHKERRQ(info);
548:   TaoFunctionReturn(0);
549: }