Actual source code: dfp.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 Davidon-Fletcher-Powell 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
 13:   matrix-vector product version of the recursive formula given in
 14:   Equation (6.15) of Nocedal and Wright "Numerical Optimization" 2nd
 15:   edition, pg 139.

 17:   Note: Q[i] = (B_i)^{-1}*S[i] terms are computed ahead of time whenever
 18:   the matrix is updated with a new (S[i], Y[i]) pair. This allows
 19:   repeated calls of MatSolve without incurring redundant computation.

 21:   dX <- J0^{-1} * F

 23:   for i = 0,1,2,...,k
 24:     # Q[i] = (B_i)^{-1} * Y[i]
 25:     gamma = (S[i]^T F) / (Y[i]^T S[i])
 26:     zeta = (Y[i]^T dX) / (Y[i]^T Q[i])
 27:     dX <- dX + (gamma * S[i]) - (zeta * Q[i])
 28:   end
 29: */
 30: PetscErrorCode MatSolve_LMVMDFP(Mat B, Vec F, Vec dX)
 31: {
 32:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
 33:   Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;
 34:   PetscInt     i, j;
 35:   PetscScalar  yjtqi, sjtyi, ytx, stf, ytq;

 37:   PetscFunctionBegin;
 40:   VecCheckSameSize(F, 2, dX, 3);
 41:   VecCheckMatCompatible(B, dX, 3, F, 2);

 43:   if (ldfp->needQ) {
 44:     /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
 45:     for (i = 0; i <= lmvm->k; ++i) {
 46:       PetscCall(MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], ldfp->Q[i]));
 47:       /* Compute the necessary dot products */
 48:       PetscCall(VecMDot(lmvm->Y[i], i, lmvm->S, ldfp->workscalar));
 49:       for (j = 0; j < i; ++j) {
 50:         sjtyi = ldfp->workscalar[j];
 51:         PetscCall(VecDot(lmvm->Y[j], ldfp->Q[i], &yjtqi));
 52:         /* Compute the pure DFP component of the inverse application*/
 53:         PetscCall(VecAXPBYPCZ(ldfp->Q[i], -PetscRealPart(yjtqi) / ldfp->ytq[j], PetscRealPart(sjtyi) / ldfp->yts[j], 1.0, ldfp->Q[j], lmvm->S[j]));
 54:       }
 55:       PetscCall(VecDot(lmvm->Y[i], ldfp->Q[i], &ytq));
 56:       ldfp->ytq[i] = PetscRealPart(ytq);
 57:     }
 58:     ldfp->needQ = PETSC_FALSE;
 59:   }

 61:   /* Start the outer loop (i) for the recursive formula */
 62:   PetscCall(MatSymBrdnApplyJ0Inv(B, F, dX));
 63:   /* Get all the dot products we need */
 64:   PetscCall(VecMDot(F, lmvm->k + 1, lmvm->S, ldfp->workscalar));
 65:   for (i = 0; i <= lmvm->k; ++i) {
 66:     stf = ldfp->workscalar[i];
 67:     PetscCall(VecDot(lmvm->Y[i], dX, &ytx));
 68:     /* Update dX_{i+1} = (B^{-1})_{i+1} * f */
 69:     PetscCall(VecAXPBYPCZ(dX, -PetscRealPart(ytx) / ldfp->ytq[i], PetscRealPart(stf) / ldfp->yts[i], 1.0, ldfp->Q[i], lmvm->S[i]));
 70:   }
 71:   PetscFunctionReturn(PETSC_SUCCESS);
 72: }

 74: /*------------------------------------------------------------*/

 76: /*
 77:   The forward product for the approximate Jacobian is the matrix-free
 78:   implementation of the recursive formula given in Equation 6.13 of
 79:   Nocedal and Wright "Numerical Optimization" 2nd edition, pg 139.

 81:   This forward product has a two-loop form similar to the BFGS two-loop
 82:   formulation for the inverse Jacobian application. However, the S and
 83:   Y vectors have interchanged roles.

 85:   work <- X

 87:   for i = k,k-1,k-2,...,0
 88:     rho[i] = 1 / (Y[i]^T S[i])
 89:     alpha[i] = rho[i] * (Y[i]^T work)
 90:     work <- work - (alpha[i] * S[i])
 91:   end

 93:   Z <- J0 * work

 95:   for i = 0,1,2,...,k
 96:     beta = rho[i] * (S[i]^T Y)
 97:     Z <- Z + ((alpha[i] - beta) * Y[i])
 98:   end
 99: */
