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