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