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: PetscFunctionBegin;
43: VecCheckSameSize(F, 2, dX, 3);
44: VecCheckMatCompatible(B, dX, 3, F, 2);
46: /* Copy the function into the work vector for the first loop */
47: PetscCall(VecCopy(F, lbfgs->work));
49: /* Start the first loop */
50: PetscCall(PetscMalloc1(lmvm->k + 1, &alpha));
51: for (i = lmvm->k; i >= 0; --i) {
52: PetscCall(VecDot(lmvm->S[i], lbfgs->work, &stf));
53: alpha[i] = PetscRealPart(stf) / lbfgs->yts[i];
54: PetscCall(VecAXPY(lbfgs->work, -alpha[i], lmvm->Y[i]));
55: }
57: /* Invert the initial Jacobian onto the work vector (or apply scaling) */
58: PetscCall(MatSymBrdnApplyJ0Inv(B, lbfgs->work, dX));
60: /* Start the second loop */
61: for (i = 0; i <= lmvm->k; ++i) {
62: PetscCall(VecDot(lmvm->Y[i], dX, &ytx));
63: beta = PetscRealPart(ytx) / lbfgs->yts[i];
64: PetscCall(VecAXPY(dX, alpha[i] - beta, lmvm->S[i]));
65: }
66: PetscCall(PetscFree(alpha));
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: /*------------------------------------------------------------*/
72: /*
73: The forward product for the approximate Jacobian is the matrix-free
74: implementation of Equation (6.19) in Nocedal and Wright "Numerical
75: Optimization" 2nd Edition, pg 140.
77: This forward product has the same structure as the inverse Jacobian
78: application in the DFP formulation, except with S and Y exchanging
79: roles.
81: Note: P[i] = (B_i)*S[i] terms are computed ahead of time whenever
82: the matrix is updated with a new (S[i], Y[i]) pair. This allows
83: repeated calls of MatMult inside KSP solvers without unnecessarily
84: recomputing P[i] terms in expensive nested-loops.
86: Z <- J0 * X
88: for i = 0,1,2,...,k
89: P[i] <- J0 * S[i]
90: for j = 0,1,2,...,(i-1)
91: gamma = (Y[j]^T S[i]) / (Y[j]^T S[j])
92: zeta = (S[j]^ P[i]) / (S[j]^T P[j])
93: P[i] <- P[i] - (zeta * P[j]) + (gamma * Y[j])
94: end
95: gamma = (Y[i]^T X) / (Y[i]^T S[i])
96: zeta = (S[i]^T Z) / (S[i]^T P[i])
97: Z <- Z - (zeta * P[i]) + (gamma * Y[i])
98: end
99: */
100: PetscErrorCode MatMult_LMVMBFGS(Mat B, Vec X, Vec Z)
101: {
102: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
103: Mat_SymBrdn *lbfgs = (Mat_SymBrdn *)lmvm->ctx;
104: PetscInt i, j;
105: PetscScalar sjtpi, yjtsi, ytx, stz, stp;
107: PetscFunctionBegin;
108: VecCheckSameSize(X, 2, Z, 3);
109: VecCheckMatCompatible(B, X, 2, Z, 3);
111: if (lbfgs->needP) {
112: /* Pre-compute (P[i] = B_i * S[i]) */
113: for (i = 0; i <= lmvm->k; ++i) {
114: PetscCall(MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lbfgs->P[i]));
115: /* Compute the necessary dot products */
116: PetscCall(VecMDot(lmvm->S[i], i, lmvm->Y, lbfgs->workscalar));
117: for (j = 0; j < i; ++j) {
118: yjtsi = lbfgs->workscalar[j];
119: PetscCall(VecDot(lmvm->S[j], lbfgs->P[i], &sjtpi));
120: /* Compute the pure BFGS component of the forward product */
121: PetscCall(VecAXPBYPCZ(lbfgs->P[i], -PetscRealPart(sjtpi) / lbfgs->stp[j], PetscRealPart(yjtsi) / lbfgs->yts[j], 1.0, lbfgs->P[j], lmvm->Y[j]));
122: }
123: PetscCall(VecDot(lmvm->S[i], lbfgs->P[i], &stp));
124: lbfgs->stp[i] = PetscRealPart(stp);
125: }
126: lbfgs->needP = PETSC_FALSE;
127: }
129: /* Start the outer loop (i) for the recursive formula */
130: PetscCall(MatSymBrdnApplyJ0Fwd(B, X, Z));
131: /* Get all the dot products we need */
132: PetscCall(VecMDot(X, lmvm->k + 1, lmvm->Y, lbfgs->workscalar));
133: for (i = 0; i <= lmvm->k; ++i) {
134: ytx = lbfgs->workscalar[i];
135: PetscCall(VecDot(lmvm->S[i], Z, &stz));
136: /* Update Z_{i+1} = B_{i+1} * X */
137: PetscCall(VecAXPBYPCZ(Z, -PetscRealPart(stz) / lbfgs->stp[i], PetscRealPart(ytx) / lbfgs->yts[i], 1.0, lbfgs->P[i], lmvm->Y[i]));
138: }
139: PetscFunctionReturn(PETSC_SUCCESS);
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, ststmp;
152: PetscScalar curvature, ytytmp;
154: PetscFunctionBegin;
155: if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
156: if (lmvm->prev_set) {
157: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
158: PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
159: PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));
161: /* Test if the updates can be accepted */
162: PetscCall(VecDotNorm2(lmvm->Xprev, lmvm->Fprev, &curvature, &ststmp));
163: if (ststmp < lmvm->eps) curvtol = 0.0;
164: else curvtol = lmvm->eps * ststmp;
166: if (PetscRealPart(curvature) > curvtol) {
167: /* Update is good, accept it */
168: lbfgs->watchdog = 0;
169: lbfgs->needP = PETSC_TRUE;
170: old_k = lmvm->k;
171: PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
172: /* If we hit the memory limit, shift the yts, yty and sts arrays */
173: if (old_k == lmvm->k) {
174: for (i = 0; i <= lmvm->k - 1; ++i) {
175: lbfgs->yts[i] = lbfgs->yts[i + 1];
176: lbfgs->yty[i] = lbfgs->yty[i + 1];
177: lbfgs->sts[i] = lbfgs->sts[i + 1];
178: }
179: }
180: /* Update history of useful scalars */
181: PetscCall(VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &ytytmp));
182: lbfgs->yts[lmvm->k] = PetscRealPart(curvature);
183: lbfgs->yty[lmvm->k] = PetscRealPart(ytytmp);
184: lbfgs->sts[lmvm->k] = ststmp;
185: /* Compute the scalar scale if necessary */
186: if (lbfgs->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) PetscCall(MatSymBrdnComputeJ0Scalar(B));
187: } else {
188: /* Update is bad, skip it */
189: ++lmvm->nrejects;
190: ++lbfgs->watchdog;
191: }
192: } else {
193: switch (lbfgs->scale_type) {
194: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
195: dbase = (Mat_LMVM *)lbfgs->D->data;
196: dctx = (Mat_DiagBrdn *)dbase->ctx;
197: PetscCall(VecSet(dctx->invD, lbfgs->delta));
198: break;
199: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
200: lbfgs->sigma = lbfgs->delta;
201: break;
202: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
203: lbfgs->sigma = 1.0;
204: break;
205: default:
206: break;
207: }
208: }
210: /* Update the scaling */
211: if (lbfgs->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatLMVMUpdate(lbfgs->D, X, F));
213: if (lbfgs->watchdog > lbfgs->max_seq_rejects) {
214: PetscCall(MatLMVMReset(B, PETSC_FALSE));
215: if (lbfgs->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatLMVMReset(lbfgs->D, PETSC_FALSE));
216: }
218: /* Save the solution and function to be used in the next update */
219: PetscCall(VecCopy(X, lmvm->Xprev));
220: PetscCall(VecCopy(F, lmvm->Fprev));
221: lmvm->prev_set = PETSC_TRUE;
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
225: /*------------------------------------------------------------*/
227: static PetscErrorCode MatCopy_LMVMBFGS(Mat B, Mat M, MatStructure str)
228: {
229: Mat_LMVM *bdata = (Mat_LMVM *)B->data;
230: Mat_SymBrdn *bctx = (Mat_SymBrdn *)bdata->ctx;
231: Mat_LMVM *mdata = (Mat_LMVM *)M->data;
232: Mat_SymBrdn *mctx = (Mat_SymBrdn *)mdata->ctx;
233: PetscInt i;
235: PetscFunctionBegin;
236: mctx->needP = bctx->needP;
237: for (i = 0; i <= bdata->k; ++i) {
238: mctx->stp[i] = bctx->stp[i];
239: mctx->yts[i] = bctx->yts[i];
240: PetscCall(VecCopy(bctx->P[i], mctx->P[i]));
241: }
242: mctx->scale_type = bctx->scale_type;
243: mctx->alpha = bctx->alpha;
244: mctx->beta = bctx->beta;
245: mctx->rho = bctx->rho;
246: mctx->delta = bctx->delta;
247: mctx->sigma_hist = bctx->sigma_hist;
248: mctx->watchdog = bctx->watchdog;
249: mctx->max_seq_rejects = bctx->max_seq_rejects;
250: switch (bctx->scale_type) {
251: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
252: mctx->sigma = bctx->sigma;
253: break;
254: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
255: PetscCall(MatCopy(bctx->D, mctx->D, SAME_NONZERO_PATTERN));
256: break;
257: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
258: mctx->sigma = 1.0;
259: break;
260: default:
261: break;
262: }
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: /*------------------------------------------------------------*/
268: static PetscErrorCode MatReset_LMVMBFGS(Mat B, PetscBool destructive)
269: {
270: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
271: Mat_SymBrdn *lbfgs = (Mat_SymBrdn *)lmvm->ctx;
272: Mat_LMVM *dbase;
273: Mat_DiagBrdn *dctx;
275: PetscFunctionBegin;
276: lbfgs->watchdog = 0;
277: lbfgs->needP = PETSC_TRUE;
278: if (lbfgs->allocated) {
279: if (destructive) {
280: PetscCall(VecDestroy(&lbfgs->work));
281: PetscCall(PetscFree5(lbfgs->stp, lbfgs->yts, lbfgs->yty, lbfgs->sts, lbfgs->workscalar));
282: PetscCall(VecDestroyVecs(lmvm->m, &lbfgs->P));
283: switch (lbfgs->scale_type) {
284: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
285: PetscCall(MatLMVMReset(lbfgs->D, PETSC_TRUE));
286: break;
287: default:
288: break;
289: }
290: lbfgs->allocated = PETSC_FALSE;
291: } else {
292: switch (lbfgs->scale_type) {
293: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
294: lbfgs->sigma = lbfgs->delta;
295: break;
296: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
297: PetscCall(MatLMVMReset(lbfgs->D, PETSC_FALSE));
298: dbase = (Mat_LMVM *)lbfgs->D->data;
299: dctx = (Mat_DiagBrdn *)dbase->ctx;
300: PetscCall(VecSet(dctx->invD, lbfgs->delta));
301: break;
302: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
303: lbfgs->sigma = 1.0;
304: break;
305: default:
306: break;
307: }
308: }
309: }
310: PetscCall(MatReset_LMVM(B, destructive));
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: /*------------------------------------------------------------*/
316: static PetscErrorCode MatAllocate_LMVMBFGS(Mat B, Vec X, Vec F)
317: {
318: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
319: Mat_SymBrdn *lbfgs = (Mat_SymBrdn *)lmvm->ctx;
321: PetscFunctionBegin;
322: PetscCall(MatAllocate_LMVM(B, X, F));
323: if (!lbfgs->allocated) {
324: PetscCall(VecDuplicate(X, &lbfgs->work));
325: PetscCall(PetscMalloc5(lmvm->m, &lbfgs->stp, lmvm->m, &lbfgs->yts, lmvm->m, &lbfgs->yty, lmvm->m, &lbfgs->sts, lmvm->m, &lbfgs->workscalar));
326: if (lmvm->m > 0) PetscCall(VecDuplicateVecs(X, lmvm->m, &lbfgs->P));
327: switch (lbfgs->scale_type) {
328: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
329: PetscCall(MatLMVMAllocate(lbfgs->D, X, F));
330: break;
331: default:
332: break;
333: }
334: lbfgs->allocated = PETSC_TRUE;
335: }
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: /*------------------------------------------------------------*/
341: static PetscErrorCode MatDestroy_LMVMBFGS(Mat B)
342: {
343: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
344: Mat_SymBrdn *lbfgs = (Mat_SymBrdn *)lmvm->ctx;
346: PetscFunctionBegin;
347: if (lbfgs->allocated) {
348: PetscCall(VecDestroy(&lbfgs->work));
349: PetscCall(PetscFree5(lbfgs->stp, lbfgs->yts, lbfgs->yty, lbfgs->sts, lbfgs->workscalar));
350: PetscCall(VecDestroyVecs(lmvm->m, &lbfgs->P));
351: lbfgs->allocated = PETSC_FALSE;
352: }
353: PetscCall(MatDestroy(&lbfgs->D));
354: PetscCall(PetscFree(lmvm->ctx));
355: PetscCall(MatDestroy_LMVM(B));
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: /*------------------------------------------------------------*/
361: static PetscErrorCode MatSetUp_LMVMBFGS(Mat B)
362: {
363: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
364: Mat_SymBrdn *lbfgs = (Mat_SymBrdn *)lmvm->ctx;
365: PetscInt n, N;
367: PetscFunctionBegin;
368: PetscCall(MatSetUp_LMVM(B));
369: lbfgs->max_seq_rejects = lmvm->m / 2;
370: if (!lbfgs->allocated) {
371: PetscCall(VecDuplicate(lmvm->Xprev, &lbfgs->work));
372: PetscCall(PetscMalloc5(lmvm->m, &lbfgs->stp, lmvm->m, &lbfgs->yts, lmvm->m, &lbfgs->yty, lmvm->m, &lbfgs->sts, lmvm->m, &lbfgs->workscalar));
373: if (lmvm->m > 0) PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbfgs->P));
374: switch (lbfgs->scale_type) {
375: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
376: PetscCall(MatGetLocalSize(B, &n, &n));
377: PetscCall(MatGetSize(B, &N, &N));
378: PetscCall(MatSetSizes(lbfgs->D, n, n, N, N));
379: PetscCall(MatSetUp(lbfgs->D));
380: break;
381: default:
382: break;
383: }
384: lbfgs->allocated = PETSC_TRUE;
385: }
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: /*------------------------------------------------------------*/
391: static PetscErrorCode MatSetFromOptions_LMVMBFGS(Mat B, PetscOptionItems *PetscOptionsObject)
392: {
393: PetscFunctionBegin;
394: PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
395: PetscOptionsHeadBegin(PetscOptionsObject, "L-BFGS method for approximating SPD Jacobian actions (MATLMVMBFGS)");
396: PetscCall(MatSetFromOptions_LMVMSymBrdn_Private(B, PetscOptionsObject));
397: PetscOptionsHeadEnd();
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
401: /*------------------------------------------------------------*/
403: PetscErrorCode MatCreate_LMVMBFGS(Mat B)
404: {
405: Mat_LMVM *lmvm;
406: Mat_SymBrdn *lbfgs;
408: PetscFunctionBegin;
409: PetscCall(MatCreate_LMVMSymBrdn(B));
410: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMBFGS));
411: B->ops->setup = MatSetUp_LMVMBFGS;
412: B->ops->destroy = MatDestroy_LMVMBFGS;
413: B->ops->setfromoptions = MatSetFromOptions_LMVMBFGS;
414: B->ops->solve = MatSolve_LMVMBFGS;
416: lmvm = (Mat_LMVM *)B->data;
417: lmvm->ops->allocate = MatAllocate_LMVMBFGS;
418: lmvm->ops->reset = MatReset_LMVMBFGS;
419: lmvm->ops->update = MatUpdate_LMVMBFGS;
420: lmvm->ops->mult = MatMult_LMVMBFGS;
421: lmvm->ops->copy = MatCopy_LMVMBFGS;
423: lbfgs = (Mat_SymBrdn *)lmvm->ctx;
424: lbfgs->needQ = PETSC_FALSE;
425: lbfgs->phi = 0.0;
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: /*------------------------------------------------------------*/
431: /*@
432: MatCreateLMVMBFGS - Creates a limited-memory Broyden-Fletcher-Goldfarb-Shano (BFGS)
433: matrix used for approximating Jacobians. L-BFGS is symmetric positive-definite by
434: construction, and is commonly used to approximate Hessians in optimization
435: problems.
437: To use the L-BFGS matrix with other vector types, the matrix must be
438: created using `MatCreate()` and `MatSetType()`, followed by `MatLMVMAllocate()`.
439: This ensures that the internal storage and work vectors are duplicated from the
440: correct type of vector.
442: Collective
444: Input Parameters:
445: + comm - MPI communicator
446: . n - number of local rows for storage vectors
447: - N - global size of the storage vectors
449: Output Parameter:
450: . B - the matrix
452: Options Database Keys:
453: + -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
454: . -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
455: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
456: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
457: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
458: - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling
460: Level: intermediate
462: Note:
463: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
464: paradigm instead of this routine directly.
466: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMBFGS`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
467: `MatCreateLMVMBrdn()`, `MatCreateLMVMBadBrdn()`, `MatCreateLMVMSymBrdn()`
468: @*/
469: PetscErrorCode MatCreateLMVMBFGS(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
470: {
471: PetscFunctionBegin;
472: PetscCall(KSPInitializePackage());
473: PetscCall(MatCreate(comm, B));
474: PetscCall(MatSetSizes(*B, n, n, N, N));
475: PetscCall(MatSetType(*B, MATLMVMBFGS));
476: PetscCall(MatSetUp(*B));
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }