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