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