Actual source code: ntl.c

  1: #include <../src/tao/unconstrained/impls/ntl/ntlimpl.h>

  3: #include <petscksp.h>

  5: #define NTL_INIT_CONSTANT      0
  6: #define NTL_INIT_DIRECTION     1
  7: #define NTL_INIT_INTERPOLATION 2
  8: #define NTL_INIT_TYPES         3

 10: #define NTL_UPDATE_REDUCTION     0
 11: #define NTL_UPDATE_INTERPOLATION 1
 12: #define NTL_UPDATE_TYPES         2

 14: static const char *NTL_INIT[64] = {"constant", "direction", "interpolation"};

 16: static const char *NTL_UPDATE[64] = {"reduction", "interpolation"};

 18: /* Implements Newton's Method with a trust-region, line-search approach for
 19:    solving unconstrained minimization problems.  A More'-Thuente line search
 20:    is used to guarantee that the bfgs preconditioner remains positive
 21:    definite. */

 23: #define NTL_NEWTON          0
 24: #define NTL_BFGS            1
 25: #define NTL_SCALED_GRADIENT 2
 26: #define NTL_GRADIENT        3

 28: static PetscErrorCode TaoSolve_NTL(Tao tao)
 29: {
 30:   TAO_NTL                     *tl = (TAO_NTL *)tao->data;
 31:   KSPType                      ksp_type;
 32:   PetscBool                    is_nash, is_stcg, is_gltr, is_bfgs, is_jacobi, is_symmetric, sym_set;
 33:   KSPConvergedReason           ksp_reason;
 34:   PC                           pc;
 35:   TaoLineSearchConvergedReason ls_reason;

 37:   PetscReal fmin, ftrial, prered, actred, kappa, sigma;
 38:   PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
 39:   PetscReal f, fold, gdx, gnorm;
 40:   PetscReal step = 1.0;

 42:   PetscReal norm_d = 0.0;
 43:   PetscInt  stepType;
 44:   PetscInt  its;

 46:   PetscInt bfgsUpdates = 0;
 47:   PetscInt needH;

 49:   PetscInt i_max = 5;
 50:   PetscInt j_max = 1;
 51:   PetscInt i, j, n, N;

 53:   PetscInt tr_reject;

 55:   PetscFunctionBegin;
 56:   if (tao->XL || tao->XU || tao->ops->computebounds) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by ntl algorithm\n"));

 58:   PetscCall(KSPGetType(tao->ksp, &ksp_type));
 59:   PetscCall(PetscStrcmp(ksp_type, KSPNASH, &is_nash));
 60:   PetscCall(PetscStrcmp(ksp_type, KSPSTCG, &is_stcg));
 61:   PetscCall(PetscStrcmp(ksp_type, KSPGLTR, &is_gltr));
 62:   PetscCheck(is_nash || is_stcg || is_gltr, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "TAO_NTR requires nash, stcg, or gltr for the KSP");

 64:   /* Initialize the radius and modify if it is too large or small */
 65:   tao->trust = tao->trust0;
 66:   tao->trust = PetscMax(tao->trust, tl->min_radius);
 67:   tao->trust = PetscMin(tao->trust, tl->max_radius);

 69:   /* Allocate the vectors needed for the BFGS approximation */
 70:   PetscCall(KSPGetPC(tao->ksp, &pc));
 71:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs));
 72:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi));
 73:   if (is_bfgs) {
 74:     tl->bfgs_pre = pc;
 75:     PetscCall(PCLMVMGetMatLMVM(tl->bfgs_pre, &tl->M));
 76:     PetscCall(VecGetLocalSize(tao->solution, &n));
 77:     PetscCall(VecGetSize(tao->solution, &N));
 78:     PetscCall(MatSetSizes(tl->M, n, n, N, N));
 79:     PetscCall(MatLMVMAllocate(tl->M, tao->solution, tao->gradient));
 80:     PetscCall(MatIsSymmetricKnown(tl->M, &sym_set, &is_symmetric));
 81:     PetscCheck(sym_set && is_symmetric, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric.");
 82:   } else if (is_jacobi) PetscCall(PCJacobiSetUseAbs(pc, PETSC_TRUE));

 84:   /* Check convergence criteria */
 85:   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
 86:   PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
 87:   PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
 88:   needH = 1;

 90:   tao->reason = TAO_CONTINUE_ITERATING;
 91:   PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
 92:   PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
 93:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
 94:   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);

 96:   /* Initialize trust-region radius */
 97:   switch (tl->init_type) {
 98:   case NTL_INIT_CONSTANT:
 99:     /* Use the initial radius specified */
100:     break;

102:   case NTL_INIT_INTERPOLATION:
103:     /* Use the initial radius specified */
104:     max_radius = 0.0;

106:     for (j = 0; j < j_max; ++j) {
107:       fmin  = f;
108:       sigma = 0.0;

110:       if (needH) {
111:         PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
112:         needH = 0;
113:       }

115:       for (i = 0; i < i_max; ++i) {
116:         PetscCall(VecCopy(tao->solution, tl->W));
117:         PetscCall(VecAXPY(tl->W, -tao->trust / gnorm, tao->gradient));

119:         PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));
120:         if (PetscIsInfOrNanReal(ftrial)) {
121:           tau = tl->gamma1_i;
122:         } else {
123:           if (ftrial < fmin) {
124:             fmin  = ftrial;
125:             sigma = -tao->trust / gnorm;
126:           }

128:           PetscCall(MatMult(tao->hessian, tao->gradient, tao->stepdirection));
129:           PetscCall(VecDot(tao->gradient, tao->stepdirection, &prered));

131:           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
132:           actred = f - ftrial;
133:           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
134:             kappa = 1.0;
135:           } else {
136:             kappa = actred / prered;
137:           }

139:           tau_1   = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust + (1.0 - tl->theta_i) * prered - actred);
140:           tau_2   = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust - (1.0 + tl->theta_i) * prered + actred);
141:           tau_min = PetscMin(tau_1, tau_2);
142:           tau_max = PetscMax(tau_1, tau_2);

