Actual source code: lmvm.c

  1: /*$Id$*/

  3: #include "lmvm.h"

  5: #define LMM_BFGS                0
  6: #define LMM_SCALED_GRADIENT     1
  7: #define LMM_GRADIENT            2

 11: static int TaoSolve_LMVM(TAO_SOLVER tao, void *solver)
 12: {
 13:   TAO_LMVM *lm = (TAO_LMVM *)solver;
 14:   TaoVec *X, *G = lm->G, *D = lm->D, *W = lm->W;
 15:   TaoVec *Xold = lm->Xold, *Gold = lm->Gold;
 16:   TaoLMVMMat *M = lm->M;

 18:   TaoTerminateReason reason;
 19:   TaoTruth success;

 21:   double f, f_full, fold, gdx, gnorm;
 22:   double step = 1.0;

 24:   double delta;

 26:   int info;
 27:   TaoInt stepType;
 28:   TaoInt iter = 0, status = 0;
 29:   TaoInt bfgsUpdates = 0;

 31:   TaoFunctionBegin;

 33:   // Get vectors we will need
 34:   info = TaoGetSolution(tao, &X); CHKERRQ(info);

 36:   // Check convergence criteria
 37:   info = TaoComputeFunctionGradient(tao, X, &f, G); CHKERRQ(info);
 38:   info = G->Norm2(&gnorm); CHKERRQ(info);
 39:   if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
 40:     SETERRQ(1, "User provided compute function generated Inf or NaN");
 41:   }

 43:   info = TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason); CHKERRQ(info);
 44:   if (reason != TAO_CONTINUE_ITERATING) {
 45:     TaoFunctionReturn(0);
 46:   }

 48:   // Set initial scaling for the function
 49:   if (f != 0.0) {
 50:     delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
 51:   }
 52:   else {
 53:     delta = 2.0 / (gnorm*gnorm);
 54:   }
 55:   info = M->SetDelta(delta); CHKERRQ(info);

 57:   // Set counter for gradient/reset steps
 58:   lm->bfgs = 0;
 59:   lm->sgrad = 0;
 60:   lm->grad = 0;

 62:   // Have not converged; continue with Newton method
 63:   while (reason == TAO_CONTINUE_ITERATING) {
 64:     // Compute direction
 65:     info = M->Update(X, G); CHKERRQ(info);
 66:     info = M->Solve(G, D, &success); CHKERRQ(info);
 67:     ++bfgsUpdates;

 69:     // Check for success (descent direction)
 70:     info = D->Dot(G, &gdx); CHKERRQ(info);
 71:     if ((gdx <= 0.0) || TaoInfOrNaN(gdx)) {
 72:       // Step is not descent or direction produced not a number
 73:       // We can assert bfgsUpdates > 1 in this case because
 74:       // the first solve produces the scaled gradient direction,
 75:       // which is guaranteed to be descent
 76:       //
 77:       // Use steepest descent direction (scaled)
 78:       ++lm->grad;

 80:       if (f != 0.0) {
 81:         delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
 82:       }
 83:       else {
 84:         delta = 2.0 / (gnorm*gnorm);
 85:       }
 86:       info = M->SetDelta(delta); CHKERRQ(info);
 87:       info = M->Reset(); CHKERRQ(info);
 88:       info = M->Update(X, G); CHKERRQ(info);
 89:       info = M->Solve(G, D, &success); CHKERRQ(info);

 91:       // On a reset, the direction cannot be not a number; it is a 
 92:       // scaled gradient step.  No need to check for this condition.
 93:       // info = D->Norm2(&dnorm); CHKERRQ(info);
 94:       // if (TaoInfOrNaN(dnorm)) {
 95:       //   SETERRQ(1, "Direction generated Not-a-Number");
 96:       // }

 98:       bfgsUpdates = 1;
 99:       ++lm->sgrad;
100:       stepType = LMM_SCALED_GRADIENT;
101:     }
102:     else {
103:       if (1 == bfgsUpdates) {
104:         // The first BFGS direction is always the scaled gradient
105:         ++lm->sgrad;
106:         stepType = LMM_SCALED_GRADIENT;
107:       }
108:       else {
109:         ++lm->bfgs;
110:         stepType = LMM_BFGS;
111:       }
112:     }
113:     info = D->Negate(); CHKERRQ(info);
114:     
115:     // Perform the linesearch
116:     fold = f;
117:     info = Xold->CopyFrom(X); CHKERRQ(info);
118:     info = Gold->CopyFrom(G); CHKERRQ(info);

120:     step = 1.0;
121:     info = TaoLineSearchApply(tao, X, G, D, W, &f, &f_full, &step, &status); CHKERRQ(info);

123:     while (status && stepType != LMM_GRADIENT) {
124:       // Linesearch failed
125:       // Reset factors and use scaled gradient step
126:       f = fold;
127:       info = X->CopyFrom(Xold); CHKERRQ(info);
128:       info = G->CopyFrom(Gold); CHKERRQ(info);
129:         
130:       switch(stepType) {
131:       case LMM_BFGS:
132:         // Failed to obtain acceptable iterate with BFGS step
133:         // Attempt to use the scaled gradient direction

135:         if (f != 0.0) {
136:           delta = 2.0 * TaoAbsDouble(f) / (gnorm*gnorm);
137:         }
138:         else {
139:           delta = 2.0 / (gnorm*gnorm);
140:         }
141:         info = M->SetDelta(delta); CHKERRQ(info);
142:         info = M->Reset(); CHKERRQ(info);
143:         info = M->Update(X, G); CHKERRQ(info);
144:         info = M->Solve(G, D, &success); CHKERRQ(info);

146:         // On a reset, the direction cannot be not a number; it is a 
147:         // scaled gradient step.  No need to check for this condition.
148:         // info = D->Norm2(&dnorm); CHKERRQ(info);
149:         // if (TaoInfOrNaN(dnorm)) {
150:         //   SETERRQ(1, "Direction generated Not-a-Number");
151:         // }
152:   
153:         bfgsUpdates = 1;
154:         ++lm->sgrad;
155:         stepType = LMM_SCALED_GRADIENT;
156:         break;

158:       case LMM_SCALED_GRADIENT:
159:         // The scaled gradient step did not produce a new iterate;
160:         // attemp to use the gradient direction.
161:         // Need to make sure we are not using a different diagonal scaling
162:         info = M->SetDelta(1.0); CHKERRQ(info);
163:         info = M->Reset(); CHKERRQ(info);
164:         info = M->Update(X, G); CHKERRQ(info);
165:         info = M->Solve(G, D, &success); CHKERRQ(info);

167:         bfgsUpdates = 1;
168:         ++lm->grad;
169:         stepType = LMM_GRADIENT;
170:         break;
171:       }
172:       info = D->Negate(); CHKERRQ(info);
173:         
174:       // Perform the linesearch
175:       step = 1.0;
176:       info = TaoLineSearchApply(tao, X, G, D, W, &f, &f_full, &step, &status); CHKERRQ(info);
177:     }

179:     if (status) {
180:       // Failed to find an improving point
181:       f = fold;
182:       info = X->CopyFrom(Xold); CHKERRQ(info);
183:       info = G->CopyFrom(Gold); CHKERRQ(info);
184:       step = 0.0;
185:     }

187:     // Check for termination
188:     info = G->Norm2(&gnorm); CHKERRQ(info);
189:     if (TaoInfOrNaN(f) || TaoInfOrNaN(gnorm)) {
190:       SETERRQ(1, "User provided compute function generated Inf or NaN");
191:     }
192:     info = TaoMonitor(tao, ++iter, f, gnorm, 0.0, step, &reason); CHKERRQ(info);
193:   }
194:   TaoFunctionReturn(0);
195: }

