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