144:           if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tl->mu1_i) {
145:             /* Great agreement */
146:             max_radius = PetscMax(max_radius, tao->trust);

148:             if (tau_max < 1.0) {
149:               tau = tl->gamma3_i;
150:             } else if (tau_max > tl->gamma4_i) {
151:               tau = tl->gamma4_i;
152:             } else if (tau_1 >= 1.0 && tau_1 <= tl->gamma4_i && tau_2 < 1.0) {
153:               tau = tau_1;
154:             } else if (tau_2 >= 1.0 && tau_2 <= tl->gamma4_i && tau_1 < 1.0) {
155:               tau = tau_2;
156:             } else {
157:               tau = tau_max;
158:             }
159:           } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tl->mu2_i) {
160:             /* Good agreement */
161:             max_radius = PetscMax(max_radius, tao->trust);

163:             if (tau_max < tl->gamma2_i) {
164:               tau = tl->gamma2_i;
165:             } else if (tau_max > tl->gamma3_i) {
166:               tau = tl->gamma3_i;
167:             } else {
168:               tau = tau_max;
169:             }
170:           } else {
171:             /* Not good agreement */
172:             if (tau_min > 1.0) {
173:               tau = tl->gamma2_i;
174:             } else if (tau_max < tl->gamma1_i) {
175:               tau = tl->gamma1_i;
176:             } else if ((tau_min < tl->gamma1_i) && (tau_max >= 1.0)) {
177:               tau = tl->gamma1_i;
178:             } else if ((tau_1 >= tl->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1_i) || (tau_2 >= 1.0))) {
179:               tau = tau_1;
180:             } else if ((tau_2 >= tl->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1_i) || (tau_2 >= 1.0))) {
181:               tau = tau_2;
182:             } else {
183:               tau = tau_max;
184:             }
185:           }
186:         }
187:         tao->trust = tau * tao->trust;
188:       }

190:       if (fmin < f) {
191:         f = fmin;
192:         PetscCall(VecAXPY(tao->solution, sigma, tao->gradient));
193:         PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));

195:         PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
196:         PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
197:         needH = 1;

199:         PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
200:         PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
201:         PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
202:         if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
203:       }
204:     }
205:     tao->trust = PetscMax(tao->trust, max_radius);

207:     /* Modify the radius if it is too large or small */
208:     tao->trust = PetscMax(tao->trust, tl->min_radius);
209:     tao->trust = PetscMin(tao->trust, tl->max_radius);
210:     break;

212:   default:
213:     /* Norm of the first direction will initialize radius */
214:     tao->trust = 0.0;
215:     break;
216:   }

218:   /* Set counter for gradient/reset steps */
219:   tl->ntrust = 0;
220:   tl->newt   = 0;
221:   tl->bfgs   = 0;
222:   tl->grad   = 0;

224:   /* Have not converged; continue with Newton method */
225:   while (tao->reason == TAO_CONTINUE_ITERATING) {
226:     /* Call general purpose update function */
227:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
228:     ++tao->niter;
229:     tao->ksp_its = 0;
230:     /* Compute the Hessian */
231:     if (needH) PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));

