Actual source code: bnk.h

      1: /*
      2: Context for bounded Newton-Krylov type optimization algorithms
      3: */

      5: #pragma once
      6: #include <petsc/private/taoimpl.h>
      7: #include <../src/tao/bound/impls/bncg/bncg.h>

      9: typedef struct {
     10:   /* Function pointer for hessian evaluation
     11:      NOTE: This is necessary so that quasi-Newton-Krylov methods can "evaluate"
     12:      a quasi-Newton approximation while full Newton-Krylov methods call-back to
     13:      the application's Hessian */
     14:   PetscErrorCode (*computehessian)(Tao);
     15:   PetscErrorCode (*computestep)(Tao, PetscBool, KSPConvergedReason *, PetscInt *);

     17:   /* Embedded TAOBNCG */
     18:   Tao       bncg;
     19:   TAO_BNCG *bncg_ctx;
     20:   PetscInt  max_cg_its, tot_cg_its;
     21:   Vec       bncg_sol;

     23:   /* Allocated vectors */
     24:   Vec W, Xwork, Gwork, Xold, Gold;
     25:   Vec unprojected_gradient, unprojected_gradient_old;

     27:   /* Unallocated matrices and vectors */
     28:   Mat H_inactive, Hpre_inactive;
     29:   Vec X_inactive, G_inactive, inactive_work, active_work;
     30:   IS  inactive_idx, active_idx, active_lower, active_upper, active_fixed;

     32:   /* Scalar values for the solution and step */
     33:   PetscReal fold, f, gnorm, dnorm;

     35:   /* Parameters for active set estimation */
     36:   PetscReal as_tol;
     37:   PetscReal as_step;

     39:   /* BFGS preconditioner data */
     40:   PC  bfgs_pre;
     41:   Mat M;
     42:   Vec Diag_min, Diag_max;

     44:   /* Parameters when updating the perturbation added to the Hessian matrix
     45:      according to the following scheme:

     47:      pert = sval;

     49:      do until convergence
     50:        shift Hessian by pert
     51:        solve Newton system

     53:        if (linear solver failed or did not compute a descent direction)
     54:          use steepest descent direction and increase perturbation

     56:          if (0 == pert)
     57:            initialize perturbation
     58:            pert = min(imax, max(imin, imfac * norm(G)))
     59:          else
     60:            increase perturbation
     61:            pert = min(pmax, max(pgfac * pert, pmgfac * norm(G)))
     62:          fi
     63:        else
     64:          use linear solver direction and decrease perturbation

     66:          pert = min(psfac * pert, pmsfac * norm(G))
     67:          if (pert < pmin)
     68:            pert = 0
     69:          fi
     70:        fi

     72:        perform line search
     73:        function and gradient evaluation
     74:        check convergence
     75:      od
     76:   */
     77:   PetscReal sval; /*  Starting perturbation value, default zero */

     79:   PetscReal imin;  /*  Minimum perturbation added during initialization  */
     80:   PetscReal imax;  /*  Maximum perturbation added during initialization */
     81:   PetscReal imfac; /*  Merit function factor during initialization */

     83:   PetscReal pert;   /*  Current perturbation value */
     84:   PetscReal pmin;   /*  Minimum perturbation value */
     85:   PetscReal pmax;   /*  Maximum perturbation value */
     86:   PetscReal pgfac;  /*  Perturbation growth factor */
     87:   PetscReal psfac;  /*  Perturbation shrink factor */
     88:   PetscReal pmgfac; /*  Merit function growth factor */
     89:   PetscReal pmsfac; /*  Merit function shrink factor */

     91:   /* Parameters when updating the trust-region radius based on steplength
     92:      if   step < nu1            (very bad step)
     93:        radius = omega1 * min(norm(d), radius)
     94:      elif step < nu2            (bad step)
     95:        radius = omega2 * min(norm(d), radius)
     96:      elif step < nu3            (okay step)
     97:        radius = omega3 * radius;
     98:      elif step < nu4            (good step)
     99:        radius = max(omega4 * norm(d), radius)
    100:      else                       (very good step)
    101:        radius = max(omega5 * norm(d), radius)
    102:      fi
    103:   */
    104:   PetscReal nu1; /*  used to compute trust-region radius */
    105:   PetscReal nu2; /*  used to compute trust-region radius */
    106:   PetscReal nu3; /*  used to compute trust-region radius */
    107:   PetscReal nu4; /*  used to compute trust-region radius */

    109:   PetscReal omega1; /*  factor used for trust-region update */
    110:   PetscReal omega2; /*  factor used for trust-region update */
    111:   PetscReal omega3; /*  factor used for trust-region update */
    112:   PetscReal omega4; /*  factor used for trust-region update */
    113:   PetscReal omega5; /*  factor used for trust-region update */

    115:   /* Parameters when updating the trust-region radius based on reduction

    117:      kappa = ared / pred
    118:      if   kappa < eta1          (very bad step)
    119:        radius = alpha1 * min(norm(d), radius)
    120:      elif kappa < eta2          (bad step)
    121:        radius = alpha2 * min(norm(d), radius)
    122:      elif kappa < eta3          (okay step)
    123:        radius = alpha3 * radius;
    124:      elif kappa < eta4          (good step)
    125:        radius = max(alpha4 * norm(d), radius)
    126:      else                       (very good step)
    127:        radius = max(alpha5 * norm(d), radius)
    128:      fi
    129:   */
    130:   PetscReal eta1; /*  used to compute trust-region radius */
    131:   PetscReal eta2; /*  used to compute trust-region radius */
    132:   PetscReal eta3; /*  used to compute trust-region radius */
    133:   PetscReal eta4; /*  used to compute trust-region radius */

    135:   PetscReal alpha1; /*  factor used for trust-region update */
    136:   PetscReal alpha2; /*  factor used for trust-region update */
    137:   PetscReal alpha3; /*  factor used for trust-region update */
    138:   PetscReal alpha4; /*  factor used for trust-region update */
    139:   PetscReal alpha5; /*  factor used for trust-region update */

    141:   /* Parameters when updating the trust-region radius based on interpolation

    143:      kappa = ared / pred
    144:      if   kappa >= 1.0 - mu1    (very good step)
    145:        choose tau in [gamma3, gamma4]
    146:        radius = max(tau * norm(d), radius)
    147:      elif kappa >= 1.0 - mu2    (good step)
    148:        choose tau in [gamma2, gamma3]
    149:        if (tau >= 1.0)
    150:          radius = max(tau * norm(d), radius)
    151:        else
    152:          radius = tau * min(norm(d), radius)
    153:        fi
    154:      else                       (bad step)
    155:        choose tau in [gamma1, 1.0]
    156:        radius = tau * min(norm(d), radius)
    157:      fi
    158:   */
    159:   PetscReal mu1; /*  used for model agreement in interpolation */
    160:   PetscReal mu2; /*  used for model agreement in interpolation */

    162:   PetscReal gamma1; /*  factor used for interpolation */
    163:   PetscReal gamma2; /*  factor used for interpolation */
    164:   PetscReal gamma3; /*  factor used for interpolation */
    165:   PetscReal gamma4; /*  factor used for interpolation */

    167:   PetscReal theta; /*  factor used for interpolation */

    169:   /*  Parameters when initializing trust-region radius based on interpolation */
    170:   PetscReal mu1_i; /*  used for model agreement in interpolation */
    171:   PetscReal mu2_i; /*  used for model agreement in interpolation */

    173:   PetscReal gamma1_i; /*  factor used for interpolation */
    174:   PetscReal gamma2_i; /*  factor used for interpolation */
    175:   PetscReal gamma3_i; /*  factor used for interpolation */
    176:   PetscReal gamma4_i; /*  factor used for interpolation */

    178:   PetscReal theta_i; /*  factor used for interpolation */

    180:   /*  Other parameters */
    181:   PetscReal min_radius; /*  lower bound on initial radius value */
    182:   PetscReal max_radius; /*  upper bound on trust region radius */
    183:   PetscReal epsilon;    /*  tolerance used when computing ared/pred */
    184:   PetscReal dmin, dmax; /*  upper and lower bounds for the Hessian diagonal vector */

    186:   PetscInt newt;  /*  Newton directions attempted */
    187:   PetscInt bfgs;  /*  BFGS directions attempted */
    188:   PetscInt sgrad; /*  Scaled gradient directions attempted */
    189:   PetscInt grad;  /*  Gradient directions attempted */

    191:   PetscInt as_type;         /*  Active set estimation method */
    192:   PetscInt bfgs_scale_type; /*  Scaling matrix to used for the bfgs preconditioner */
    193:   PetscInt init_type;       /*  Trust-region initialization method */
    194:   PetscInt update_type;     /*  Trust-region update method */

    196:   /* Trackers for KSP solution type and convergence reasons */
    197:   PetscInt  ksp_atol;
    198:   PetscInt  ksp_rtol;
    199:   PetscInt  ksp_ctol;
    200:   PetscInt  ksp_negc;
    201:   PetscInt  ksp_dtol;
    202:   PetscInt  ksp_iter;
    203:   PetscInt  ksp_othr;
    204:   PetscBool resetksp;

    206:   /* Implementation specific context */
    207:   void *ctx;
    208: } TAO_BNK;

    210: #define BNK_NEWTON          0
    211: #define BNK_BFGS            1
    212: #define BNK_SCALED_GRADIENT 2
    213: #define BNK_GRADIENT        3

    215: #define BNK_INIT_CONSTANT      0
    216: #define BNK_INIT_DIRECTION     1
    217: #define BNK_INIT_INTERPOLATION 2
    218: #define BNK_INIT_TYPES         3

    220: #define BNK_UPDATE_STEP          0
    221: #define BNK_UPDATE_REDUCTION     1
    222: #define BNK_UPDATE_INTERPOLATION 2
    223: #define BNK_UPDATE_TYPES         3

    225: #define BNK_AS_NONE      0
    226: #define BNK_AS_BERTSEKAS 1
    227: #define BNK_AS_TYPES     2

    229: PETSC_INTERN PetscErrorCode TaoCreate_BNK(Tao);
    230: PETSC_INTERN PetscErrorCode TaoSetUp_BNK(Tao);
    231: PETSC_INTERN PetscErrorCode TaoSetFromOptions_BNK(Tao, PetscOptionItems);
    232: PETSC_INTERN PetscErrorCode TaoDestroy_BNK(Tao);
    233: PETSC_INTERN PetscErrorCode TaoView_BNK(Tao, PetscViewer);

    235: PETSC_INTERN PetscErrorCode TaoSolve_BNLS(Tao);
    236: PETSC_INTERN PetscErrorCode TaoSolve_BNTR(Tao);
    237: PETSC_INTERN PetscErrorCode TaoSolve_BNTL(Tao);

    239: PETSC_INTERN PetscErrorCode TaoBNKPreconBFGS(PC, Vec, Vec);
    240: PETSC_INTERN PetscErrorCode TaoBNKInitialize(Tao, PetscInt, PetscBool *);
    241: PETSC_INTERN PetscErrorCode TaoBNKEstimateActiveSet(Tao, PetscInt);
    242: PETSC_INTERN PetscErrorCode TaoBNKComputeHessian(Tao);
    243: PETSC_INTERN PetscErrorCode TaoBNKBoundStep(Tao, PetscInt, Vec);
    244: PETSC_INTERN PetscErrorCode TaoBNKTakeCGSteps(Tao, PetscBool *);
    245: PETSC_INTERN PetscErrorCode TaoBNKComputeStep(Tao, PetscBool, KSPConvergedReason *, PetscInt *);
    246: PETSC_INTERN PetscErrorCode TaoBNKRecomputePred(Tao, Vec, PetscReal *);
    247: PETSC_INTERN PetscErrorCode TaoBNKSafeguardStep(Tao, KSPConvergedReason, PetscInt *);
    248: PETSC_INTERN PetscErrorCode TaoBNKPerformLineSearch(Tao, PetscInt *, PetscReal *, TaoLineSearchConvergedReason *);
    249: PETSC_INTERN PetscErrorCode TaoBNKUpdateTrustRadius(Tao, PetscReal, PetscReal, PetscInt, PetscInt, PetscBool *);
    250: PETSC_INTERN PetscErrorCode TaoBNKAddStepCounts(Tao, PetscInt);