199: static int TaoSetUp_LMVM(TAO_SOLVER tao, void *solver)
200: {
201:   TAO_LMVM *lm = (TAO_LMVM *)solver;
202:   TaoVec *X;
203:   int info;

205:   TaoFunctionBegin;

207:   info = TaoGetSolution(tao, &X); CHKERRQ(info);
208:   info = X->Clone(&lm->G); CHKERRQ(info);
209:   info = X->Clone(&lm->D); CHKERRQ(info);
210:   info = X->Clone(&lm->W); CHKERRQ(info);

212:   // Create vectors we will need for linesearch
213:   info = X->Clone(&lm->Xold); CHKERRQ(info);
214:   info = X->Clone(&lm->Gold); CHKERRQ(info);

216:   info = TaoSetLagrangianGradientVector(tao, lm->G); CHKERRQ(info);
217:   info = TaoSetStepDirectionVector(tao, lm->D); CHKERRQ(info);

219:   // Create matrix for the limited memory approximation
220:   lm->M = new TaoLMVMMat(X);

222:   info = TaoCheckFG(tao); CHKERRQ(info);
223:   TaoFunctionReturn(0);
224: }

226: /* ---------------------------------------------------------- */
229: static int TaoSetDown_LMVM(TAO_SOLVER tao, void *solver)
230: {
231:   TAO_LMVM *lm = (TAO_LMVM *)solver;
232:   int info;

234:   TaoFunctionBegin;
235:   info = TaoVecDestroy(lm->G); CHKERRQ(info);
236:   info = TaoVecDestroy(lm->D); CHKERRQ(info);
237:   info = TaoVecDestroy(lm->W); CHKERRQ(info);

239:   info = TaoVecDestroy(lm->Xold); CHKERRQ(info);
240:   info = TaoVecDestroy(lm->Gold); CHKERRQ(info);

242:   info = TaoMatDestroy(lm->M); CHKERRQ(info);

244:   info = TaoSetLagrangianGradientVector(tao, 0); CHKERRQ(info);
245:   info = TaoSetStepDirectionVector(tao, 0); CHKERRQ(info);
246:   TaoFunctionReturn(0);
247: }

