Actual source code: linesearchsecant.c

      1: #include <petsc/private/linesearchimpl.h>
      2: #include <petsc/private/snesimpl.h>
      3: #include <petscsnes.h>

      5: static PetscErrorCode SNESLineSearchApply_Secant(SNESLineSearch linesearch)
      6: {
      7:   PetscBool        changed_y, changed_w;
      8:   Vec              X;
      9:   Vec              F;
     10:   Vec              Y;
     11:   Vec              W;
     12:   SNES             snes;
     13:   PetscReal        gnorm;
     14:   PetscReal        ynorm;
     15:   PetscReal        xnorm;
     16:   PetscReal        minlambda, maxlambda, atol, ltol;
     17:   PetscViewer      monitor;
     18:   PetscReal        lambda, lambda_old, lambda_mid, lambda_update, delLambda;
     19:   PetscReal        fnrm, fnrm_old, fnrm_mid;
     20:   PetscReal        delFnrm, delFnrm_old, del2Fnrm;
     21:   PetscInt         i, max_it;
     22:   SNESObjectiveFn *objective;

     24:   PetscFunctionBegin;
     25:   PetscCall(SNESLineSearchGetVecs(linesearch, &X, &F, &Y, &W, NULL));
     26:   PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &gnorm, &ynorm));
     27:   PetscCall(SNESLineSearchGetLambda(linesearch, &lambda));
     28:   PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
     29:   PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED));
     30:   PetscCall(SNESLineSearchGetTolerances(linesearch, &minlambda, &maxlambda, NULL, &atol, &ltol, &max_it));
     31:   PetscCall(SNESLineSearchGetDefaultMonitor(linesearch, &monitor));

     33:   PetscCall(SNESGetObjective(snes, &objective, NULL));

     35:   /* precheck */
     36:   PetscCall(SNESLineSearchPreCheck(linesearch, X, Y, &changed_y));
     37:   lambda_old = 0.0;
     38:   if (!objective) {
     39:     fnrm_old = gnorm * gnorm;
     40:   } else {
     41:     PetscCall(SNESComputeObjective(snes, X, &fnrm_old));
     42:     SNESLineSearchCheckObjectiveDomainError(snes, fnrm_old);
     43:   }
     44:   lambda_mid = 0.5 * (lambda + lambda_old);

     46:   for (i = 0; i < max_it; i++) {
     47:     /* repeatedly cut lambda until the norm or objective function is not infinity or NaN or lambda is too small */
     48:     while (PETSC_TRUE) {
     49:       PetscCall(VecWAXPY(W, -lambda_mid, Y, X));
     50:       if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
     51:       if (!objective) {
     52:         /* compute the norm at the midpoint */
     53:         PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
     54:         if (linesearch->ops->vinorm) {
     55:           fnrm_mid = gnorm;
     56:           PetscCall((*linesearch->ops->vinorm)(snes, F, W, &fnrm_mid));
     57:         } else {
     58:           PetscCall(VecNorm(F, NORM_2, &fnrm_mid));
     59:         }

     61:         /* compute the norm at the new endpoint */
     62:         PetscCall(VecWAXPY(W, -lambda, Y, X));
     63:         if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
     64:         PetscCall((*linesearch->ops->snesfunc)(snes, W, F));
     65:         if (linesearch->ops->vinorm) {
     66:           fnrm = gnorm;
     67:           PetscCall((*linesearch->ops->vinorm)(snes, F, W, &fnrm));
     68:         } else {
     69:           PetscCall(VecNorm(F, NORM_2, &fnrm));
     70:         }
     71:         fnrm_mid = fnrm_mid * fnrm_mid;
     72:         fnrm     = fnrm * fnrm;
     73:       } else {
     74:         /* compute the objective at the midpoint */
     75:         PetscCall(SNESComputeObjective(snes, W, &fnrm_mid));

     77:         /* compute the objective at the new endpoint */
     78:         PetscCall(VecWAXPY(W, -lambda, Y, X));
     79:         PetscCall(SNESComputeObjective(snes, W, &fnrm));
     80:       }

     82:       /* if new endpoint is viable, exit */
     83:       if (!PetscIsInfOrNanReal(fnrm)) break;
     84:       if (monitor) {
     85:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
     86:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: objective function at lambda = %g is infinity or NaN, cutting lambda\n", (double)lambda));
     87:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
     88:       }

     90:       /* if smallest allowable lambda gives NaN or Inf, exit line search */
     91:       if (lambda <= minlambda) {
     92:         if (monitor) {
     93:           PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
     94:           PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: objective function at lambda = %g <= lambda_min = %g is infinity or NaN, can not further cut lambda\n", (double)lambda, (double)lambda));
     95:           PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
     96:         }
     97:         PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_REDUCT));
     98:         PetscFunctionReturn(PETSC_SUCCESS);
     99:       }

    101:       /* forbid the search from ever going back to the "failed" length that generates infinity or NaN */
    102:       maxlambda = .95 * lambda;

    104:       /* shrink lambda towards the previous one which was viable */
    105:       lambda     = .5 * (lambda + lambda_old);
    106:       lambda_mid = .5 * (lambda + lambda_old);
    107:     }

    109:     /* monitoring output */
    110:     if (monitor) {
    111:       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    112:       if (!objective) {
    113:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: lambdas = [%g, %g, %g], fnorms = [%g, %g, %g]\n", (double)lambda, (double)lambda_mid, (double)lambda_old, (double)PetscSqrtReal(fnrm), (double)PetscSqrtReal(fnrm_mid), (double)PetscSqrtReal(fnrm_old)));
    114:       } else {
    115:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: lambdas = [%g, %g, %g], obj = [%g, %g, %g]\n", (double)lambda, (double)lambda_mid, (double)lambda_old, (double)fnrm, (double)fnrm_mid, (double)fnrm_old));
    116:       }
    117:       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    118:     }

    120:     /* determine change of lambda */
    121:     delLambda = lambda - lambda_old;

    123:     /* check change of lambda tolerance */
    124:     if (PetscAbsReal(delLambda) < ltol) {
    125:       if (monitor) {
    126:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    127:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: abs(delLambda) = %g < ltol = %g\n", (double)PetscAbsReal(delLambda), (double)ltol));
    128:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    129:       }
    130:       break;
    131:     }

    133:     /* compute f'() at the end points using second order one sided differencing */
    134:     delFnrm     = (3. * fnrm - 4. * fnrm_mid + 1. * fnrm_old) / delLambda;
    135:     delFnrm_old = (-3. * fnrm_old + 4. * fnrm_mid - 1. * fnrm) / delLambda;

    137:     /* compute f''() at the midpoint using centered differencing */
    138:     del2Fnrm = (delFnrm - delFnrm_old) / delLambda;

    140:     /* check absolute tolerance */
    141:     if (PetscAbsReal(delFnrm) <= atol) {
    142:       if (monitor) {
    143:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    144:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: abs(delFnrm) = %g <= atol = %g\n", (double)PetscAbsReal(delFnrm), (double)atol));
    145:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    146:       }
    147:       break;
    148:     }

    150:     /* compute the secant (Newton) update -- always go downhill */
    151:     if (del2Fnrm > 0.) lambda_update = lambda - delFnrm / del2Fnrm;
    152:     else if (del2Fnrm < 0.) lambda_update = lambda + delFnrm / del2Fnrm;
    153:     else {
    154:       if (monitor) {
    155:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    156:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: del2Fnrm = 0\n"));
    157:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    158:       }
    159:       break;
    160:     }

    162:     /* prevent secant method from stepping out of bounds */
    163:     if (lambda_update < minlambda) lambda_update = 0.5 * (lambda + lambda_old);
    164:     if (lambda_update > maxlambda) {
    165:       lambda_update = maxlambda;
    166:       break;
    167:     }

    169:     /* if lambda is NaN or Inf, do not accept update but exit with previous lambda */
    170:     if (PetscIsInfOrNanReal(lambda_update)) {
    171:       if (monitor) {
    172:         PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    173:         PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: lambda_update is NaN or Inf\n"));
    174:         PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    175:       }
    176:       PetscCall(SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_FAILED_NANORINF));
    177:       break;
    178:     }

    180:     /* update the endpoints and the midpoint of the bracketed secant region */
    181:     lambda_old = lambda;
    182:     lambda     = lambda_update;
    183:     fnrm_old   = fnrm;
    184:     lambda_mid = 0.5 * (lambda + lambda_old);

    186:     if ((i == max_it - 1) && monitor) {
    187:       PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    188:       PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search: reached maximum number of iterations!\n"));
    189:       PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    190:     }
    191:   }

    193:   /* construct the solution */
    194:   PetscCall(VecWAXPY(W, -lambda, Y, X));
    195:   if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));

    197:   /* postcheck */
    198:   PetscCall(SNESLineSearchSetLambda(linesearch, lambda));
    199:   PetscCall(SNESLineSearchPostCheck(linesearch, X, Y, W, &changed_y, &changed_w));
    200:   if (changed_y) {
    201:     if (!changed_w) PetscCall(VecWAXPY(W, -lambda, Y, X));
    202:     if (linesearch->ops->viproject) PetscCall((*linesearch->ops->viproject)(snes, W));
    203:   }
    204:   PetscCall(VecCopy(W, X));
    205:   PetscCall((*linesearch->ops->snesfunc)(snes, X, F));

    207:   PetscCall(SNESLineSearchComputeNorms(linesearch));
    208:   PetscCall(SNESLineSearchGetNorms(linesearch, NULL, &gnorm, NULL));
    209:   SNESLineSearchCheckFunctionDomainError(snes, linesearch, gnorm);
    210:   if (monitor) {
    211:     PetscCall(PetscViewerASCIIAddTab(monitor, ((PetscObject)linesearch)->tablevel));
    212:     PetscCall(PetscViewerASCIIPrintf(monitor, "    Line search terminated: lambda = %g, fnorm = %g\n", (double)lambda, (double)gnorm));
    213:     PetscCall(PetscViewerASCIISubtractTab(monitor, ((PetscObject)linesearch)->tablevel));
    214:   }
    215:   PetscFunctionReturn(PETSC_SUCCESS);
    216: }

    218: /*MC
    219:    SNESLINESEARCHSECANT - Secant search in the L2 norm of the function or the objective function.

    221:    Attempts to solve $ \min_{\lambda} f(x_k + \lambda Y_k) $ using the secant method with the initial bracketing of $ \lambda $ between [0,damping].
    222:    $f(x_k + \lambda Y_k)$ is either the squared L2-norm of the function $||F(x_k + \lambda Y_k)||_2^2$,
    223:    or the objective function $G(x_k + \lambda Y_k)$ if it is provided with `SNESSetObjective()`
    224:    Differences of $f()$ are used to approximate the first and second derivative of $f()$ with respect to
    225:    $\lambda$, $f'()$ and $f''()$.

    227:    When an objective function is provided $f(w)$ is the objective function otherwise $f(w) = ||F(w)||^2$.
    228:    $x$ is the current step and $y$ is the search direction.

    230:    Options Database Keys:
    231: +  -snes_linesearch_max_it 1         - maximum number of iterations within the line search
    232: .  -snes_linesearch_damping 1.0      - initial `lambda` on entry to the line search
    233: .  -snes_linesearch_minlambda 1e\-12 - minimum allowable `lambda`
    234: .  -snes_linesearch_maxlambda 1.0    - maximum `lambda` (scaling of solution update) allowed
    235: .  -snes_linesearch_atol 1e\-15      - absolute tolerance for the secant method $ f'() < atol $
    236: -  -snes_linesearch_ltol 1e\-8       - minimum absolute change in `lambda` allowed

    238:    Level: advanced

    240: .seealso: [](ch_snes), `SNESLINESEARCHBT`, `SNESLINESEARCHCP`, `SNESLineSearch`, `SNESLineSearchType`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
    241: M*/
    242: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_Secant(SNESLineSearch linesearch)
    243: {
    244:   PetscFunctionBegin;
    245:   linesearch->ops->apply          = SNESLineSearchApply_Secant;
    246:   linesearch->ops->destroy        = NULL;
    247:   linesearch->ops->setfromoptions = NULL;
    248:   linesearch->ops->reset          = NULL;
    249:   linesearch->ops->view           = NULL;
    250:   linesearch->ops->setup          = NULL;

    252:   linesearch->max_it = 1;
    253:   PetscFunctionReturn(PETSC_SUCCESS);
    254: }