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: }