233:     if (tl->bfgs_pre) {
234:       /* Update the limited memory preconditioner */
235:       PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
236:       ++bfgsUpdates;
237:     }
238:     PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre));
239:     /* Solve the Newton system of equations */
240:     PetscCall(KSPCGSetRadius(tao->ksp, tl->max_radius));
241:     PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
242:     PetscCall(KSPGetIterationNumber(tao->ksp, &its));
243:     tao->ksp_its += its;
244:     tao->ksp_tot_its += its;
245:     PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));

247:     if (0.0 == tao->trust) {
248:       /* Radius was uninitialized; use the norm of the direction */
249:       if (norm_d > 0.0) {
250:         tao->trust = norm_d;

252:         /* Modify the radius if it is too large or small */
253:         tao->trust = PetscMax(tao->trust, tl->min_radius);
254:         tao->trust = PetscMin(tao->trust, tl->max_radius);
255:       } else {
256:         /* The direction was bad; set radius to default value and re-solve
257:            the trust-region subproblem to get a direction */
258:         tao->trust = tao->trust0;

260:         /* Modify the radius if it is too large or small */
261:         tao->trust = PetscMax(tao->trust, tl->min_radius);
262:         tao->trust = PetscMin(tao->trust, tl->max_radius);

264:         PetscCall(KSPCGSetRadius(tao->ksp, tl->max_radius));
265:         PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
266:         PetscCall(KSPGetIterationNumber(tao->ksp, &its));
267:         tao->ksp_its += its;
268:         tao->ksp_tot_its += its;
269:         PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));

271:         PetscCheck(norm_d != 0.0, PetscObjectComm((PetscObject)tao), PETSC_ERR_PLIB, "Initial direction zero");
272:       }
273:     }

275:     PetscCall(VecScale(tao->stepdirection, -1.0));
276:     PetscCall(KSPGetConvergedReason(tao->ksp, &ksp_reason));
277:     if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tl->bfgs_pre)) {
278:       /* Preconditioner is numerically indefinite; reset the
279:          approximate if using BFGS preconditioning. */
280:       PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
281:       PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
282:       bfgsUpdates = 1;
283:     }

285:     /* Check trust-region reduction conditions */
286:     tr_reject = 0;
287:     if (NTL_UPDATE_REDUCTION == tl->update_type) {
288:       /* Get predicted reduction */
289:       PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
290:       if (prered >= 0.0) {
291:         /* The predicted reduction has the wrong sign.  This cannot
292:            happen in infinite precision arithmetic.  Step should
293:            be rejected! */
294:         tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
295:         tr_reject  = 1;
296:       } else {
297:         /* Compute trial step and function value */
298:         PetscCall(VecCopy(tao->solution, tl->W));
299:         PetscCall(VecAXPY(tl->W, 1.0, tao->stepdirection));
300:         PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));

302:         if (PetscIsInfOrNanReal(ftrial)) {
303:           tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
304:           tr_reject  = 1;
305:         } else {
306:           /* Compute and actual reduction */
307:           actred = f - ftrial;
308:           prered = -prered;
309:           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
310:             kappa = 1.0;
311:           } else {
312:             kappa = actred / prered;
313:           }

315:           /* Accept of reject the step and update radius */
316:           if (kappa < tl->eta1) {
317:             /* Reject the step */
318:             tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
319:             tr_reject  = 1;
320:           } else {
321:             /* Accept the step */
322:             if (kappa < tl->eta2) {
323:               /* Marginal bad step */
324:               tao->trust = tl->alpha2 * PetscMin(tao->trust, norm_d);
325:             } else if (kappa < tl->eta3) {
326:               /* Reasonable step */
327:               tao->trust = tl->alpha3 * tao->trust;
328:             } else if (kappa < tl->eta4) {
329:               /* Good step */
330:               tao->trust = PetscMax(tl->alpha4 * norm_d, tao->trust);
331:             } else {
332:               /* Very good step */
333:               tao->trust = PetscMax(tl->alpha5 * norm_d, tao->trust);
334:             }
335:           }
336:         }
337:       }
338:     } else {
339:       /* Get predicted reduction */
340:       PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
341:       if (prered >= 0.0) {
342:         /* The predicted reduction has the wrong sign.  This cannot
343:            happen in infinite precision arithmetic.  Step should
344:            be rejected! */
345:         tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
346:         tr_reject  = 1;
347:       } else {
348:         PetscCall(VecCopy(tao->solution, tl->W));
349:         PetscCall(VecAXPY(tl->W, 1.0, tao->stepdirection));
350:         PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));
351:         if (PetscIsInfOrNanReal(ftrial)) {
352:           tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
353:           tr_reject  = 1;
354:         } else {
355:           PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx));