100: PetscErrorCode MatMult_LMVMDFP(Mat B, Vec X, Vec Z)
101: {
102:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
103:   Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;
104:   PetscInt     i;
105:   PetscReal   *alpha, beta;
106:   PetscScalar  ytx, stz;

108:   PetscFunctionBegin;
109:   /* Copy the function into the work vector for the first loop */
110:   PetscCall(VecCopy(X, ldfp->work));

112:   /* Start the first loop */
113:   PetscCall(PetscMalloc1(lmvm->k + 1, &alpha));
114:   for (i = lmvm->k; i >= 0; --i) {
115:     PetscCall(VecDot(lmvm->Y[i], ldfp->work, &ytx));
116:     alpha[i] = PetscRealPart(ytx) / ldfp->yts[i];
117:     PetscCall(VecAXPY(ldfp->work, -alpha[i], lmvm->S[i]));
118:   }

120:   /* Apply the forward product with initial Jacobian */
121:   PetscCall(MatSymBrdnApplyJ0Fwd(B, ldfp->work, Z));

123:   /* Start the second loop */
124:   for (i = 0; i <= lmvm->k; ++i) {
125:     PetscCall(VecDot(lmvm->S[i], Z, &stz));
126:     beta = PetscRealPart(stz) / ldfp->yts[i];
127:     PetscCall(VecAXPY(Z, alpha[i] - beta, lmvm->Y[i]));
128:   }
129:   PetscCall(PetscFree(alpha));
130:   PetscFunctionReturn(PETSC_SUCCESS);
131: }

133: /*------------------------------------------------------------*/

135: static PetscErrorCode MatUpdate_LMVMDFP(Mat B, Vec X, Vec F)
136: {
137:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
138:   Mat_SymBrdn  *ldfp = (Mat_SymBrdn *)lmvm->ctx;
139:   Mat_LMVM     *dbase;
140:   Mat_DiagBrdn *dctx;
141:   PetscInt      old_k, i;
142:   PetscReal     curvtol, ststmp;
143:   PetscScalar   curvature, ytytmp;

145:   PetscFunctionBegin;
146:   if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
147:   if (lmvm->prev_set) {
148:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
149:     PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
150:     PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));

152:     /* Test if the updates can be accepted */
153:     PetscCall(VecDotNorm2(lmvm->Xprev, lmvm->Fprev, &curvature, &ststmp));
154:     if (ststmp < lmvm->eps) curvtol = 0.0;
155:     else curvtol = lmvm->eps * ststmp;

