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