Actual source code: tron.c
1: /*$Id: tron.c 1.143 05/05/10 15:21:37-05:00 sarich@zorak.(none) $*/
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);
34: info = TaoDestroyLinearSolver(tao);CHKERRQ(info);
36: info = TaoSetLagrangianGradientVector(tao,0);CHKERRQ(info);
37: info = TaoSetStepDirectionVector(tao,0);CHKERRQ(info);
38: info = TaoSetVariableBounds(tao,0,0);CHKERRQ(info);
40: TaoFunctionReturn(0);
41: }
43: /*------------------------------------------------------------*/
46: static int TaoSetOptions_TRON(TAO_SOLVER tao, void*solver)
47: {
48: TAO_TRON *tron = (TAO_TRON *)solver;
49: int info,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);
65: /* info = TaoSetLinearSolverOptions(tao);CHKERRQ(info); */
67: TaoFunctionReturn(0);
68: }
70: /*------------------------------------------------------------*/
73: static int TaoView_TRON(TAO_SOLVER tao,void*solver)
74: {
75: TAO_TRON *tron = (TAO_TRON *)solver;
76: int info;
78: TaoFunctionBegin;
79: /*
80: info = TaoPrintf1(tao," Variables, Total: %d,",tron->n);
81: info = TaoPrintf3(tao,"Free: %d, Binding: %d \n",
82: tron->n_free, tron->n - tron->n_free,
83: tron->n_bind);CHKERRQ(info);
84: info = TaoPrintf1(tao," Equal lower bound: %d,",
85: tron->n_lower);CHKERRQ(info);
86: info = TaoPrintf1(tao," Equal upper bound: %d \n",
87: tron->n_upper);CHKERRQ(info);
88: */
89: info = TaoPrintInt(tao," Total PG its: %d,",tron->total_gp_its);CHKERRQ(info);
90: info = TaoPrintDouble(tao," PG tolerance: %4.3f \n",tron->pg_ftol);CHKERRQ(info);
91: info = TaoLineSearchView(tao);CHKERRQ(info);
92: info = TaoPrintStatement(tao," Linear Solver minimizes quadratic over Trust Region: \n");CHKERRQ(info);
94: TaoFunctionReturn(0);
95: }
98: /* ---------------------------------------------------------- */
101: static int TaoSetUp_TRON(TAO_SOLVER tao, void*solver){
103: int info;
104: TAO_TRON *tron = (TAO_TRON *)solver;
105: TaoVec* X;
106: TaoMat *HH;
107: TaoIndexSet *TIS;
108: TaoLinearSolver *ksp;
110: TaoFunctionBegin;
111: info = TaoGetSolution(tao,&tron->X);CHKERRQ(info); X=tron->X;
112: info = TaoGetHessian(tao,&tron->H);CHKERRQ(info); HH=tron->H;
114: /* Allocate some arrays */
115: info = X->Clone(&tron->DX); CHKERRQ(info);
116: info = X->Clone(&tron->X_New); CHKERRQ(info);
117: info = X->Clone(&tron->G_New); CHKERRQ(info);
118: info = X->Clone(&tron->Work); CHKERRQ(info);
119: info = X->Clone(&tron->DXFree); CHKERRQ(info);
120: info = X->Clone(&tron->R); CHKERRQ(info);
121: info = X->Clone(&tron->G); CHKERRQ(info);
122: info = X->Clone(&tron->PG); CHKERRQ(info);
123: info = X->Clone(&tron->XL); CHKERRQ(info);
124: info = X->Clone(&tron->XU); CHKERRQ(info);
126: info = TaoSetLagrangianGradientVector(tao,tron->PG);CHKERRQ(info);
127: info = TaoSetStepDirectionVector(tao,tron->DX);CHKERRQ(info);
128: info = TaoSetVariableBounds(tao,tron->XL,tron->XU);CHKERRQ(info);
130: info = X->GetDimension(&tron->n); CHKERRQ(info);
131:
132: info = X->CreateIndexSet(&tron->Free_Local); CHKERRQ(info);
133: info = tron->Free_Local->Duplicate(&tron->TT); CHKERRQ(info);
135: TIS=tron->Free_Local;
136: info = tron->H->CreateReducedMatrix(TIS,TIS,&tron->Hsub); CHKERRQ(info);
137: info = TaoCreateLinearSolver(tao,HH,220,0); CHKERRQ(info);
138: info = TaoGetLinearSolver(tao,&ksp); CHKERRQ(info);
139: // info = TaoCreateLinearSolver(tao,tron->Hsub,220,&tron->kspMarker); CHKERRQ(info);
141: info = ksp->SetTolerances(1e-2,1e-30,1e30,tron->n);CHKERRQ(info);
143: info = TaoCheckFGH(tao);CHKERRQ(info);
145: TaoFunctionReturn(0);
146: }
152: static int TaoSolve_TRON(TAO_SOLVER tao, void*solver){
154: TAO_TRON *tron = (TAO_TRON *)solver;
155: int info,lsflag,iter=0;
156: TaoTerminateReason reason;
157: TaoTruth optimal_face=TAO_FALSE,success;
158: double prered,actred,delta,f,f_new,rhok,gnorm,gdx,xdiff,stepsize;
159: TaoVec *XU, *XL;
160: TaoVec *X, *G;
161: TaoVec *PG=tron->PG;
162: TaoVec *R=tron->R, *DXFree=tron->DXFree;
163: TaoVec *X_New=tron->X_New, *G_New=tron->G_New;
164: TaoVec *DX=tron->DX, *Work=tron->Work;
165: TaoMat *H, *Hsub=tron->Hsub;
166: TaoIndexSet *Free_Local = tron->Free_Local, *TIS=tron->TT;
168: TaoFunctionBegin;
170: /* Check if upper bound greater than lower bound. */
171: info = TaoGetSolution(tao,&X);CHKERRQ(info); tron->X=X;
172: info = TaoGetGradient(tao,&G);CHKERRQ(info);
173: info = TaoGetVariableBounds(tao,&XL,&XU);CHKERRQ(info);
174: info = TaoEvaluateVariableBounds(tao,XL,XU); CHKERRQ(info);
175: info = TaoGetHessian(tao,&H);CHKERRQ(info); tron->H=H;
177: info = TaoGetInitialTrustRegionRadius(tao,&tron->delta); CHKERRQ(info);
178: tron->pgstepsize=1.0;
180: /* Project the current point onto the feasible set */
181: info = X->Median(XL,X,XU); CHKERRQ(info);
182:
183: info = TaoComputeMeritFunctionGradient(tao,X,&tron->f,G);CHKERRQ(info);
184: info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
185:
186: /* Project the gradient and calculate the norm */
187: // info = G_New->CopyFrom(G);CHKERRQ(info);
188: info = PG->BoundGradientProjection(G,XL,X,XU);CHKERRQ(info);
189: info = PG->Norm2(&tron->gnorm); CHKERRQ(info);
191: if (tron->delta <= 0.0){
192: tron->delta=TaoMax(tron->gnorm*tron->gnorm,1.0);
193: // tron->delta = TAO_INFINITY;
194: }
196: tron->stepsize=tron->delta;
198: info = TaoMonitor(tao,iter++,tron->f,tron->gnorm,0.0,tron->delta,&reason);
199: CHKERRQ(info);
201: while (reason==TAO_CONTINUE_ITERATING){
202:
203: info = TaoGradProjections(tao,tron); CHKERRQ(info);
205: info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
206: info = Free_Local->GetSize(&tron->n_free); CHKERRQ(info);
207: f=tron->f; delta=tron->delta; gnorm=tron->gnorm;
209: if (tron->n_free > 0){
210:
211: info = TaoComputeHessian(tao,X,H);CHKERRQ(info);
213: /* Create a reduced linear system */
214: info = R->SetReducedVec(G,Free_Local);CHKERRQ(info);
215: info = R->Negate(); CHKERRQ(info);
216: info = DXFree->SetReducedVec(DX,Free_Local);CHKERRQ(info);
217: info = DXFree->SetToZero(); CHKERRQ(info);
219: info = Hsub->SetReducedMatrix(H,Free_Local,Free_Local);CHKERRQ(info);
221: info = TaoPreLinearSolve(tao,Hsub);CHKERRQ(info);
223: while (1){
225: /* Approximately solve the reduced linear system */
226: // info = Hsub->minQuadraticTrustRegion(R,DXFree,delta,&success); CHKERRQ(info);
227: info = TaoMinQPTrust(tao, Hsub, R, DXFree, delta, &success); CHKERRQ(info);
229: info=DX->SetToZero(); CHKERRQ(info);
230: info=DX->ReducedXPY(DXFree,Free_Local);CHKERRQ(info);
232: info = G->Dot(DX,&gdx); CHKERRQ(info);
233: info = PetscLogInfo((tao,"Expected decrease in function value: %14.12e\n",gdx)); CHKERRQ(info);
235: stepsize=1.0; f_new=f;
236: info = X_New->CopyFrom(X); CHKERRQ(info);
237: info = G_New->CopyFrom(G); CHKERRQ(info);
238:
239: info = TaoLineSearchApply(tao,X_New,G_New,DX,Work,
240: &f_new,&stepsize,&lsflag);
241: CHKERRQ(info);
242: info = H->Multiply(DX,Work); CHKERRQ(info);
243: info = Work->Aypx(0.5,G); CHKERRQ(info);
244: info = Work->Dot(DX,&prered); CHKERRQ(info);
245: actred = f_new - f;
246:
247: if (actred<0) rhok=TaoAbsScalar(-actred/prered);
248: else rhok=0.0;
250: /* Compare actual improvement to the quadratic model */
251: if (rhok > tron->eta1) { /* Accept the point */
253: info = DX->Waxpby(1.0,X_New,-1.0, X); CHKERRQ(info);
254: info = DX->Norm2(&xdiff); CHKERRQ(info);
255: xdiff*=stepsize;
257: /* Adjust trust region size */
258: if (rhok < tron->eta2 ){
259: delta = TaoMin(xdiff,delta)*tron->sigma1;
260: } else if (rhok > tron->eta4 ){
261: delta= TaoMin(xdiff,delta)*tron->sigma3;
262: } else if (rhok > tron->eta3 ){
263: delta=TaoMin(xdiff,delta)*tron->sigma2;
264: }
266: info = PG->BoundGradientProjection(G_New,XL,X_New,XU);
267: CHKERRQ(info);
268: info = PG->Norm2(&gnorm); CHKERRQ(info);
269:
270: info = TaoMonitor(tao,iter,f_new,gnorm,0.0,delta,&reason);CHKERRQ(info);
272: info = TronCheckOptimalFace(X_New,XL,XU,G_New,PG, Free_Local, TIS,
273: &optimal_face); CHKERRQ(info);
275: if (stepsize < 1 || optimal_face==TAO_FALSE || reason!=TAO_CONTINUE_ITERATING ){
276: f=f_new;
277: info = X->CopyFrom(X_New); CHKERRQ(info);
278: info = G->CopyFrom(G_New); CHKERRQ(info);
279: break;
280: }
281: } else {
282: delta /= 4.0;
283: info = TaoMonitor(tao,iter,f,gnorm,0.0,delta,&reason);CHKERRQ(info);
284: if (reason!=TAO_CONTINUE_ITERATING ){ break; }
285: }
287: } /* end linear solve loop */
288:
289:
290: } else {
291:
292: actred=0;
293: info = Work->BoundGradientProjection(G,XL,X,XU);
294: CHKERRQ(info);
295: info = Work->Norm2(&gnorm); CHKERRQ(info);
296: /* if there were no free variables, no cg method */
297: info = TaoMonitor(tao,iter,f,gnorm,0.0,delta,&reason); CHKERRQ(info);
299: }
301: tron->f=f;tron->gnorm=gnorm; tron->actred=actred; tron->delta=delta;
302: if (reason!=TAO_CONTINUE_ITERATING) break;
303: iter++;
304:
305: } /* END MAIN LOOP */
307: TaoFunctionReturn(0);
308: }
313: static int TaoGradProjections(TAO_SOLVER tao,TAO_TRON *tron)
314: {
315: int i,lsflag=0,info;
316: TaoTruth sameface=TAO_FALSE;
317: double actred=-1.0,actred_max=0.0;
318: double f_new;
319: TaoVec *DX=tron->DX,*XL=tron->XL,*XU=tron->XU,*Work=tron->Work;
320: TaoVec *X=tron->X,*G=tron->G;
321: TaoIndexSet *TT1=tron->Free_Local, *TT2=tron->TT, *TT3;
322: /*
323: The gradient and function value passed into and out of this
324: routine should be current and correct.
325:
326: The free, active, and binding variables should be already identified
327: */
328:
329: TaoFunctionBegin;
330:
331: info = TaoGetSolution(tao,&X);CHKERRQ(info);
332: info = TaoGetGradient(tao,&G);CHKERRQ(info);
333: info = TaoGetVariableBounds(tao,&XL,&XU);CHKERRQ(info);
335: info = TT1->WhichBetween(XL,X,XU); CHKERRQ(info);
337: for (i=0;i<tron->maxgpits;i++){
339: if ( -actred <= (tron->pg_ftol)*actred_max) break;
340:
341: tron->gp_iterates++; tron->total_gp_its++;
342: f_new=tron->f;
344: info = DX->ScaleCopyFrom(-1.0,G); CHKERRQ(info);
346: info = TaoLineSearchApply(tao,X,G,DX,Work,
347: &f_new,&tron->pgstepsize,&lsflag);
348: CHKERRQ(info);
350: /* Update the iterate */
351: actred = f_new - tron->f;
352: actred_max = TaoMax(actred_max,-(f_new - tron->f));
353: tron->f = f_new;
355: info = TT2->WhichBetween(XL,X,XU); CHKERRQ(info);
356: info = TT2->IsSame(TT1,&sameface); CHKERRQ(info);
357: if (sameface==TAO_TRUE) {
358: break;
359: } else {
360: // info = TT1->WhichBetween(XL,X,XU); CHKERRQ(info);
361: TT3=TT2;
362: TT2=TT1;
363: TT1=TT3;
364: }
366: }
367:
368: TaoFunctionReturn(0);
369: }
374: static int TronCheckOptimalFace(TaoVec *X,TaoVec *XL,TaoVec*XU,TaoVec *PG,TaoVec*W,
375: TaoIndexSet*Free_Local, TaoIndexSet*TT,
376: TaoTruth *optimal)
377: {
378: int info,n_free;
379: double rr;
380: TaoTruth same;
382: TaoFunctionBegin;
383: *optimal = TAO_FALSE;
385: /* Also check to see if the active set is the same */
387: info = TT->WhichBetween(XL,X,XU); CHKERRQ(info);
388: info = Free_Local->IsSame(TT,&same); CHKERRQ(info);
389: info = Free_Local->GetSize(&n_free); CHKERRQ(info);
390: if (same == TAO_FALSE){
391: info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
392: *optimal = TAO_FALSE;
393: TaoFunctionReturn(0);
394: } else {
395: *optimal = TAO_TRUE;
396: }
398: info = W->CopyFrom(PG); CHKERRQ(info);
399: info = W->Negate(); CHKERRQ(info);
401: info = W->BoundGradientProjection(W,XL,X,XU); CHKERRQ(info);
402: info = W->Axpy(1.0,PG); CHKERRQ(info);
404: info = W->Norm2(&rr); CHKERRQ(info);
405: if (rr>0) *optimal = TAO_FALSE;
407: *optimal = TAO_FALSE;
408: /*
409: info = tron->TT->whichNonNegative(W); CHKERRQ(info);
410: info = tron->TT->GetSize(&n); CHKERRQ(info);
411: if (n==0) *optimal = TAO_TRUE;
412: */
413: TaoFunctionReturn(0);
414: }
420: int TaoDefaultMonitor_TRON(TAO_SOLVER tao,void *dummy)
421: {
422: int info;
423: int its,nfree,nbind;
424: double fct,gnorm;
425: TAO_TRON *tron;
427: TaoFunctionBegin;
428: info = TaoGetSolutionStatus(tao,&its,&fct,&gnorm,0,0,0);CHKERRQ(info);
429: info = TaoGetSolverContext(tao,"tao_tron",(void**)&tron); CHKERRQ(info);
430: if (tron){
431: nfree=tron->n_free;
432: nbind=tron->n_bind;
433: info=TaoPrintInt(tao,"iter = %d,",its); CHKERRQ(info);
434: info=TaoPrintDouble(tao," Function value: %g,",fct); CHKERRQ(info);
435: info=TaoPrintDouble(tao," Residual: %g \n",gnorm);CHKERRQ(info);
436:
437: info=TaoPrintInt(tao," free vars = %d,",nfree); CHKERRQ(info);
438: info=TaoPrintInt(tao," binding vars = %d\n",nbind); CHKERRQ(info);
439: }
440: TaoFunctionReturn(0);
441: }
443: int TaoSetMaxGPIts(TAO_SOLVER tao, int its){
444: int info;
445: TAO_TRON *tron;
447: TaoFunctionBegin;
449: info = TaoGetSolverContext(tao,"tao_tron",(void**)&tron); CHKERRQ(info);
450: if (tron){
451: tron->maxgpits = its;
452: }
453: TaoFunctionReturn(0);
454: }
459: static int TaoGetDualVariables_TRON(TAO_SOLVER tao, TaoVec* DXL, TaoVec* DXU, void *solver){
461: TAO_TRON *tron = (TAO_TRON *) solver;
462: TaoVec *G=tron->G,*GP=tron->Work;
463: TaoVec *X,*XL,*XU;
464: int info;
466: TaoFunctionBegin;
467: info = TaoGetSolution(tao,&X); CHKERRQ(info);
468: info = TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
469: info = GP->BoundGradientProjection(G,XL,X,XU); CHKERRQ(info);
471: info = DXL->Waxpby(-1.0,G,1.0,GP); CHKERRQ(info);
472: info = DXU->SetToZero(); CHKERRQ(info);
473: info = DXL->PointwiseMaximum(DXL,DXU); CHKERRQ(info);
475: info = DXU->Waxpby(-1.0,GP,1.0,G); CHKERRQ(info);
476: info = GP->SetToZero(); CHKERRQ(info);
477: info = DXU->PointwiseMinimum(GP,DXU); CHKERRQ(info);
479: TaoFunctionReturn(0);
480: }
482: /*------------------------------------------------------------*/
486: int TaoCreate_TRON(TAO_SOLVER tao)
487: {
488: TAO_TRON *tron;
489: int info;
491: TaoFunctionBegin;
493: info = TaoNew(TAO_TRON,&tron); CHKERRQ(info);
494: info = PetscLogObjectMemory(tao,sizeof(TAO_TRON)); CHKERRQ(info);
496: info=TaoSetTaoSolveRoutine(tao,TaoSolve_TRON,(void*)tron); CHKERRQ(info);
497: info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_TRON,TaoSetDown_TRON); CHKERRQ(info);
498: info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_TRON); CHKERRQ(info);
499: info=TaoSetTaoViewRoutine(tao,TaoView_TRON); CHKERRQ(info);
500: info=TaoSetTaoDualVariablesRoutine(tao,TaoGetDualVariables_TRON); CHKERRQ(info);
502: info = TaoSetMaximumIterates(tao,500); CHKERRQ(info);
503: info = TaoSetTolerances(tao,1e-12,1e-12,0,0); CHKERRQ(info);
504: /*
505: tao->defaultmonitor = TaoDefaultMonitor_TRON;
506: */
507: info = TaoSetTrustRegionTolerance(tao,1.0e-12);CHKERRQ(info);
508: // info = TaoSetTrustRegionRadius(tao,-1.0);CHKERRQ(info);
509: info = TaoSetTrustRegionRadius(tao,1);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;
516: tron->eta1 = 0.001;
517: tron->eta2 = 0.3;
518: tron->eta3 = 0.5;
519: tron->eta4 = 0.9;
520: tron->sigma1 = 0.5;
521: tron->sigma2 = 2.0;
522: tron->sigma3 = 4.0;
524: tron->gp_iterates = 0; /* Cumulative number */
525: tron->cgits = 0; /* Current iteration */
526: tron->total_gp_its = 0;
527: tron->cg_iterates = 0;
528: tron->total_cgits = 0;
529:
530: tron->n_bind = 0;
531: tron->n_free = 0;
532: tron->n_upper = 0;
533: tron->n_lower = 0;
535: tron->DX=0;
536: tron->DXFree=0;
537: tron->R=0;
538: tron->X_New=0;
539: tron->G_New=0;
540: tron->Work=0;
541: tron->Free_Local=0;
542: tron->TT=0;
543: tron->Hsub=0;
545: // info = TaoCreateProjectedLineSearch(tao); CHKERRQ(info);
546: info = TaoCreateMoreThuenteBoundLineSearch(tao,0,0.9); CHKERRQ(info);
548: TaoFunctionReturn(0);
549: }