157:     if (PetscRealPart(curvature) > curvtol) {
158:       /* Update is good, accept it */
159:       ldfp->watchdog = 0;
160:       ldfp->needQ    = PETSC_TRUE;
161:       old_k          = lmvm->k;
162:       PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
163:       /* If we hit the memory limit, shift the yts, yty and sts arrays */
164:       if (old_k == lmvm->k) {
165:         for (i = 0; i <= lmvm->k - 1; ++i) {
166:           ldfp->yts[i] = ldfp->yts[i + 1];
167:           ldfp->yty[i] = ldfp->yty[i + 1];
168:           ldfp->sts[i] = ldfp->sts[i + 1];
169:         }
170:       }
171:       /* Update history of useful scalars */
172:       PetscCall(VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &ytytmp));
173:       ldfp->yts[lmvm->k] = PetscRealPart(curvature);
174:       ldfp->yty[lmvm->k] = PetscRealPart(ytytmp);
175:       ldfp->sts[lmvm->k] = ststmp;
176:       /* Compute the scalar scale if necessary */
177:       if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) PetscCall(MatSymBrdnComputeJ0Scalar(B));
178:     } else {
179:       /* Update is bad, skip it */
180:       ++lmvm->nrejects;
181:       ++ldfp->watchdog;
182:     }
183:   } else {
184:     switch (ldfp->scale_type) {
185:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
186:       dbase = (Mat_LMVM *)ldfp->D->data;
187:       dctx  = (Mat_DiagBrdn *)dbase->ctx;
188:       PetscCall(VecSet(dctx->invD, ldfp->delta));
189:       break;
190:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
191:       ldfp->sigma = ldfp->delta;
192:       break;
193:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
194:       ldfp->sigma = 1.0;
195:       break;
196:     default:
197:       break;
198:     }
199:   }

201:   /* Update the scaling */
202:   if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatLMVMUpdate(ldfp->D, X, F));

204:   if (ldfp->watchdog > ldfp->max_seq_rejects) {
205:     PetscCall(MatLMVMReset(B, PETSC_FALSE));
206:     if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatLMVMReset(ldfp->D, PETSC_FALSE));
207:   }

209:   /* Save the solution and function to be used in the next update */
210:   PetscCall(VecCopy(X, lmvm->Xprev));
211:   PetscCall(VecCopy(F, lmvm->Fprev));
212:   lmvm->prev_set = PETSC_TRUE;
213:   PetscFunctionReturn(PETSC_SUCCESS);
214: }

216: /*------------------------------------------------------------*/

218: static PetscErrorCode MatCopy_LMVMDFP(Mat B, Mat M, MatStructure str)
219: {
220:   Mat_LMVM    *bdata = (Mat_LMVM *)B->data;
221:   Mat_SymBrdn *bctx  = (Mat_SymBrdn *)bdata->ctx;
222:   Mat_LMVM    *mdata = (Mat_LMVM *)M->data;
223:   Mat_SymBrdn *mctx  = (Mat_SymBrdn *)mdata->ctx;
224:   PetscInt     i;

226:   PetscFunctionBegin;
227:   mctx->needQ = bctx->needQ;
228:   for (i = 0; i <= bdata->k; ++i) {
229:     mctx->ytq[i] = bctx->ytq[i];
230:     mctx->yts[i] = bctx->yts[i];
231:     PetscCall(VecCopy(bctx->Q[i], mctx->Q[i]));
232:   }
233:   mctx->scale_type      = bctx->scale_type;
234:   mctx->alpha           = bctx->alpha;
235:   mctx->beta            = bctx->beta;
236:   mctx->rho             = bctx->rho;
237:   mctx->sigma_hist      = bctx->sigma_hist;
238:   mctx->watchdog        = bctx->watchdog;
239:   mctx->max_seq_rejects = bctx->max_seq_rejects;
240:   switch (bctx->scale_type) {
241:   case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
242:     mctx->sigma = bctx->sigma;
243:     break;
244:   case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
245:     PetscCall(MatCopy(bctx->D, mctx->D, SAME_NONZERO_PATTERN));
246:     break;
247:   case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
248:     mctx->sigma = 1.0;
249:     break;
250:   default:
251:     break;
252:   }
253:   PetscFunctionReturn(PETSC_SUCCESS);
254: }

256: /*------------------------------------------------------------*/

