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