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