258: static PetscErrorCode MatReset_LMVMDFP(Mat B, PetscBool destructive)
259: {
260:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
261:   Mat_SymBrdn  *ldfp = (Mat_SymBrdn *)lmvm->ctx;
262:   Mat_LMVM     *dbase;
263:   Mat_DiagBrdn *dctx;

265:   PetscFunctionBegin;
266:   ldfp->watchdog = 0;
267:   ldfp->needQ    = PETSC_TRUE;
268:   if (ldfp->allocated) {
269:     if (destructive) {
270:       PetscCall(VecDestroy(&ldfp->work));
271:       PetscCall(PetscFree5(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts, ldfp->workscalar));
272:       PetscCall(VecDestroyVecs(lmvm->m, &ldfp->Q));
273:       switch (ldfp->scale_type) {
274:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
275:         PetscCall(MatLMVMReset(ldfp->D, PETSC_TRUE));
276:         break;
277:       default:
278:         break;
279:       }
280:       ldfp->allocated = PETSC_FALSE;
281:     } else {
282:       switch (ldfp->scale_type) {
283:       case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
284:         ldfp->sigma = ldfp->delta;
285:         break;
286:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
287:         PetscCall(MatLMVMReset(ldfp->D, PETSC_FALSE));
288:         dbase = (Mat_LMVM *)ldfp->D->data;
289:         dctx  = (Mat_DiagBrdn *)dbase->ctx;
290:         PetscCall(VecSet(dctx->invD, ldfp->delta));
291:         break;
292:       case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
293:         ldfp->sigma = 1.0;
294:         break;
295:       default:
296:         break;
297:       }
298:     }
299:   }
300:   PetscCall(MatReset_LMVM(B, destructive));
301:   PetscFunctionReturn(PETSC_SUCCESS);
302: }

304: /*------------------------------------------------------------*/

306: static PetscErrorCode MatAllocate_LMVMDFP(Mat B, Vec X, Vec F)
307: {
308:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
309:   Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;

311:   PetscFunctionBegin;
312:   PetscCall(MatAllocate_LMVM(B, X, F));
313:   if (!ldfp->allocated) {
314:     PetscCall(VecDuplicate(X, &ldfp->work));
315:     PetscCall(PetscMalloc5(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts, lmvm->m, &ldfp->workscalar));
316:     if (lmvm->m > 0) PetscCall(VecDuplicateVecs(X, lmvm->m, &ldfp->Q));
317:     switch (ldfp->scale_type) {
318:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
319:       PetscCall(MatLMVMAllocate(ldfp->D, X, F));
320:       break;
321:     default:
322:       break;
323:     }
324:     ldfp->allocated = PETSC_TRUE;
325:   }
326:   PetscFunctionReturn(PETSC_SUCCESS);
327: }

329: /*------------------------------------------------------------*/

331: static PetscErrorCode MatDestroy_LMVMDFP(Mat B)
332: {
333:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
334:   Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;

336:   PetscFunctionBegin;
337:   if (ldfp->allocated) {
338:     PetscCall(VecDestroy(&ldfp->work));
339:     PetscCall(PetscFree5(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts, ldfp->workscalar));
340:     PetscCall(VecDestroyVecs(lmvm->m, &ldfp->Q));
341:     ldfp->allocated = PETSC_FALSE;
342:   }
343:   PetscCall(MatDestroy(&ldfp->D));
344:   PetscCall(PetscFree(lmvm->ctx));
345:   PetscCall(MatDestroy_LMVM(B));
346:   PetscFunctionReturn(PETSC_SUCCESS);
347: }

349: /*------------------------------------------------------------*/

351: static PetscErrorCode MatSetUp_LMVMDFP(Mat B)
352: {
353:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
354:   Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;
355:   PetscInt     n, N;

357:   PetscFunctionBegin;
358:   PetscCall(MatSetUp_LMVM(B));
359:   if (!ldfp->allocated) {
360:     PetscCall(VecDuplicate(lmvm->Xprev, &ldfp->work));
361:     PetscCall(PetscMalloc5(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts, lmvm->m, &ldfp->workscalar));
362:     if (lmvm->m > 0) PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &ldfp->Q));
363:     switch (ldfp->scale_type) {
364:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
365:       PetscCall(MatGetLocalSize(B, &n, &n));
366:       PetscCall(MatGetSize(B, &N, &N));
367:       PetscCall(MatSetSizes(ldfp->D, n, n, N, N));
368:       PetscCall(MatSetUp(ldfp->D));
369:       break;
370:     default:
371:       break;
372:     }
373:     ldfp->allocated = PETSC_TRUE;
374:   }
375:   PetscFunctionReturn(PETSC_SUCCESS);
376: }

