xref: /petsc/src/tao/bound/impls/tron/tron.c (revision cf1aed2ce99d23e50336629af3ca8cf096900abb) 
1 #include <../src/tao/bound/impls/tron/tron.h>
2 #include <petsc-private/kspimpl.h>
3 #include <petsc-private/matimpl.h>
4 #include <../src/tao/matrix/submatfree.h>
5 
6 
7 /* TRON Routines */
8 static PetscErrorCode TronGradientProjections(Tao,TAO_TRON*);
9 /*------------------------------------------------------------*/
10 #undef __FUNCT__
11 #define __FUNCT__ "TaoDestroy_TRON"
12 static PetscErrorCode TaoDestroy_TRON(Tao tao)
13 {
14   TAO_TRON       *tron = (TAO_TRON *)tao->data;
15   PetscErrorCode ierr;
16 
17   PetscFunctionBegin;
18   ierr = VecDestroy(&tron->X_New);CHKERRQ(ierr);
19   ierr = VecDestroy(&tron->G_New);CHKERRQ(ierr);
20   ierr = VecDestroy(&tron->Work);CHKERRQ(ierr);
21   ierr = VecDestroy(&tron->DXFree);CHKERRQ(ierr);
22   ierr = VecDestroy(&tron->R);CHKERRQ(ierr);
23   ierr = VecDestroy(&tron->diag);CHKERRQ(ierr);
24   ierr = VecScatterDestroy(&tron->scatter);CHKERRQ(ierr);
25   ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr);
26   ierr = MatDestroy(&tron->H_sub);CHKERRQ(ierr);
27   ierr = MatDestroy(&tron->Hpre_sub);CHKERRQ(ierr);
28   ierr = PetscFree(tao->data);
29   PetscFunctionReturn(0);
30 }
31 
32 /*------------------------------------------------------------*/
33 #undef __FUNCT__
34 #define __FUNCT__ "TaoSetFromOptions_TRON"
35 static PetscErrorCode TaoSetFromOptions_TRON(Tao tao)
36 {
37   TAO_TRON       *tron = (TAO_TRON *)tao->data;
38   PetscErrorCode ierr;
39   PetscBool      flg;
40 
41   PetscFunctionBegin;
42   ierr = PetscOptionsHead("Newton Trust Region Method for bound constrained optimization");CHKERRQ(ierr);
43   ierr = PetscOptionsInt("-tao_tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg);CHKERRQ(ierr);
44   ierr = PetscOptionsTail();CHKERRQ(ierr);
45   ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
46   ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
47   PetscFunctionReturn(0);
48 }
49 
50 /*------------------------------------------------------------*/
51 #undef __FUNCT__
52 #define __FUNCT__ "TaoView_TRON"
53 static PetscErrorCode TaoView_TRON(Tao tao, PetscViewer viewer)
54 {
55   TAO_TRON         *tron = (TAO_TRON *)tao->data;
56   PetscBool        isascii;
57   PetscErrorCode   ierr;
58 
59   PetscFunctionBegin;
60   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
61   if (isascii) {
62     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
63     ierr = PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);CHKERRQ(ierr);
64     ierr = PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)tron->pg_ftol);CHKERRQ(ierr);
65     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
66   }
67   PetscFunctionReturn(0);
68 }
69 
70 
71 /* ---------------------------------------------------------- */
72 #undef __FUNCT__
73 #define __FUNCT__ "TaoSetup_TRON"
74 static PetscErrorCode TaoSetup_TRON(Tao tao)
75 {
76   PetscErrorCode ierr;
77   TAO_TRON       *tron = (TAO_TRON *)tao->data;
78 
79   PetscFunctionBegin;
80 
81   /* Allocate some arrays */
82   ierr = VecDuplicate(tao->solution, &tron->diag);CHKERRQ(ierr);
83   ierr = VecDuplicate(tao->solution, &tron->X_New);CHKERRQ(ierr);
84   ierr = VecDuplicate(tao->solution, &tron->G_New);CHKERRQ(ierr);
85   ierr = VecDuplicate(tao->solution, &tron->Work);CHKERRQ(ierr);
86   ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr);
87   ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr);
88   if (!tao->XL) {
89       ierr = VecDuplicate(tao->solution, &tao->XL);CHKERRQ(ierr);
90       ierr = VecSet(tao->XL, PETSC_NINFINITY);CHKERRQ(ierr);
91   }
92   if (!tao->XU) {
93       ierr = VecDuplicate(tao->solution, &tao->XU);CHKERRQ(ierr);
94       ierr = VecSet(tao->XU, PETSC_INFINITY);CHKERRQ(ierr);
95   }
96   PetscFunctionReturn(0);
97 }
98 
99 
100 
101 #undef __FUNCT__
102 #define __FUNCT__ "TaoSolve_TRON"
103 static PetscErrorCode TaoSolve_TRON(Tao tao)
104 {
105   TAO_TRON                     *tron = (TAO_TRON *)tao->data;
106   PetscErrorCode               ierr;
107   PetscInt                     iter=0,its;
108   TaoConvergedReason           reason = TAO_CONTINUE_ITERATING;
109   TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
110   PetscReal                    prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize;
111 
112   PetscFunctionBegin;
113   tron->pgstepsize=1.0;
114   tao->trust = tao->trust0;
115   /*   Project the current point onto the feasible set */
116   ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr);
117   ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr);
118   ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr);
119 
120   ierr = TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);CHKERRQ(ierr);
121   ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr);
122 
123   ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr);
124 
125   /* Project the gradient and calculate the norm */
126   ierr = VecBoundGradientProjection(tao->gradient,tao->solution, tao->XL, tao->XU, tao->gradient);CHKERRQ(ierr);
127   ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr);
128 
129   if (PetscIsInfOrNanReal(tron->f) || PetscIsInfOrNanReal(tron->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf pr NaN");
130   if (tao->trust <= 0) {
131     tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0);
132   }
133 
134   tron->stepsize=tao->trust;
135   ierr = TaoMonitor(tao, iter, tron->f, tron->gnorm, 0.0, tron->stepsize, &reason);CHKERRQ(ierr);
136   while (reason==TAO_CONTINUE_ITERATING){
137 
138     ierr = TronGradientProjections(tao,tron);CHKERRQ(ierr);
139     f=tron->f; delta=tao->trust;
140     tron->n_free_last = tron->n_free;
141     ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr);
142 
143     ierr = ISGetSize(tron->Free_Local, &tron->n_free);CHKERRQ(ierr);
144 
145     /* If no free variables */
146     if (tron->n_free == 0) {
147       actred=0;
148       PetscInfo(tao,"No free variables in tron iteration.");
149       break;
150 
151     }
152     /* use free_local to mask/submat gradient, hessian, stepdirection */
153     ierr = TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R);CHKERRQ(ierr);
154     ierr = TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree);CHKERRQ(ierr);
155     ierr = VecSet(tron->DXFree,0.0);CHKERRQ(ierr);
156     ierr = VecScale(tron->R, -1.0);CHKERRQ(ierr);
157     ierr = TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub);CHKERRQ(ierr);
158     if (tao->hessian == tao->hessian_pre) {
159       ierr = MatDestroy(&tron->Hpre_sub);CHKERRQ(ierr);
160       ierr = PetscObjectReference((PetscObject)(tron->H_sub));CHKERRQ(ierr);
161       tron->Hpre_sub = tron->H_sub;
162     } else {
163       ierr = TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub);CHKERRQ(ierr);
164     }
165     ierr = KSPReset(tao->ksp);CHKERRQ(ierr);
166     ierr = KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub);CHKERRQ(ierr);
167     while (1) {
168 
169       /* Approximately solve the reduced linear system */
170       ierr = KSPSTCGSetRadius(tao->ksp,delta);CHKERRQ(ierr);
171 
172       ierr = KSPSolve(tao->ksp, tron->R, tron->DXFree);CHKERRQ(ierr);
173       ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
174       tao->ksp_its+=its;
175       ierr = VecSet(tao->stepdirection,0.0);CHKERRQ(ierr);
176 
177       /* Add dxfree matrix to compute step direction vector */
178       ierr = VecISAXPY(tao->stepdirection,tron->Free_Local,1.0,tron->DXFree);CHKERRQ(ierr);
179       if (0) {
180         PetscReal rhs,stepnorm;
181         ierr = VecNorm(tron->R,NORM_2,&rhs);CHKERRQ(ierr);
182         ierr = VecNorm(tron->DXFree,NORM_2,&stepnorm);CHKERRQ(ierr);
183         ierr = PetscPrintf(PETSC_COMM_WORLD,"|rhs|=%g\t|s|=%g\n",(double)rhs,(double)stepnorm);CHKERRQ(ierr);
184       }
185 
186 
187       ierr = VecDot(tao->gradient, tao->stepdirection, &gdx);CHKERRQ(ierr);
188       ierr = PetscInfo1(tao,"Expected decrease in function value: %14.12e\n",(double)gdx);CHKERRQ(ierr);
189 
190       ierr = VecCopy(tao->solution, tron->X_New);CHKERRQ(ierr);
191       ierr = VecCopy(tao->gradient, tron->G_New);CHKERRQ(ierr);
192 
193       stepsize=1.0;f_new=f;
194 
195       ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);CHKERRQ(ierr);
196       ierr = TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);CHKERRQ(ierr);CHKERRQ(ierr);
197       ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
198 
199       ierr = MatMult(tao->hessian, tao->stepdirection, tron->Work);CHKERRQ(ierr);
200       ierr = VecAYPX(tron->Work, 0.5, tao->gradient);CHKERRQ(ierr);
201       ierr = VecDot(tao->stepdirection, tron->Work, &prered);CHKERRQ(ierr);
202       actred = f_new - f;
203       if (actred<0) {
204         rhok=PetscAbs(-actred/prered);
205       } else {
206         rhok=0.0;
207       }
208 
209       /* Compare actual improvement to the quadratic model */
210       if (rhok > tron->eta1) { /* Accept the point */
211         /* d = x_new - x */
212         ierr = VecCopy(tron->X_New, tao->stepdirection);CHKERRQ(ierr);
213         ierr = VecAXPY(tao->stepdirection, -1.0, tao->solution);CHKERRQ(ierr);
214 
215         ierr = VecNorm(tao->stepdirection, NORM_2, &xdiff);CHKERRQ(ierr);
216         xdiff *= stepsize;
217 
218         /* Adjust trust region size */
219         if (rhok < tron->eta2 ){
220           delta = PetscMin(xdiff,delta)*tron->sigma1;
221         } else if (rhok > tron->eta4 ){
222           delta= PetscMin(xdiff,delta)*tron->sigma3;
223         } else if (rhok > tron->eta3 ){
224           delta=PetscMin(xdiff,delta)*tron->sigma2;
225         }
226         ierr = VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient);CHKERRQ(ierr);
227         if (tron->Free_Local) {
228           ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr);
229         }
230         ierr = VecWhichBetween(tao->XL, tron->X_New, tao->XU, &tron->Free_Local);CHKERRQ(ierr);
231         f=f_new;
232         ierr = VecNorm(tao->gradient,NORM_2,&tron->gnorm);CHKERRQ(ierr);
233         ierr = VecCopy(tron->X_New, tao->solution);CHKERRQ(ierr);
234         ierr = VecCopy(tron->G_New, tao->gradient);CHKERRQ(ierr);
235         break;
236       }
237       else if (delta <= 1e-30) {
238         break;
239       }
240       else {
241         delta /= 4.0;
242       }
243     } /* end linear solve loop */
244 
245 
246     tron->f=f; tron->actred=actred; tao->trust=delta;
247     iter++;
248     ierr = TaoMonitor(tao, iter, tron->f, tron->gnorm, 0.0, delta, &reason);CHKERRQ(ierr);
249   }  /* END MAIN LOOP  */
250 
251   PetscFunctionReturn(0);
252 }
253 
254 
255 #undef __FUNCT__
256 #define __FUNCT__ "TronGradientProjections"
257 static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron)
258 {
259   PetscErrorCode                 ierr;
260   PetscInt                       i;
261   TaoLineSearchConvergedReason ls_reason;
262   PetscReal                      actred=-1.0,actred_max=0.0;
263   PetscReal                      f_new;
264   /*
265      The gradient and function value passed into and out of this
266      routine should be current and correct.
267 
268      The free, active, and binding variables should be already identified
269   */
270   PetscFunctionBegin;
271   if (tron->Free_Local) {
272     ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr);
273   }
274   ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr);
275 
276   for (i=0;i<tron->maxgpits;i++){
277 
278     if ( -actred <= (tron->pg_ftol)*actred_max) break;
279 
280     tron->gp_iterates++; tron->total_gp_its++;
281     f_new=tron->f;
282 
283     ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr);
284     ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
285     ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);CHKERRQ(ierr);
286     ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection,
287                               &tron->pgstepsize, &ls_reason);CHKERRQ(ierr);
288     ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
289 
290 
291     /* Update the iterate */
292     actred = f_new - tron->f;
293     actred_max = PetscMax(actred_max,-(f_new - tron->f));
294     tron->f = f_new;
295     if (tron->Free_Local) {
296       ierr = ISDestroy(&tron->Free_Local);CHKERRQ(ierr);
297     }
298     ierr = VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);CHKERRQ(ierr);
299   }
300 
301   PetscFunctionReturn(0);
302 }
303 
304 #undef __FUNCT__
305 #define __FUNCT__ "TaoComputeDual_TRON"
306 static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU) {
307 
308   TAO_TRON       *tron = (TAO_TRON *)tao->data;
309   PetscErrorCode ierr;
310 
311   PetscFunctionBegin;
312   PetscValidHeaderSpecific(tao,TAO_CLASSID,1);
313   PetscValidHeaderSpecific(DXL,VEC_CLASSID,2);
314   PetscValidHeaderSpecific(DXU,VEC_CLASSID,3);
315   if (!tron->Work || !tao->gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n");
316 
317   ierr = VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work);CHKERRQ(ierr);
318   ierr = VecCopy(tron->Work,DXL);CHKERRQ(ierr);
319   ierr = VecAXPY(DXL,-1.0,tao->gradient);CHKERRQ(ierr);
320   ierr = VecSet(DXU,0.0);CHKERRQ(ierr);
321   ierr = VecPointwiseMax(DXL,DXL,DXU);CHKERRQ(ierr);
322 
323   ierr = VecCopy(tao->gradient,DXU);CHKERRQ(ierr);
324   ierr = VecAXPY(DXU,-1.0,tron->Work);CHKERRQ(ierr);
325   ierr = VecSet(tron->Work,0.0);CHKERRQ(ierr);
326   ierr = VecPointwiseMin(DXU,tron->Work,DXU);CHKERRQ(ierr);
327   PetscFunctionReturn(0);
328 }
329 
330 /*------------------------------------------------------------*/
331 /*MC
332   TAOTRON - The TRON algorithm is an active-set Newton trust region method
333   for bound-constrained minimization.
334 
335   Options Database Keys:
336 + -tao_tron_maxgpits - maximum number of gradient projections per TRON iterate
337 - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets
338 
339   Level: beginner
340 M*/
341 EXTERN_C_BEGIN
342 #undef __FUNCT__
343 #define __FUNCT__ "TaoCreate_TRON"
344 PetscErrorCode TaoCreate_TRON(Tao tao)
345 {
346   TAO_TRON       *tron;
347   PetscErrorCode ierr;
348   const char     *morethuente_type = TAOLINESEARCHMT;
349 
350   PetscFunctionBegin;
351   tao->ops->setup = TaoSetup_TRON;
352   tao->ops->solve = TaoSolve_TRON;
353   tao->ops->view = TaoView_TRON;
354   tao->ops->setfromoptions = TaoSetFromOptions_TRON;
355   tao->ops->destroy = TaoDestroy_TRON;
356   tao->ops->computedual = TaoComputeDual_TRON;
357 
358   ierr = PetscNewLog(tao,&tron);CHKERRQ(ierr);
359 
360   tao->max_it = 50;
361 #if defined(PETSC_USE_REAL_SINGLE)
362   tao->fatol = 1e-5;
363   tao->frtol = 1e-5;
364   tao->steptol = 1e-6;
365 #else
366   tao->fatol = 1e-10;
367   tao->frtol = 1e-10;
368   tao->steptol = 1e-12;
369 #endif
370   tao->data = (void*)tron;
371   tao->trust0       = 1.0;
372 
373   /* Initialize pointers and variables */
374   tron->n            = 0;
375   tron->maxgpits     = 3;
376   tron->pg_ftol      = 0.001;
377 
378   tron->eta1         = 1.0e-4;
379   tron->eta2         = 0.25;
380   tron->eta3         = 0.50;
381   tron->eta4         = 0.90;
382 
383   tron->sigma1       = 0.5;
384   tron->sigma2       = 2.0;
385   tron->sigma3       = 4.0;
386 
387   tron->gp_iterates  = 0; /* Cumulative number */
388   tron->total_gp_its = 0;
389   tron->n_free       = 0;
390 
391   tron->DXFree=NULL;
392   tron->R=NULL;
393   tron->X_New=NULL;
394   tron->G_New=NULL;
395   tron->Work=NULL;
396   tron->Free_Local=NULL;
397   tron->H_sub=NULL;
398   tron->Hpre_sub=NULL;
399   tao->subset_type = TAO_SUBSET_SUBVEC;
400 
401   ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr);
402   ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type);CHKERRQ(ierr);
403   ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr);
404 
405   ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp);CHKERRQ(ierr);
406   ierr = KSPSetType(tao->ksp,KSPSTCG);CHKERRQ(ierr);
407   PetscFunctionReturn(0);
408 }
409 EXTERN_C_END
410