Actual source code: ntr.c

      1: /*$Id: ntr.c 1.98 05/05/10 15:21:39-05:00 sarich@zorak.(none) $*/

      3: /* Note that this implementation will soon be changed.  It currently
      4: accesses data in qcg.h -- this is bad! */
      5: #include "tao_solver.h"

      7: #ifdef TAO_USE_PETSC

      9: #include "ntr.h"                /*I "tao_solver.h" I*/
     10: #include "petscksp.h"
     11: #include "src/ksp/ksp/kspimpl.h"
     12: #include "src/ksp/ksp/impls/qcg/qcg.h"
     13: #include "src/petsctao/linearsolver/taolinearsolver_petsc.h"

     15: static int TaoDestroy_NTR(TAO_SOLVER,void *);

     17: /*
     18:     TaoSolve_NTR - Implements Newton's Method with a trust region approach
     19:     for solving unconstrained minimization problems.

     21:     The basic algorithm is taken from MINPACK-2 (dstrn).

     23:     TaoSolve_NTR computes a local minimizer of a twice differentiable function
     24:     f  by applying a trust region variant of Newton's method.  At each stage
     25:     of the algorithm, we us the prconditioned conjugate gradient method to
     26:     determine an approximate minimizer of the quadratic equation

     28:              q(s) = g^T * s + .5 * s^T * H * s

     30:     subject to the Euclidean norm trust region constraint

     32:              || D * s || <= delta,

     34:     where delta is the trust region radius and D is a scaling matrix.
     35:     Here g is the gradient and H is the Hessian matrix.

     37:     Note:  TaoSolve_NTR MUST use the iterative solver KSPQCG; thus, we
     38:            set KSPQCG in this routine regardless of what the user may have
     39:            previously specified.
     40: */
     43: static int TaoSolve_NTR(TAO_SOLVER tao,void *ntrptr)
     44: {
     45:   TAO_NTR            *neP = (TAO_NTR *)ntrptr;
     46:   int                info;
     47:   double             xnorm, gnorm, max_val, ftrial, f, gdx,delta, zero = 0.0;
     48:   TaoTruth           newton;
     49:   /* Vec              X, G,  S, Xtrial; */
     50:   TaoVec             *xx, *gg, *ss=neP->S, *xtrial=neP->Xtrial;
     51:   TaoMat             *HH;
     52:   KSP                pksp;
     53:   KSP_QCG            *qcgP;
     54:   TaoTerminateReason reason=TAO_CONTINUE_ITERATING;
     55:   KSPConvergedReason kreason;
     56:   TaoTruth           success;
     57:   int                iter=0;
     58:   TaoLinearSolver *tls;

     60:   TaoFunctionBegin;

     62:   info=TaoGetSolution(tao,&xx);CHKERRQ(info);
     63:   info=TaoGetGradient(tao,&gg);CHKERRQ(info);
     64:   info=TaoGetHessian(tao,&HH);CHKERRQ(info);

     66:   info = TaoGetInitialTrustRegionRadius(tao,&delta);CHKERRQ(info);  /* trust region radius */

     68:   info = xx->Norm2(&xnorm);CHKERRQ(info);              /* xnorm = || xx || */
     69:   info = TaoComputeMeritFunctionGradient(tao,xx,&f,gg);CHKERRQ(info);
     70:   info = gg->Norm2(&gnorm);CHKERRQ(info);               /* gnorm = || gg || */

     72:   info = TaoGetLinearSolver(tao,&tls); CHKERRQ(info);
     73:   pksp=((TaoLinearSolverPetsc *)tls)->GetKSP();

     75:   info = KSPSetType(pksp,KSPQCG);CHKERRQ(info);
     76:   info = PetscLogInfo((tao,"TaoSolve_NTR: setting KSPType = KSPQCG\n")); CHKERRQ(info);
     77:   qcgP = (KSP_QCG *) pksp->data;

     79:   while (1) {

     81:     info = TaoMonitor(tao,iter++,f,gnorm,0.0,delta,&reason);CHKERRQ(info);
     82:     if (reason!=TAO_CONTINUE_ITERATING) break;

     84:     newton        = TAO_FALSE;
     85:     neP->success  = 0;
     86:     info = TaoComputeHessian(tao,xx,HH);CHKERRQ(info);

     88:     if (iter == 1) {                        /* Initialize delta */
     89:       if (delta <= 0) {
     90:         if (xnorm > zero) delta = neP->factor1*xnorm;
     91:         else {
     92:           info = TaoGetInitialTrustRegionRadius(tao,&delta);CHKERRQ(info);
     93:         }
     94:         info = HH->Norm1(&max_val);
     95:         if (info == PETSC_ERR_SUP) {
     96:           info = PetscLogInfo((tao,"TaoSolve_NTR: Initial delta computed without matrix norm info\n")); CHKERRQ(info);
     97:         } else {
     98:           if (TaoAbsDouble(max_val)<1.e-14)SETERRQ(PETSC_ERR_PLIB,"Hessian norm is too small");
     99:           delta = TaoMax(delta,gnorm/max_val);
    100:         }
    101:         delta = TAO_INFINITY;
    102:       } else {
    103:         info = TaoGetInitialTrustRegionRadius(tao,&delta);CHKERRQ(info);
    104:       }
    105:     }
    106:     do {
    107:       /* Minimize the quadratic to compute the step s */
    108:       qcgP->delta = delta;

    110:       info = TaoPreLinearSolve(tao,HH);CHKERRQ(info);
    111:       //      info = TaoLinearSolve(tao,HH,gg,ss,&success);CHKERRQ(info);
    112:       info = TaoMinQPTrust(tao, HH, gg, ss, delta, &success); CHKERRQ(info);

    114:       info = KSPGetConvergedReason(pksp,&kreason);CHKERRQ(info);
    115:       info = ss->Dot(gg,&gdx); CHKERRQ(info);
    116:       //      kreason=KSP_DIVERGED_ITS;

    118:       if (kreason != KSP_CONVERGED_QCG_NEG_CURVE && kreason != KSP_CONVERGED_QCG_CONSTRAINED) {
    119:         newton = TAO_TRUE;  /* truncated Newton step */
    120:       } else if (gdx<=0){
    121:         if ((int)kreason < 0) SETERRQ(1,"Failure in SLESSolve");
    122:       }
    123:       info=PetscLogInfo((tao,"TaoSolve_NTR: %d: ltsnrm=%g, delta=%g, q=%g \n",
    124:                  iter-1, qcgP->ltsnrm, delta, qcgP->quadratic));CHKERRQ(info);

    126:       info = xtrial->Waxpby(1.0,xx,-1.0,ss);CHKERRQ(info);

    128:       info = xtrial->Norm2(&xnorm);CHKERRQ(info);
    129:                                                           /* ftrial = f(Xtrial) */
    130:       info = TaoComputeFunction(tao,xtrial,&ftrial);CHKERRQ(info);

    132:       /* Compute the function reduction and the step size */
    133:       neP->actred = f - ftrial;
    134:       neP->prered = -qcgP->quadratic;

    136:       /* Adjust delta for the first Newton step */
    137:       if ((iter == 1) && (newton)) delta = TaoMin(delta,qcgP->ltsnrm);

    139:       if (neP->actred < neP->eta1 * neP->prered) {  /* Unsuccessful step */

    141:          info = PetscLogInfo((tao,"TaoSolve_NTR: Rejecting step\n")); CHKERRQ(info);

    143:          /* If iterate is Newton step, reduce delta to current step length */
    144:          if (newton) {
    145:            delta = qcgP->ltsnrm;
    146:            newton = TAO_FALSE;
    147:          }
    148:          delta /= 4.0;

    150:       } else {          /* Successful iteration; adjust trust radius */

    152:         neP->success = 1;
    153:         info = PetscLogInfo((tao,"TaoSolve_NTR: Accepting step\n")); CHKERRQ(info);
    154:         if (newton) {
    155:            delta = sqrt(qcgP->ltsnrm*delta);
    156:            if (neP->actred < neP->eta2 * neP->prered) delta /= 2.0;
    157:         } else if (neP->actred < neP->eta2 * neP->prered)
    158:            delta /= delta;
    159:         else if ((neP->actred >= neP->eta3 * neP->prered) &&
    160:            (neP->actred < neP->eta4 * neP->prered))
    161:            delta *= 2.0;
    162:         else if (neP->actred >= neP->eta4 * neP->prered)
    163:            delta *= 4.0;
    164:         else neP->sflag = 1;
    165:       }

    167:       neP->delta = delta;
    168:       neP->xnorm=xnorm;

    170:     } while (!neP->success);

    172:     /* Question:  If (!neP->success && break), then last step was rejected,
    173:        but convergence was detected.  Should this update really happen? */
    174:     info = xx->CopyFrom(xtrial);CHKERRQ(info);

    176:     /* Note:  At last iteration, the gradient evaluation is unnecessary */
    177:     /* Actually, already have the function, so just comput fjunk for now */
    178:     info = TaoComputeFunctionGradient(tao,xx,&f,gg);CHKERRQ(info);

    180:     info = gg->Norm2(&gnorm);CHKERRQ(info);
    181:     /* need to set step variable */

    183:   }

    185:   TaoFunctionReturn(0);
    186: }
    187: /*------------------------------------------------------------*/
    190: static int TaoSetUp_NTR(TAO_SOLVER tao, void *solver)
    191: {
    192:   int     n,info;
    193:   TaoVec  *xx;
    194:   TaoMat *HH;
    195:   TaoLinearSolver *ksp;
    196:   TAO_NTR *ctx = (TAO_NTR *)solver;

    198:   TaoFunctionBegin;
    199:   info = TaoGetSolution(tao,&xx);CHKERRQ(info);

    201:   info=xx->Clone(&ctx->S); CHKERRQ(info);
    202:   info=xx->Clone(&ctx->Xtrial); CHKERRQ(info);
    203:   info=xx->Clone(&ctx->G); CHKERRQ(info);

    205:   info = TaoSetLagrangianGradientVector(tao,ctx->G);CHKERRQ(info);
    206:   info = TaoSetStepDirectionVector(tao,ctx->S);CHKERRQ(info);

    208:   info = TaoGetHessian(tao,&HH); CHKERRQ(info);
    209:   info = TaoCreateLinearSolver(tao,HH,220,0); CHKERRQ(info);
    210:   info = TaoGetLinearSolver(tao,&ksp); CHKERRQ(info);
    211:   info = xx->GetDimension(&n);CHKERRQ(info);
    212:   info = ksp->SetTolerances(1e-8,1e-30,1e30,n);CHKERRQ(info);

    214:   info = TaoCheckFGH(tao);CHKERRQ(info);

    216:   TaoFunctionReturn(0);
    217: }
    218: /*------------------------------------------------------------*/
    221: static int TaoDestroy_NTR(TAO_SOLVER tao,void *solver)
    222: {
    223:   TAO_NTR *ctx = (TAO_NTR *)solver;
    224:   int  info;

    226:   TaoFunctionBegin;

    228:   info=TaoVecDestroy(ctx->G);CHKERRQ(info);
    229:   info=TaoVecDestroy(ctx->S);CHKERRQ(info);
    230:   info=TaoVecDestroy(ctx->Xtrial);CHKERRQ(info);
    231:   info = TaoDestroyLinearSolver(tao);CHKERRQ(info);

    233:   info = TaoSetLagrangianGradientVector(tao,0);CHKERRQ(info);
    234:   info = TaoSetStepDirectionVector(tao,0);CHKERRQ(info);

    236:   TaoFunctionReturn(0);
    237: }
    238: /*------------------------------------------------------------*/
    241: static int TaoSetOptions_NTR(TAO_SOLVER tao, void*solver)
    242: {
    243:   TAO_NTR    *ctx = (TAO_NTR *)solver;
    244:   int        info;

    246:   TaoFunctionBegin;

    248:   info = TaoOptionsHead("Newton trust region method for unconstrained optimization");CHKERRQ(info);
    249:   info = TaoOptionDouble("-tao_ntr_eta1","step is unsuccessful if actual reduction < eta1 * predicted reduction","",ctx->eta1,&ctx->eta1,0);CHKERRQ(info);

    251:   info = TaoOptionDouble("-tao_ntr_eta2","","",ctx->eta2,&ctx->eta2,0);CHKERRQ(info);
    252:   info = TaoOptionDouble("-tao_ntr_eta3","","",ctx->eta3,&ctx->eta3,0);CHKERRQ(info);
    253:   info = TaoOptionDouble("-tao_ntr_eta4","","",ctx->eta3,&ctx->eta3,0);CHKERRQ(info);
    254:   info = TaoOptionDouble("-tao_ntr_factor1","","",ctx->factor1,&ctx->factor1,0);CHKERRQ(info);

    256:   /* set the special KSP monitor for matrix-free application */
    257:   /*
    258:   info = TaoSetLinearSolverOptions(tao); CHKERRQ(info);
    259:   info = OptionsHasName(tao->prefix,"-tao_mf_ksp_monitor",&flg); CHKERRQ(info);
    260:   if (flg) {
    261:   SLES sles;
    262:     KSP ksp;
    263:     info = SLESGetKSP(sles,&ksp);CHKERRQ(info);
    264:     info = KSPSetMonitor(ksp,MatSNESMFKSPMonitor,TAO_NULL,0); CHKERRQ(info);
    265:   }
    266:   */

    268:   info = TaoOptionsTail();CHKERRQ(info);
    269:   TaoFunctionReturn(0);
    270: }
    271: /*------------------------------------------------------------*/

    275: static int TaoView_NTR(TAO_SOLVER tao,void*solver)
    276: {
    277:   TAO_NTR         *tr = (TAO_NTR *)solver;
    278:   int             info;

    280:   TaoFunctionBegin;

    282:   info = TaoPrintStatement(tao,"  Linear solver minimizes quadratic over Trust Region: \n");CHKERRQ(info);

    284:   info = TaoPrintDouble(tao,"  eta1=%g,",tr->eta1);CHKERRQ(info);
    285:   info = TaoPrintDouble(tao," eta1=%g,",tr->eta2);CHKERRQ(info);
    286:   info = TaoPrintDouble(tao," eta3=%g,",tr->eta3);CHKERRQ(info);
    287:   info = TaoPrintDouble(tao," eta4=%g\n",tr->eta4);CHKERRQ(info);
    288:   info = TaoPrintDouble(tao,"  factor1=%g\n",tr->factor1);CHKERRQ(info);

    290:   TaoFunctionReturn(0);
    291: }
    292: /*------------------------------------------------------------*/
    296: int TaoCreate_NTR(TAO_SOLVER tao)
    297: {
    298:   TAO_NTR         *neP;
    299:   int             info;

    301:   TaoFunctionBegin;

    303:   info = TaoNew(TAO_NTR,&neP); CHKERRQ(info);
    304:   info = PetscLogObjectMemory(tao,sizeof(TAO_NTR)); CHKERRQ(info);

    306:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_NTR,(void*)neP); CHKERRQ(info);
    307:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_NTR,TaoDestroy_NTR); CHKERRQ(info);
    308:   info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_NTR); CHKERRQ(info);
    309:   info=TaoSetTaoViewRoutine(tao,TaoView_NTR); CHKERRQ(info);

    311:   info = TaoSetMaximumIterates(tao,50); CHKERRQ(info);
    312:   info = TaoSetTolerances(tao,1e-10,1e-10,0,0); CHKERRQ(info);
    313:   info = TaoSetTrustRegionRadius(tao,1);CHKERRQ(info);

    315:   //  tao->converged        = TaoConverged_Default;

    317:   info = TaoSetTrustRegionTolerance(tao,1.0e-12);CHKERRQ(info);

    319:   neP->delta                 = 0.0;
    320:   neP->eta1                = 1.0e-4;
    321:   neP->eta2                = 0.25;
    322:   neP->eta3                = 0.50;
    323:   neP->eta4                = 0.90;
    324:   neP->factor1                = 1.0e-8;
    325:   neP->actred                = 0.0;
    326:   neP->prered                = 0.0;
    327:   neP->success                = 0;
    328:   neP->sflag                = 0;
    329:   neP->xnorm            = 1.0e+20;
    330:   /* Set default preconditioner to be Jacobi, to override KSP default. */
    331:   /* This implementation currently requires a symmetric preconditioner. */

    333:   /* info = KSPGetPC(tao->ksp,&pc);CHKERRQ(info);
    334:   info = PCSetType(pc,PCJACOBI);CHKERRQ(info); */

    336:   TaoFunctionReturn(0);
    337: }

    341: #endif