357:           actred = f - ftrial;
358:           prered = -prered;
359:           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
360:             kappa = 1.0;
361:           } else {
362:             kappa = actred / prered;
363:           }

365:           tau_1   = tl->theta * gdx / (tl->theta * gdx - (1.0 - tl->theta) * prered + actred);
366:           tau_2   = tl->theta * gdx / (tl->theta * gdx + (1.0 + tl->theta) * prered - actred);
367:           tau_min = PetscMin(tau_1, tau_2);
368:           tau_max = PetscMax(tau_1, tau_2);

370:           if (kappa >= 1.0 - tl->mu1) {
371:             /* Great agreement; accept step and update radius */
372:             if (tau_max < 1.0) {
373:               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
374:             } else if (tau_max > tl->gamma4) {
375:               tao->trust = PetscMax(tao->trust, tl->gamma4 * norm_d);
376:             } else {
377:               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
378:             }
379:           } else if (kappa >= 1.0 - tl->mu2) {
380:             /* Good agreement */

382:             if (tau_max < tl->gamma2) {
383:               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
384:             } else if (tau_max > tl->gamma3) {
385:               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
386:             } else if (tau_max < 1.0) {
387:               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
388:             } else {
389:               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
390:             }
391:           } else {
392:             /* Not good agreement */
393:             if (tau_min > 1.0) {
394:               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
395:             } else if (tau_max < tl->gamma1) {
396:               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
397:             } else if ((tau_min < tl->gamma1) && (tau_max >= 1.0)) {
398:               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
399:             } else if ((tau_1 >= tl->gamma1) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1) || (tau_2 >= 1.0))) {
400:               tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
401:             } else if ((tau_2 >= tl->gamma1) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1) || (tau_2 >= 1.0))) {
402:               tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
403:             } else {
404:               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
405:             }
406:             tr_reject = 1;
407:           }
408:         }
409:       }
410:     }

412:     if (tr_reject) {
413:       /* The trust-region constraints rejected the step.  Apply a linesearch.
414:          Check for descent direction. */
415:       PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
416:       if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
417:         /* Newton step is not descent or direction produced Inf or NaN */

419:         if (!tl->bfgs_pre) {
420:           /* We don't have the bfgs matrix around and updated
421:              Must use gradient direction in this case */
422:           PetscCall(VecCopy(tao->gradient, tao->stepdirection));
423:           PetscCall(VecScale(tao->stepdirection, -1.0));
424:           ++tl->grad;
425:           stepType = NTL_GRADIENT;
426:         } else {
427:           /* Attempt to use the BFGS direction */
428:           PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
429:           PetscCall(VecScale(tao->stepdirection, -1.0));

431:           /* Check for success (descent direction) */
432:           PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
433:           if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
434:             /* BFGS direction is not descent or direction produced not a number
435:                We can assert bfgsUpdates > 1 in this case because
436:                the first solve produces the scaled gradient direction,
437:                which is guaranteed to be descent */

439:             /* Use steepest descent direction (scaled) */
440:             PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
441:             PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
442:             PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
443:             PetscCall(VecScale(tao->stepdirection, -1.0));

445:             bfgsUpdates = 1;
446:             ++tl->grad;
447:             stepType = NTL_GRADIENT;
448:           } else {
449:             PetscCall(MatLMVMGetUpdateCount(tl->M, &bfgsUpdates));
450:             if (1 == bfgsUpdates) {
451:               /* The first BFGS direction is always the scaled gradient */
452:               ++tl->grad;
453:               stepType = NTL_GRADIENT;
454:             } else {
455:               ++tl->bfgs;
456:               stepType = NTL_BFGS;
457:             }
458:           }
459:         }
460:       } else {
461:         /* Computed Newton step is descent */
462:         ++tl->newt;
463:         stepType = NTL_NEWTON;
464:       }

466:       /* Perform the linesearch */
467:       fold = f;
468:       PetscCall(VecCopy(tao->solution, tl->Xold));
469:       PetscCall(VecCopy(tao->gradient, tl->Gold));

471:       step = 1.0;
472:       PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason));
473:       PetscCall(TaoAddLineSearchCounts(tao));