378: /*------------------------------------------------------------*/

380: static PetscErrorCode MatSetFromOptions_LMVMDFP(Mat B, PetscOptionItems *PetscOptionsObject)
381: {
382:   PetscFunctionBegin;
383:   PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
384:   PetscOptionsHeadBegin(PetscOptionsObject, "DFP method for approximating SPD Jacobian actions (MATLMVMDFP)");
385:   PetscCall(MatSetFromOptions_LMVMSymBrdn_Private(B, PetscOptionsObject));
386:   PetscOptionsHeadEnd();
387:   PetscFunctionReturn(PETSC_SUCCESS);
388: }

390: /*------------------------------------------------------------*/

392: PetscErrorCode MatCreate_LMVMDFP(Mat B)
393: {
394:   Mat_LMVM    *lmvm;
395:   Mat_SymBrdn *ldfp;

397:   PetscFunctionBegin;
398:   PetscCall(MatCreate_LMVMSymBrdn(B));
399:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDFP));
400:   B->ops->setup          = MatSetUp_LMVMDFP;
401:   B->ops->destroy        = MatDestroy_LMVMDFP;
402:   B->ops->setfromoptions = MatSetFromOptions_LMVMDFP;
403:   B->ops->solve          = MatSolve_LMVMDFP;

405:   lmvm                = (Mat_LMVM *)B->data;
406:   lmvm->ops->allocate = MatAllocate_LMVMDFP;
407:   lmvm->ops->reset    = MatReset_LMVMDFP;
408:   lmvm->ops->update   = MatUpdate_LMVMDFP;
409:   lmvm->ops->mult     = MatMult_LMVMDFP;
410:   lmvm->ops->copy     = MatCopy_LMVMDFP;

412:   ldfp        = (Mat_SymBrdn *)lmvm->ctx;
413:   ldfp->needP = PETSC_FALSE;
414:   ldfp->phi   = 1.0;
415:   PetscFunctionReturn(PETSC_SUCCESS);
416: }

418: /*------------------------------------------------------------*/

420: /*@
421:    MatCreateLMVMDFP - Creates a limited-memory Davidon-Fletcher-Powell (DFP) matrix
422:    used for approximating Jacobians. L-DFP is symmetric positive-definite by
423:    construction, and is the dual of L-BFGS where Y and S vectors swap roles.

425:    To use the L-DFP matrix with other vector types, the matrix must be
426:    created using `MatCreate()` and `MatSetType()`, followed by `MatLMVMAllocate()`.
427:    This ensures that the internal storage and work vectors are duplicated from the
428:    correct type of vector.

430:    Collective

432:    Input Parameters:
433: +  comm - MPI communicator
434: .  n - number of local rows for storage vectors
435: -  N - global size of the storage vectors

437:    Output Parameter:
438: .  B - the matrix

440:    Options Database Keys:
441: +   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
442: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
443: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
444: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
445: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
446: -   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

448:    Level: intermediate

450:    Note:
451:    It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
452:    paradigm instead of this routine directly.

454: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMDFP`, `MatCreateLMVMBFGS()`, `MatCreateLMVMSR1()`,
455:           `MatCreateLMVMBrdn()`, `MatCreateLMVMBadBrdn()`, `MatCreateLMVMSymBrdn()`
456: @*/
457: PetscErrorCode MatCreateLMVMDFP(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
458: {
459:   PetscFunctionBegin;
460:   PetscCall(MatCreate(comm, B));
461:   PetscCall(MatSetSizes(*B, n, n, N, N));
462:   PetscCall(MatSetType(*B, MATLMVMDFP));
463:   PetscCall(MatSetUp(*B));
464:   PetscFunctionReturn(PETSC_SUCCESS);
465: }