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