475:       while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NTL_GRADIENT) { /* Linesearch failed */
476:         /* Linesearch failed */
477:         f = fold;
478:         PetscCall(VecCopy(tl->Xold, tao->solution));
479:         PetscCall(VecCopy(tl->Gold, tao->gradient));

481:         switch (stepType) {
482:         case NTL_NEWTON:
483:           /* Failed to obtain acceptable iterate with Newton step */

485:           if (tl->bfgs_pre) {
486:             /* We don't have the bfgs matrix around and being updated
487:                Must use gradient direction in this case */
488:             PetscCall(VecCopy(tao->gradient, tao->stepdirection));
489:             ++tl->grad;
490:             stepType = NTL_GRADIENT;
491:           } else {
492:             /* Attempt to use the BFGS direction */
493:             PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));

495:             /* Check for success (descent direction) */
496:             PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
497:             if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
498:               /* BFGS direction is not descent or direction produced
499:                  not a number.  We can assert bfgsUpdates > 1 in this case
500:                  Use steepest descent direction (scaled) */
501:               PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
502:               PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
503:               PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));

505:               bfgsUpdates = 1;
506:               ++tl->grad;
507:               stepType = NTL_GRADIENT;
508:             } else {
509:               PetscCall(MatLMVMGetUpdateCount(tl->M, &bfgsUpdates));
510:               if (1 == bfgsUpdates) {
511:                 /* The first BFGS direction is always the scaled gradient */
512:                 ++tl->grad;
513:                 stepType = NTL_GRADIENT;
514:               } else {
515:                 ++tl->bfgs;
516:                 stepType = NTL_BFGS;
517:               }
518:             }
519:           }
520:           break;

522:         case NTL_BFGS:
523:           /* Can only enter if pc_type == NTL_PC_BFGS
524:              Failed to obtain acceptable iterate with BFGS step
525:              Attempt to use the scaled gradient direction */
526:           PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
527:           PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
528:           PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));

530:           bfgsUpdates = 1;
531:           ++tl->grad;
532:           stepType = NTL_GRADIENT;
533:           break;
534:         }
535:         PetscCall(VecScale(tao->stepdirection, -1.0));

537:         /* This may be incorrect; linesearch has values for stepmax and stepmin
538:            that should be reset. */
539:         step = 1.0;
540:         PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason));
541:         PetscCall(TaoAddLineSearchCounts(tao));
542:       }

544:       if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) {
545:         /* Failed to find an improving point */
546:         f = fold;
547:         PetscCall(VecCopy(tl->Xold, tao->solution));
548:         PetscCall(VecCopy(tl->Gold, tao->gradient));
549:         tao->trust  = 0.0;
550:         step        = 0.0;
551:         tao->reason = TAO_DIVERGED_LS_FAILURE;
552:         break;
553:       } else if (stepType == NTL_NEWTON) {
554:         if (step < tl->nu1) {
555:           /* Very bad step taken; reduce radius */
556:           tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
557:         } else if (step < tl->nu2) {
558:           /* Reasonably bad step taken; reduce radius */
559:           tao->trust = tl->omega2 * PetscMin(norm_d, tao->trust);
560:         } else if (step < tl->nu3) {
561:           /* Reasonable step was taken; leave radius alone */
562:           if (tl->omega3 < 1.0) {
563:             tao->trust = tl->omega3 * PetscMin(norm_d, tao->trust);
564:           } else if (tl->omega3 > 1.0) {
565:             tao->trust = PetscMax(tl->omega3 * norm_d, tao->trust);
566:           }
567:         } else if (step < tl->nu4) {
568:           /* Full step taken; increase the radius */
569:           tao->trust = PetscMax(tl->omega4 * norm_d, tao->trust);
570:         } else {
571:           /* More than full step taken; increase the radius */
572:           tao->trust = PetscMax(tl->omega5 * norm_d, tao->trust);
573:         }
574:       } else {
575:         /* Newton step was not good; reduce the radius */
576:         tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
577:       }
578:     } else {
579:       /* Trust-region step is accepted */
580:       PetscCall(VecCopy(tl->W, tao->solution));
581:       f = ftrial;
582:       PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));
583:       ++tl->ntrust;
584:     }

586:     /* The radius may have been increased; modify if it is too large */
587:     tao->trust = PetscMin(tao->trust, tl->max_radius);

589:     /* Check for converged */
590:     PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
591:     PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number");
592:     needH = 1;

594:     PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
595:     PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
596:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
597:   }
598:   PetscFunctionReturn(PETSC_SUCCESS);
599: }

601: /* ---------------------------------------------------------- */
602: static PetscErrorCode TaoSetUp_NTL(Tao tao)
603: {
604:   TAO_NTL *tl = (TAO_NTL *)tao->data;

606:   PetscFunctionBegin;
607:   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
608:   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
609:   if (!tl->W) PetscCall(VecDuplicate(tao->solution, &tl->W));
610:   if (!tl->Xold) PetscCall(VecDuplicate(tao->solution, &tl->Xold));
611:   if (!tl->Gold) PetscCall(VecDuplicate(tao->solution, &tl->Gold));
612:   tl->bfgs_pre = NULL;
613:   tl->M        = NULL;
614:   PetscFunctionReturn(PETSC_SUCCESS);
615: }

