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