249: /*------------------------------------------------------------*/
252: static int TaoSetOptions_LMVM(TAO_SOLVER tao, void *solver)
253: {
254:   int info;

256:   TaoFunctionBegin;
257:   info = TaoOptionsHead("Limited-memory variable-metric method for unconstrained optimization"); CHKERRQ(info);
258:   info = TaoLineSearchSetFromOptions(tao); CHKERRQ(info);
259:   info = TaoOptionsTail(); CHKERRQ(info);
260:   TaoFunctionReturn(0);
261: }

263: /*------------------------------------------------------------*/
266: static int TaoView_LMVM(TAO_SOLVER tao, void *solver)
267: {
268:   TAO_LMVM *lm = (TAO_LMVM *)solver;
269:   int info;

271:   TaoFunctionBegin;
272:   info = TaoPrintInt(tao, "  Rejected matrix updates: %d\n", lm->M->GetRejects()); CHKERRQ(info);
273:   info = TaoPrintInt(tao, "  BFGS steps: %d\n", lm->bfgs); CHKERRQ(info);
274:   info = TaoPrintInt(tao, "  Scaled gradient steps: %d\n", lm->sgrad); CHKERRQ(info);
275:   info = TaoPrintInt(tao, "  Gradient steps: %d\n", lm->grad); CHKERRQ(info);
276:   info = TaoLineSearchView(tao); CHKERRQ(info);
277:   TaoFunctionReturn(0);
278: }

280: /* ---------------------------------------------------------- */

285: int TaoCreate_LMVM(TAO_SOLVER tao)
286: {
287:   TAO_LMVM *lm;
288:   int info;

290:   TaoFunctionBegin;

292:   info = TaoNew(TAO_LMVM, &lm); CHKERRQ(info);
293:   info = PetscLogObjectMemory(tao, sizeof(TAO_LMVM)); CHKERRQ(info);

295:   info = TaoSetTaoSolveRoutine(tao, TaoSolve_LMVM, (void *)lm); CHKERRQ(info);
296:   info = TaoSetTaoSetUpDownRoutines(tao, TaoSetUp_LMVM, TaoSetDown_LMVM); CHKERRQ(info);
297:   info = TaoSetTaoOptionsRoutine(tao, TaoSetOptions_LMVM); CHKERRQ(info);
298:   info = TaoSetTaoViewRoutine(tao, TaoView_LMVM); CHKERRQ(info);

300:   info = TaoSetMaximumIterates(tao, 2000); CHKERRQ(info);
301:   info = TaoSetMaximumFunctionEvaluations(tao, 4000); CHKERRQ(info);
302:   info = TaoSetTolerances(tao, 1e-4, 1e-4, 0, 0); CHKERRQ(info);
303:   
304:   info = TaoCreateMoreThuenteLineSearch(tao, 0, 0); CHKERRQ(info);
305:   TaoFunctionReturn(0);
306: }

309: // Todd: do not delete; they are needed for the component version
310: // of the code.

312: /* ---------------------------------------------------------- */
315: int TaoLMVMGetX0(TAO_SOLVER tao, TaoVec *x0)
316: {
317:   TAO_LMVM *lm;
318:   int info;

320:   TaoFunctionBegin;
321:   info=TaoGetSolverContext(tao, "tao_lmvm", (void **)&lm); CHKERRQ(info);
322:   if (lm && lm->M) {
323:     info=lm->M->GetX0(x0); CHKERRQ(info);
324:   }
325:   TaoFunctionReturn(0);
326: }

328: /* ---------------------------------------------------------- */
331: int TaoInitializeLMVMmatrix(TAO_SOLVER tao, TaoVec *HV)
332: {
333:   TAO_LMVM *lm;
334:   int info;
335:   
336:   TaoFunctionBegin;
337:   info = TaoGetSolverContext(tao, "tao_lmvm", (void **)&lm); CHKERRQ(info);
338:   if (lm && lm->M) {
339:     info = lm->M->InitialApproximation(HV); CHKERRQ(info);
340:   }
341:   TaoFunctionReturn(0);
342: }