617: /*------------------------------------------------------------*/
618: static PetscErrorCode TaoDestroy_NTL(Tao tao)
619: {
620:   TAO_NTL *tl = (TAO_NTL *)tao->data;

622:   PetscFunctionBegin;
623:   if (tao->setupcalled) {
624:     PetscCall(VecDestroy(&tl->W));
625:     PetscCall(VecDestroy(&tl->Xold));
626:     PetscCall(VecDestroy(&tl->Gold));
627:   }
628:   PetscCall(KSPDestroy(&tao->ksp));
629:   PetscCall(PetscFree(tao->data));
630:   PetscFunctionReturn(PETSC_SUCCESS);
631: }

633: /*------------------------------------------------------------*/
634: static PetscErrorCode TaoSetFromOptions_NTL(Tao tao, PetscOptionItems *PetscOptionsObject)
635: {
636:   TAO_NTL *tl = (TAO_NTL *)tao->data;

638:   PetscFunctionBegin;
639:   PetscOptionsHeadBegin(PetscOptionsObject, "Newton trust region with line search method for unconstrained optimization");
640:   PetscCall(PetscOptionsEList("-tao_ntl_init_type", "radius initialization type", "", NTL_INIT, NTL_INIT_TYPES, NTL_INIT[tl->init_type], &tl->init_type, NULL));
641:   PetscCall(PetscOptionsEList("-tao_ntl_update_type", "radius update type", "", NTL_UPDATE, NTL_UPDATE_TYPES, NTL_UPDATE[tl->update_type], &tl->update_type, NULL));
642:   PetscCall(PetscOptionsReal("-tao_ntl_eta1", "poor steplength; reduce radius", "", tl->eta1, &tl->eta1, NULL));
643:   PetscCall(PetscOptionsReal("-tao_ntl_eta2", "reasonable steplength; leave radius alone", "", tl->eta2, &tl->eta2, NULL));
644:   PetscCall(PetscOptionsReal("-tao_ntl_eta3", "good steplength; increase radius", "", tl->eta3, &tl->eta3, NULL));
645:   PetscCall(PetscOptionsReal("-tao_ntl_eta4", "excellent steplength; greatly increase radius", "", tl->eta4, &tl->eta4, NULL));
646:   PetscCall(PetscOptionsReal("-tao_ntl_alpha1", "", "", tl->alpha1, &tl->alpha1, NULL));
647:   PetscCall(PetscOptionsReal("-tao_ntl_alpha2", "", "", tl->alpha2, &tl->alpha2, NULL));
648:   PetscCall(PetscOptionsReal("-tao_ntl_alpha3", "", "", tl->alpha3, &tl->alpha3, NULL));
649:   PetscCall(PetscOptionsReal("-tao_ntl_alpha4", "", "", tl->alpha4, &tl->alpha4, NULL));
650:   PetscCall(PetscOptionsReal("-tao_ntl_alpha5", "", "", tl->alpha5, &tl->alpha5, NULL));
651:   PetscCall(PetscOptionsReal("-tao_ntl_nu1", "poor steplength; reduce radius", "", tl->nu1, &tl->nu1, NULL));
652:   PetscCall(PetscOptionsReal("-tao_ntl_nu2", "reasonable steplength; leave radius alone", "", tl->nu2, &tl->nu2, NULL));
653:   PetscCall(PetscOptionsReal("-tao_ntl_nu3", "good steplength; increase radius", "", tl->nu3, &tl->nu3, NULL));
654:   PetscCall(PetscOptionsReal("-tao_ntl_nu4", "excellent steplength; greatly increase radius", "", tl->nu4, &tl->nu4, NULL));
655:   PetscCall(PetscOptionsReal("-tao_ntl_omega1", "", "", tl->omega1, &tl->omega1, NULL));
656:   PetscCall(PetscOptionsReal("-tao_ntl_omega2", "", "", tl->omega2, &tl->omega2, NULL));
657:   PetscCall(PetscOptionsReal("-tao_ntl_omega3", "", "", tl->omega3, &tl->omega3, NULL));
658:   PetscCall(PetscOptionsReal("-tao_ntl_omega4", "", "", tl->omega4, &tl->omega4, NULL));
659:   PetscCall(PetscOptionsReal("-tao_ntl_omega5", "", "", tl->omega5, &tl->omega5, NULL));
660:   PetscCall(PetscOptionsReal("-tao_ntl_mu1_i", "", "", tl->mu1_i, &tl->mu1_i, NULL));
661:   PetscCall(PetscOptionsReal("-tao_ntl_mu2_i", "", "", tl->mu2_i, &tl->mu2_i, NULL));
662:   PetscCall(PetscOptionsReal("-tao_ntl_gamma1_i", "", "", tl->gamma1_i, &tl->gamma1_i, NULL));
663:   PetscCall(PetscOptionsReal("-tao_ntl_gamma2_i", "", "", tl->gamma2_i, &tl->gamma2_i, NULL));
664:   PetscCall(PetscOptionsReal("-tao_ntl_gamma3_i", "", "", tl->gamma3_i, &tl->gamma3_i, NULL));
665:   PetscCall(PetscOptionsReal("-tao_ntl_gamma4_i", "", "", tl->gamma4_i, &tl->gamma4_i, NULL));
666:   PetscCall(PetscOptionsReal("-tao_ntl_theta_i", "", "", tl->theta_i, &tl->theta_i, NULL));
667:   PetscCall(PetscOptionsReal("-tao_ntl_mu1", "", "", tl->mu1, &tl->mu1, NULL));
668:   PetscCall(PetscOptionsReal("-tao_ntl_mu2", "", "", tl->mu2, &tl->mu2, NULL));
669:   PetscCall(PetscOptionsReal("-tao_ntl_gamma1", "", "", tl->gamma1, &tl->gamma1, NULL));
670:   PetscCall(PetscOptionsReal("-tao_ntl_gamma2", "", "", tl->gamma2, &tl->gamma2, NULL));
671:   PetscCall(PetscOptionsReal("-tao_ntl_gamma3", "", "", tl->gamma3, &tl->gamma3, NULL));
672:   PetscCall(PetscOptionsReal("-tao_ntl_gamma4", "", "", tl->gamma4, &tl->gamma4, NULL));
673:   PetscCall(PetscOptionsReal("-tao_ntl_theta", "", "", tl->theta, &tl->theta, NULL));
674:   PetscCall(PetscOptionsReal("-tao_ntl_min_radius", "lower bound on initial radius", "", tl->min_radius, &tl->min_radius, NULL));
675:   PetscCall(PetscOptionsReal("-tao_ntl_max_radius", "upper bound on radius", "", tl->max_radius, &tl->max_radius, NULL));
676:   PetscCall(PetscOptionsReal("-tao_ntl_epsilon", "tolerance used when computing actual and predicted reduction", "", tl->epsilon, &tl->epsilon, NULL));
677:   PetscOptionsHeadEnd();
678:   PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
679:   PetscCall(KSPSetFromOptions(tao->ksp));
680:   PetscFunctionReturn(PETSC_SUCCESS);
681: }

