Actual source code: bfgs.c
1: #include <../src/ksp/ksp/utils/lmvm/symbrdn/symbrdn.h>
2: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>
4: /*
5: Limited-memory Broyden-Fletcher-Goldfarb-Shano method for approximating both
6: the forward product and inverse application of a Jacobian.
7: */
9: /*------------------------------------------------------------*/
11: /*
12: The solution method (approximate inverse Jacobian application) is adapted
13: from Algorithm 7.4 on page 178 of Nocedal and Wright "Numerical Optimization"
14: 2nd edition (https://doi.org/10.1007/978-0-387-40065-5). The initial inverse
15: Jacobian application falls back onto the gamma scaling recommended in equation
16: (7.20) if the user has not provided any estimation of the initial Jacobian or
17: its inverse.
19: work <- F
21: for i = k,k-1,k-2,...,0
22: rho[i] = 1 / (Y[i]^T S[i])
23: alpha[i] = rho[i] * (S[i]^T work)
24: Fwork <- work - (alpha[i] * Y[i])
25: end
27: dX <- J0^{-1} * work
29: for i = 0,1,2,...,k
30: beta = rho[i] * (Y[i]^T dX)
31: dX <- dX + ((alpha[i] - beta) * S[i])
32: end
33: */
34: PetscErrorCode MatSolve_LMVMBFGS(Mat B, Vec F, Vec dX)
35: {
36: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
37: Mat_SymBrdn *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
38: PetscInt i;
39: PetscReal *alpha, beta;
40: PetscScalar stf, ytx;
42: VecCheckSameSize(F, 2, dX, 3);
43: VecCheckMatCompatible(B, dX, 3, F, 2);
45: /* Copy the function into the work vector for the first loop */
46: VecCopy(F, lbfgs->work);
48: /* Start the first loop */
49: PetscMalloc1(lmvm->k+1, &alpha);
50: for (i = lmvm->k; i >= 0; --i) {
51: VecDot(lmvm->S[i], lbfgs->work, &stf);
52: alpha[i] = PetscRealPart(stf)/lbfgs->yts[i];
53: VecAXPY(lbfgs->work, -alpha[i], lmvm->Y[i]);
54: }
56: /* Invert the initial Jacobian onto the work vector (or apply scaling) */
57: MatSymBrdnApplyJ0Inv(B, lbfgs->work, dX);
59: /* Start the second loop */
60: for (i = 0; i <= lmvm->k; ++i) {
61: VecDot(lmvm->Y[i], dX, &ytx);
62: beta = PetscRealPart(ytx)/lbfgs->yts[i];
63: VecAXPY(dX, alpha[i]-beta, lmvm->S[i]);
64: }
65: PetscFree(alpha);
66: return 0;
67: }
69: /*------------------------------------------------------------*/
71: /*
72: The forward product for the approximate Jacobian is the matrix-free
73: implementation of Equation (6.19) in Nocedal and Wright "Numerical
74: Optimization" 2nd Edition, pg 140.
76: This forward product has the same structure as the inverse Jacobian
77: application in the DFP formulation, except with S and Y exchanging
78: roles.
80: Note: P[i] = (B_i)*S[i] terms are computed ahead of time whenever
81: the matrix is updated with a new (S[i], Y[i]) pair. This allows
82: repeated calls of MatMult inside KSP solvers without unnecessarily
83: recomputing P[i] terms in expensive nested-loops.
85: Z <- J0 * X
87: for i = 0,1,2,...,k
88: P[i] <- J0 * S[i]
89: for j = 0,1,2,...,(i-1)
90: gamma = (Y[j]^T S[i]) / (Y[j]^T S[j])
91: zeta = (S[j]^ P[i]) / (S[j]^T P[j])
92: P[i] <- P[i] - (zeta * P[j]) + (gamma * Y[j])
93: end
94: gamma = (Y[i]^T X) / (Y[i]^T S[i])
95: zeta = (S[i]^T Z) / (S[i]^T P[i])
96: Z <- Z - (zeta * P[i]) + (gamma * Y[i])
97: end
98: */
99: PetscErrorCode MatMult_LMVMBFGS(Mat B, Vec X, Vec Z)
100: {
101: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
102: Mat_SymBrdn *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
103: PetscInt i, j;
104: PetscScalar sjtpi, yjtsi, ytx, stz, stp;
106: VecCheckSameSize(X, 2, Z, 3);
107: VecCheckMatCompatible(B, X, 2, Z, 3);
109: if (lbfgs->needP) {
110: /* Pre-compute (P[i] = B_i * S[i]) */
111: for (i = 0; i <= lmvm->k; ++i) {
112: MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lbfgs->P[i]);
113: for (j = 0; j <= i-1; ++j) {
114: /* Compute the necessary dot products */
115: VecDotBegin(lmvm->S[j], lbfgs->P[i], &sjtpi);
116: VecDotBegin(lmvm->Y[j], lmvm->S[i], &yjtsi);
117: VecDotEnd(lmvm->S[j], lbfgs->P[i], &sjtpi);
118: VecDotEnd(lmvm->Y[j], lmvm->S[i], &yjtsi);
119: /* Compute the pure BFGS component of the forward product */
120: VecAXPBYPCZ(lbfgs->P[i], -PetscRealPart(sjtpi)/lbfgs->stp[j], PetscRealPart(yjtsi)/lbfgs->yts[j], 1.0, lbfgs->P[j], lmvm->Y[j]);
121: }
122: VecDot(lmvm->S[i], lbfgs->P[i], &stp);
123: lbfgs->stp[i] = PetscRealPart(stp);
124: }
125: lbfgs->needP = PETSC_FALSE;
126: }
128: /* Start the outer loop (i) for the recursive formula */
129: MatSymBrdnApplyJ0Fwd(B, X, Z);
130: for (i = 0; i <= lmvm->k; ++i) {
131: /* Get all the dot products we need */
132: VecDotBegin(lmvm->S[i], Z, &stz);
133: VecDotBegin(lmvm->Y[i], X, &ytx);
134: VecDotEnd(lmvm->S[i], Z, &stz);
135: VecDotEnd(lmvm->Y[i], X, &ytx);
136: /* Update Z_{i+1} = B_{i+1} * X */
137: VecAXPBYPCZ(Z, -PetscRealPart(stz)/lbfgs->stp[i], PetscRealPart(ytx)/lbfgs->yts[i], 1.0, lbfgs->P[i], lmvm->Y[i]);
138: }
139: return 0;
140: }
142: /*------------------------------------------------------------*/
144: static PetscErrorCode MatUpdate_LMVMBFGS(Mat B, Vec X, Vec F)
145: {
146: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
147: Mat_SymBrdn *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
148: Mat_LMVM *dbase;
149: Mat_DiagBrdn *dctx;
150: PetscInt old_k, i;
151: PetscReal curvtol;
152: PetscScalar curvature, ytytmp, ststmp;
154: if (!lmvm->m) return 0;
155: if (lmvm->prev_set) {
156: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
157: VecAYPX(lmvm->Xprev, -1.0, X);
158: VecAYPX(lmvm->Fprev, -1.0, F);
159: /* Test if the updates can be accepted */
160: VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
161: VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
162: VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
163: VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
164: if (PetscRealPart(ststmp) < lmvm->eps) {
165: curvtol = 0.0;
166: } else {
167: curvtol = lmvm->eps * PetscRealPart(ststmp);
168: }
169: if (PetscRealPart(curvature) > curvtol) {
170: /* Update is good, accept it */
171: lbfgs->watchdog = 0;
172: lbfgs->needP = PETSC_TRUE;
173: old_k = lmvm->k;
174: MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
175: /* If we hit the memory limit, shift the yts, yty and sts arrays */
176: if (old_k == lmvm->k) {
177: for (i = 0; i <= lmvm->k-1; ++i) {
178: lbfgs->yts[i] = lbfgs->yts[i+1];
179: lbfgs->yty[i] = lbfgs->yty[i+1];
180: lbfgs->sts[i] = lbfgs->sts[i+1];
181: }
182: }
183: /* Update history of useful scalars */
184: VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &ytytmp);
185: lbfgs->yts[lmvm->k] = PetscRealPart(curvature);
186: lbfgs->yty[lmvm->k] = PetscRealPart(ytytmp);
187: lbfgs->sts[lmvm->k] = PetscRealPart(ststmp);
188: /* Compute the scalar scale if necessary */
189: if (lbfgs->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) {
190: MatSymBrdnComputeJ0Scalar(B);
191: }
192: } else {
193: /* Update is bad, skip it */
194: ++lmvm->nrejects;
195: ++lbfgs->watchdog;
196: }
197: } else {
198: switch (lbfgs->scale_type) {
199: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
200: dbase = (Mat_LMVM*)lbfgs->D->data;
201: dctx = (Mat_DiagBrdn*)dbase->ctx;
202: VecSet(dctx->invD, lbfgs->delta);
203: break;
204: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
205: lbfgs->sigma = lbfgs->delta;
206: break;
207: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
208: lbfgs->sigma = 1.0;
209: break;
210: default:
211: break;
212: }
213: }
215: /* Update the scaling */
216: if (lbfgs->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
217: MatLMVMUpdate(lbfgs->D, X, F);
218: }
220: if (lbfgs->watchdog > lbfgs->max_seq_rejects) {
221: MatLMVMReset(B, PETSC_FALSE);
222: if (lbfgs->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
223: MatLMVMReset(lbfgs->D, PETSC_FALSE);
224: }
225: }
227: /* Save the solution and function to be used in the next update */
228: VecCopy(X, lmvm->Xprev);
229: VecCopy(F, lmvm->Fprev);
230: lmvm->prev_set = PETSC_TRUE;
231: return 0;
232: }
234: /*------------------------------------------------------------*/
236: static PetscErrorCode MatCopy_LMVMBFGS(Mat B, Mat M, MatStructure str)
237: {
238: Mat_LMVM *bdata = (Mat_LMVM*)B->data;
239: Mat_SymBrdn *bctx = (Mat_SymBrdn*)bdata->ctx;
240: Mat_LMVM *mdata = (Mat_LMVM*)M->data;
241: Mat_SymBrdn *mctx = (Mat_SymBrdn*)mdata->ctx;
242: PetscInt i;
244: mctx->needP = bctx->needP;
245: for (i=0; i<=bdata->k; ++i) {
246: mctx->stp[i] = bctx->stp[i];
247: mctx->yts[i] = bctx->yts[i];
248: VecCopy(bctx->P[i], mctx->P[i]);
249: }
250: mctx->scale_type = bctx->scale_type;
251: mctx->alpha = bctx->alpha;
252: mctx->beta = bctx->beta;
253: mctx->rho = bctx->rho;
254: mctx->delta = bctx->delta;
255: mctx->sigma_hist = bctx->sigma_hist;
256: mctx->watchdog = bctx->watchdog;
257: mctx->max_seq_rejects = bctx->max_seq_rejects;
258: switch (bctx->scale_type) {
259: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
260: mctx->sigma = bctx->sigma;
261: break;
262: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
263: MatCopy(bctx->D, mctx->D, SAME_NONZERO_PATTERN);
264: break;
265: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
266: mctx->sigma = 1.0;
267: break;
268: default:
269: break;
270: }
271: return 0;
272: }
274: /*------------------------------------------------------------*/
276: static PetscErrorCode MatReset_LMVMBFGS(Mat B, PetscBool destructive)
277: {
278: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
279: Mat_SymBrdn *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
280: Mat_LMVM *dbase;
281: Mat_DiagBrdn *dctx;
283: lbfgs->watchdog = 0;
284: lbfgs->needP = PETSC_TRUE;
285: if (lbfgs->allocated) {
286: if (destructive) {
287: VecDestroy(&lbfgs->work);
288: PetscFree4(lbfgs->stp, lbfgs->yts, lbfgs->yty, lbfgs->sts);
289: VecDestroyVecs(lmvm->m, &lbfgs->P);
290: switch (lbfgs->scale_type) {
291: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
292: MatLMVMReset(lbfgs->D, PETSC_TRUE);
293: break;
294: default:
295: break;
296: }
297: lbfgs->allocated = PETSC_FALSE;
298: } else {
299: switch (lbfgs->scale_type) {
300: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
301: lbfgs->sigma = lbfgs->delta;
302: break;
303: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
304: MatLMVMReset(lbfgs->D, PETSC_FALSE);
305: dbase = (Mat_LMVM*)lbfgs->D->data;
306: dctx = (Mat_DiagBrdn*)dbase->ctx;
307: VecSet(dctx->invD, lbfgs->delta);
308: break;
309: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
310: lbfgs->sigma = 1.0;
311: break;
312: default:
313: break;
314: }
315: }
316: }
317: MatReset_LMVM(B, destructive);
318: return 0;
319: }
321: /*------------------------------------------------------------*/
323: static PetscErrorCode MatAllocate_LMVMBFGS(Mat B, Vec X, Vec F)
324: {
325: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
326: Mat_SymBrdn *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
328: MatAllocate_LMVM(B, X, F);
329: if (!lbfgs->allocated) {
330: VecDuplicate(X, &lbfgs->work);
331: PetscMalloc4(lmvm->m, &lbfgs->stp, lmvm->m, &lbfgs->yts, lmvm->m, &lbfgs->yty, lmvm->m, &lbfgs->sts);
332: if (lmvm->m > 0) {
333: VecDuplicateVecs(X, lmvm->m, &lbfgs->P);
334: }
335: switch (lbfgs->scale_type) {
336: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
337: MatLMVMAllocate(lbfgs->D, X, F);
338: break;
339: default:
340: break;
341: }
342: lbfgs->allocated = PETSC_TRUE;
343: }
344: return 0;
345: }
347: /*------------------------------------------------------------*/
349: static PetscErrorCode MatDestroy_LMVMBFGS(Mat B)
350: {
351: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
352: Mat_SymBrdn *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
354: if (lbfgs->allocated) {
355: VecDestroy(&lbfgs->work);
356: PetscFree4(lbfgs->stp, lbfgs->yts, lbfgs->yty, lbfgs->sts);
357: VecDestroyVecs(lmvm->m, &lbfgs->P);
358: lbfgs->allocated = PETSC_FALSE;
359: }
360: MatDestroy(&lbfgs->D);
361: PetscFree(lmvm->ctx);
362: MatDestroy_LMVM(B);
363: return 0;
364: }
366: /*------------------------------------------------------------*/
368: static PetscErrorCode MatSetUp_LMVMBFGS(Mat B)
369: {
370: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
371: Mat_SymBrdn *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
372: PetscInt n, N;
374: MatSetUp_LMVM(B);
375: lbfgs->max_seq_rejects = lmvm->m/2;
376: if (!lbfgs->allocated) {
377: VecDuplicate(lmvm->Xprev, &lbfgs->work);
378: PetscMalloc4(lmvm->m, &lbfgs->stp, lmvm->m, &lbfgs->yts, lmvm->m, &lbfgs->yty, lmvm->m, &lbfgs->sts);
379: if (lmvm->m > 0) {
380: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbfgs->P);
381: }
382: switch (lbfgs->scale_type) {
383: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
384: MatGetLocalSize(B, &n, &n);
385: MatGetSize(B, &N, &N);
386: MatSetSizes(lbfgs->D, n, n, N, N);
387: MatSetUp(lbfgs->D);
388: break;
389: default:
390: break;
391: }
392: lbfgs->allocated = PETSC_TRUE;
393: }
394: return 0;
395: }
397: /*------------------------------------------------------------*/
399: static PetscErrorCode MatSetFromOptions_LMVMBFGS(PetscOptionItems *PetscOptionsObject, Mat B)
400: {
401: MatSetFromOptions_LMVM(PetscOptionsObject, B);
402: PetscOptionsHead(PetscOptionsObject,"L-BFGS method for approximating SPD Jacobian actions (MATLMVMBFGS)");
403: MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionsObject, B);
404: PetscOptionsTail();
405: return 0;
406: }
408: /*------------------------------------------------------------*/
410: PetscErrorCode MatCreate_LMVMBFGS(Mat B)
411: {
412: Mat_LMVM *lmvm;
413: Mat_SymBrdn *lbfgs;
415: MatCreate_LMVMSymBrdn(B);
416: PetscObjectChangeTypeName((PetscObject)B, MATLMVMBFGS);
417: B->ops->setup = MatSetUp_LMVMBFGS;
418: B->ops->destroy = MatDestroy_LMVMBFGS;
419: B->ops->setfromoptions = MatSetFromOptions_LMVMBFGS;
420: B->ops->solve = MatSolve_LMVMBFGS;
422: lmvm = (Mat_LMVM*)B->data;
423: lmvm->ops->allocate = MatAllocate_LMVMBFGS;
424: lmvm->ops->reset = MatReset_LMVMBFGS;
425: lmvm->ops->update = MatUpdate_LMVMBFGS;
426: lmvm->ops->mult = MatMult_LMVMBFGS;
427: lmvm->ops->copy = MatCopy_LMVMBFGS;
429: lbfgs = (Mat_SymBrdn*)lmvm->ctx;
430: lbfgs->needQ = PETSC_FALSE;
431: lbfgs->phi = 0.0;
432: return 0;
433: }
435: /*------------------------------------------------------------*/
437: /*@
438: MatCreateLMVMBFGS - Creates a limited-memory Broyden-Fletcher-Goldfarb-Shano (BFGS)
439: matrix used for approximating Jacobians. L-BFGS is symmetric positive-definite by
440: construction, and is commonly used to approximate Hessians in optimization
441: problems.
443: The provided local and global sizes must match the solution and function vectors
444: used with MatLMVMUpdate() and MatSolve(). The resulting L-BFGS matrix will have
445: storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
446: parallel. To use the L-BFGS matrix with other vector types, the matrix must be
447: created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
448: This ensures that the internal storage and work vectors are duplicated from the
449: correct type of vector.
451: Collective
453: Input Parameters:
454: + comm - MPI communicator, set to PETSC_COMM_SELF
455: . n - number of local rows for storage vectors
456: - N - global size of the storage vectors
458: Output Parameter:
459: . B - the matrix
461: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
462: paradigm instead of this routine directly.
464: Options Database Keys:
465: + -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
466: . -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
467: . -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
468: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
469: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
470: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
471: - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling
473: Level: intermediate
475: .seealso: MatCreate(), MATLMVM, MATLMVMBFGS, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
476: MatCreateLMVMBrdn(), MatCreateLMVMBadBrdn(), MatCreateLMVMSymBrdn()
477: @*/
478: PetscErrorCode MatCreateLMVMBFGS(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
479: {
480: MatCreate(comm, B);
481: MatSetSizes(*B, n, n, N, N);
482: MatSetType(*B, MATLMVMBFGS);
483: MatSetUp(*B);
484: return 0;
485: }