683: /*------------------------------------------------------------*/
684: static PetscErrorCode TaoView_NTL(Tao tao, PetscViewer viewer)
685: {
686:   TAO_NTL  *tl = (TAO_NTL *)tao->data;
687:   PetscBool isascii;

689:   PetscFunctionBegin;
690:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
691:   if (isascii) {
692:     PetscCall(PetscViewerASCIIPushTab(viewer));
693:     PetscCall(PetscViewerASCIIPrintf(viewer, "Trust-region steps: %" PetscInt_FMT "\n", tl->ntrust));
694:     PetscCall(PetscViewerASCIIPrintf(viewer, "Newton search steps: %" PetscInt_FMT "\n", tl->newt));
695:     PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS search steps: %" PetscInt_FMT "\n", tl->bfgs));
696:     PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient search steps: %" PetscInt_FMT "\n", tl->grad));
697:     PetscCall(PetscViewerASCIIPopTab(viewer));
698:   }
699:   PetscFunctionReturn(PETSC_SUCCESS);
700: }

702: /* ---------------------------------------------------------- */
703: /*MC
704:   TAONTL - Newton's method with trust region globalization and line search fallback.
705:   At each iteration, the Newton trust region method solves the system for d
706:   and performs a line search in the d direction:

708:             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k

710:   Options Database Keys:
711: + -tao_ntl_init_type - "constant","direction","interpolation"
712: . -tao_ntl_update_type - "reduction","interpolation"
713: . -tao_ntl_min_radius - lower bound on trust region radius
714: . -tao_ntl_max_radius - upper bound on trust region radius
715: . -tao_ntl_epsilon - tolerance for accepting actual / predicted reduction
716: . -tao_ntl_mu1_i - mu1 interpolation init factor
717: . -tao_ntl_mu2_i - mu2 interpolation init factor
718: . -tao_ntl_gamma1_i - gamma1 interpolation init factor
719: . -tao_ntl_gamma2_i - gamma2 interpolation init factor
720: . -tao_ntl_gamma3_i - gamma3 interpolation init factor
721: . -tao_ntl_gamma4_i - gamma4 interpolation init factor
722: . -tao_ntl_theta_i - theta1 interpolation init factor
723: . -tao_ntl_eta1 - eta1 reduction update factor
724: . -tao_ntl_eta2 - eta2 reduction update factor
725: . -tao_ntl_eta3 - eta3 reduction update factor
726: . -tao_ntl_eta4 - eta4 reduction update factor
727: . -tao_ntl_alpha1 - alpha1 reduction update factor
728: . -tao_ntl_alpha2 - alpha2 reduction update factor
729: . -tao_ntl_alpha3 - alpha3 reduction update factor
730: . -tao_ntl_alpha4 - alpha4 reduction update factor
731: . -tao_ntl_alpha4 - alpha4 reduction update factor
732: . -tao_ntl_mu1 - mu1 interpolation update
733: . -tao_ntl_mu2 - mu2 interpolation update
734: . -tao_ntl_gamma1 - gamma1 interpolcation update
735: . -tao_ntl_gamma2 - gamma2 interpolcation update
736: . -tao_ntl_gamma3 - gamma3 interpolcation update
737: . -tao_ntl_gamma4 - gamma4 interpolation update
738: - -tao_ntl_theta - theta1 interpolation update

740:   Level: beginner
741: M*/
742: PETSC_EXTERN PetscErrorCode TaoCreate_NTL(Tao tao)
743: {
744:   TAO_NTL    *tl;
745:   const char *morethuente_type = TAOLINESEARCHMT;

747:   PetscFunctionBegin;
748:   PetscCall(PetscNew(&tl));
749:   tao->ops->setup          = TaoSetUp_NTL;
750:   tao->ops->solve          = TaoSolve_NTL;
751:   tao->ops->view           = TaoView_NTL;
752:   tao->ops->setfromoptions = TaoSetFromOptions_NTL;
753:   tao->ops->destroy        = TaoDestroy_NTL;

755:   /* Override default settings (unless already changed) */
756:   if (!tao->max_it_changed) tao->max_it = 50;
757:   if (!tao->trust0_changed) tao->trust0 = 100.0;

759:   tao->data = (void *)tl;

761:   /* Default values for trust-region radius update based on steplength */
762:   tl->nu1 = 0.25;
763:   tl->nu2 = 0.50;
764:   tl->nu3 = 1.00;
765:   tl->nu4 = 1.25;

767:   tl->omega1 = 0.25;
768:   tl->omega2 = 0.50;
769:   tl->omega3 = 1.00;
770:   tl->omega4 = 2.00;
771:   tl->omega5 = 4.00;

773:   /* Default values for trust-region radius update based on reduction */
774:   tl->eta1 = 1.0e-4;
775:   tl->eta2 = 0.25;
776:   tl->eta3 = 0.50;
777:   tl->eta4 = 0.90;

779:   tl->alpha1 = 0.25;
780:   tl->alpha2 = 0.50;
781:   tl->alpha3 = 1.00;
782:   tl->alpha4 = 2.00;
783:   tl->alpha5 = 4.00;

785:   /* Default values for trust-region radius update based on interpolation */
786:   tl->mu1 = 0.10;
787:   tl->mu2 = 0.50;

789:   tl->gamma1 = 0.25;
790:   tl->gamma2 = 0.50;
791:   tl->gamma3 = 2.00;
792:   tl->gamma4 = 4.00;

794:   tl->theta = 0.05;

796:   /* Default values for trust region initialization based on interpolation */
797:   tl->mu1_i = 0.35;
798:   tl->mu2_i = 0.50;

800:   tl->gamma1_i = 0.0625;
801:   tl->gamma2_i = 0.5;
802:   tl->gamma3_i = 2.0;
803:   tl->gamma4_i = 5.0;

805:   tl->theta_i = 0.25;

807:   /* Remaining parameters */
808:   tl->min_radius = 1.0e-10;
809:   tl->max_radius = 1.0e10;
810:   tl->epsilon    = 1.0e-6;

812:   tl->init_type   = NTL_INIT_INTERPOLATION;
813:   tl->update_type = NTL_UPDATE_REDUCTION;

815:   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
816:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
817:   PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
818:   PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
819:   PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
820:   PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
821:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
822:   PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
823:   PetscCall(KSPAppendOptionsPrefix(tao->ksp, "tao_ntl_"));
824:   PetscCall(KSPSetType(tao->ksp, KSPSTCG));
825:   PetscFunctionReturn(PETSC_SUCCESS);
826: }