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;

 39:   VecCheckSameSize(F, 2, dX, 3);
 40:   VecCheckMatCompatible(B, dX, 3, F, 2);

 42:   if (ldfp->needQ) {
 43:     /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
 44:     for (i = 0; i <= lmvm->k; ++i) {
 45:       MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], ldfp->Q[i]);
 46:       for (j = 0; j <= i-1; ++j) {
 47:         /* Compute the necessary dot products */
 48:         VecDotBegin(lmvm->Y[j], ldfp->Q[i], &yjtqi);
 49:         VecDotBegin(lmvm->S[j], lmvm->Y[i], &sjtyi);
 50:         VecDotEnd(lmvm->Y[j], ldfp->Q[i], &yjtqi);
 51:         VecDotEnd(lmvm->S[j], lmvm->Y[i], &sjtyi);
 52:         /* Compute the pure DFP component of the inverse application*/
 53:         VecAXPBYPCZ(ldfp->Q[i], -PetscRealPart(yjtqi)/ldfp->ytq[j], PetscRealPart(sjtyi)/ldfp->yts[j], 1.0, ldfp->Q[j], lmvm->S[j]);
 54:       }
 55:       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:   MatSymBrdnApplyJ0Inv(B, F, dX);
 63:   for (i = 0; i <= lmvm->k; ++i) {
 64:     /* Get all the dot products we need */
 65:     VecDotBegin(lmvm->Y[i], dX, &ytx);
 66:     VecDotBegin(lmvm->S[i], F, &stf);
 67:     VecDotEnd(lmvm->Y[i], dX, &ytx);
 68:     VecDotEnd(lmvm->S[i], F, &stf);
 69:     /* Update dX_{i+1} = (B^{-1})_{i+1} * f */
 70:     VecAXPBYPCZ(dX, -PetscRealPart(ytx)/ldfp->ytq[i], PetscRealPart(stf)/ldfp->yts[i], 1.0, ldfp->Q[i], lmvm->S[i]);
 71:   }
 72:   return 0;
 73: }

 75: /*------------------------------------------------------------*/

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

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

 86:   work <- X

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

 94:   Z <- J0 * work

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

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

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

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

123:   /* Start the second loop */
124:   for (i = 0; i <= lmvm->k; ++i) {
125:     VecDot(lmvm->S[i], Z, &stz);
126:     beta = PetscRealPart(stz)/ldfp->yts[i];
127:     VecAXPY(Z, alpha[i]-beta, lmvm->Y[i]);
128:   }
129:   PetscFree(alpha);
130:   return 0;
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;
143:   PetscScalar       curvature, ytytmp, ststmp;

145:   if (!lmvm->m) return 0;
146:   if (lmvm->prev_set) {
147:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
148:     VecAYPX(lmvm->Xprev, -1.0, X);
149:     VecAYPX(lmvm->Fprev, -1.0, F);
150:     /* Test if the updates can be accepted */
151:     VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
152:     VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
153:     VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
154:     VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
155:     if (PetscRealPart(ststmp) < lmvm->eps) {
156:       curvtol = 0.0;
157:     } else {
158:       curvtol = lmvm->eps * PetscRealPart(ststmp);
159:     }
160:     if (PetscRealPart(curvature) > curvtol) {
161:       /* Update is good, accept it */
162:       ldfp->watchdog = 0;
163:       ldfp->needQ = PETSC_TRUE;
164:       old_k = lmvm->k;
165:       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:           ldfp->yts[i] = ldfp->yts[i+1];
170:           ldfp->yty[i] = ldfp->yty[i+1];
171:           ldfp->sts[i] = ldfp->sts[i+1];
172:         }
173:       }
174:       /* Update history of useful scalars */
175:       VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &ytytmp);
176:       ldfp->yts[lmvm->k] = PetscRealPart(curvature);
177:       ldfp->yty[lmvm->k] = PetscRealPart(ytytmp);
178:       ldfp->sts[lmvm->k] = PetscRealPart(ststmp);
179:       /* Compute the scalar scale if necessary */
180:       if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) {
181:         MatSymBrdnComputeJ0Scalar(B);
182:       }
183:     } else {
184:       /* Update is bad, skip it */
185:       ++lmvm->nrejects;
186:       ++ldfp->watchdog;
187:     }
188:   } else {
189:     switch (ldfp->scale_type) {
190:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
191:       dbase = (Mat_LMVM*)ldfp->D->data;
192:       dctx = (Mat_DiagBrdn*)dbase->ctx;
193:       VecSet(dctx->invD, ldfp->delta);
194:       break;
195:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
196:       ldfp->sigma = ldfp->delta;
197:       break;
198:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
199:       ldfp->sigma = 1.0;
200:       break;
201:     default:
202:       break;
203:     }
204:   }

206:   /* Update the scaling */
207:   if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
208:     MatLMVMUpdate(ldfp->D, X, F);
209:   }

211:   if (ldfp->watchdog > ldfp->max_seq_rejects) {
212:     MatLMVMReset(B, PETSC_FALSE);
213:     if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
214:       MatLMVMReset(ldfp->D, PETSC_FALSE);
215:     }
216:   }

218:   /* Save the solution and function to be used in the next update */
219:   VecCopy(X, lmvm->Xprev);
220:   VecCopy(F, lmvm->Fprev);
221:   lmvm->prev_set = PETSC_TRUE;
222:   return 0;
223: }

225: /*------------------------------------------------------------*/

227: static PetscErrorCode MatCopy_LMVMDFP(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:   mctx->needQ = bctx->needQ;
236:   for (i=0; i<=bdata->k; ++i) {
237:     mctx->ytq[i] = bctx->ytq[i];
238:     mctx->yts[i] = bctx->yts[i];
239:     VecCopy(bctx->Q[i], mctx->Q[i]);
240:   }
241:   mctx->scale_type      = bctx->scale_type;
242:   mctx->alpha           = bctx->alpha;
243:   mctx->beta            = bctx->beta;
244:   mctx->rho             = bctx->rho;
245:   mctx->sigma_hist      = bctx->sigma_hist;
246:   mctx->watchdog        = bctx->watchdog;
247:   mctx->max_seq_rejects = bctx->max_seq_rejects;
248:   switch (bctx->scale_type) {
249:   case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
250:     mctx->sigma = bctx->sigma;
251:     break;
252:   case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
253:     MatCopy(bctx->D, mctx->D, SAME_NONZERO_PATTERN);
254:     break;
255:   case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
256:     mctx->sigma = 1.0;
257:     break;
258:   default:
259:     break;
260:   }
261:   return 0;
262: }

264: /*------------------------------------------------------------*/

266: static PetscErrorCode MatReset_LMVMDFP(Mat B, PetscBool destructive)
267: {
268:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
269:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
270:   Mat_LMVM          *dbase;
271:   Mat_DiagBrdn      *dctx;

273:   ldfp->watchdog = 0;
274:   ldfp->needQ = PETSC_TRUE;
275:   if (ldfp->allocated) {
276:     if (destructive) {
277:       VecDestroy(&ldfp->work);
278:       PetscFree4(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts);
279:       VecDestroyVecs(lmvm->m, &ldfp->Q);
280:       switch (ldfp->scale_type) {
281:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
282:         MatLMVMReset(ldfp->D, PETSC_TRUE);
283:         break;
284:       default:
285:         break;
286:       }
287:       ldfp->allocated = PETSC_FALSE;
288:     } else {
289:       switch (ldfp->scale_type) {
290:       case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
291:         ldfp->sigma = ldfp->delta;
292:         break;
293:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
294:         MatLMVMReset(ldfp->D, PETSC_FALSE);
295:         dbase = (Mat_LMVM*)ldfp->D->data;
296:         dctx = (Mat_DiagBrdn*)dbase->ctx;
297:         VecSet(dctx->invD, ldfp->delta);
298:         break;
299:       case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
300:         ldfp->sigma = 1.0;
301:         break;
302:       default:
303:         break;
304:       }
305:     }
306:   }
307:   MatReset_LMVM(B, destructive);
308:   return 0;
309: }

311: /*------------------------------------------------------------*/

313: static PetscErrorCode MatAllocate_LMVMDFP(Mat B, Vec X, Vec F)
314: {
315:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
316:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;

318:   MatAllocate_LMVM(B, X, F);
319:   if (!ldfp->allocated) {
320:     VecDuplicate(X, &ldfp->work);
321:     PetscMalloc4(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts);
322:     if (lmvm->m > 0) {
323:       VecDuplicateVecs(X, lmvm->m, &ldfp->Q);
324:     }
325:     switch (ldfp->scale_type) {
326:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
327:       MatLMVMAllocate(ldfp->D, X, F);
328:       break;
329:     default:
330:       break;
331:     }
332:     ldfp->allocated = PETSC_TRUE;
333:   }
334:   return 0;
335: }

337: /*------------------------------------------------------------*/

339: static PetscErrorCode MatDestroy_LMVMDFP(Mat B)
340: {
341:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
342:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;

344:   if (ldfp->allocated) {
345:     VecDestroy(&ldfp->work);
346:     PetscFree4(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts);
347:     VecDestroyVecs(lmvm->m, &ldfp->Q);
348:     ldfp->allocated = PETSC_FALSE;
349:   }
350:   MatDestroy(&ldfp->D);
351:   PetscFree(lmvm->ctx);
352:   MatDestroy_LMVM(B);
353:   return 0;
354: }

356: /*------------------------------------------------------------*/

358: static PetscErrorCode MatSetUp_LMVMDFP(Mat B)
359: {
360:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
361:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
362:   PetscInt          n, N;

364:   MatSetUp_LMVM(B);
365:   if (!ldfp->allocated) {
366:     VecDuplicate(lmvm->Xprev, &ldfp->work);
367:     PetscMalloc4(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts);
368:     if (lmvm->m > 0) {
369:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &ldfp->Q);
370:     }
371:     switch (ldfp->scale_type) {
372:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
373:       MatGetLocalSize(B, &n, &n);
374:       MatGetSize(B, &N, &N);
375:       MatSetSizes(ldfp->D, n, n, N, N);
376:       MatSetUp(ldfp->D);
377:       break;
378:     default:
379:       break;
380:     }
381:     ldfp->allocated = PETSC_TRUE;
382:   }
383:   return 0;
384: }

386: /*------------------------------------------------------------*/

388: static PetscErrorCode MatSetFromOptions_LMVMDFP(PetscOptionItems *PetscOptionsObject, Mat B)
389: {
390:   MatSetFromOptions_LMVM(PetscOptionsObject, B);
391:   PetscOptionsHead(PetscOptionsObject,"DFP method for approximating SPD Jacobian actions (MATLMVMDFP)");
392:   MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionsObject, B);
393:   PetscOptionsTail();
394:   return 0;
395: }

397: /*------------------------------------------------------------*/

399: PetscErrorCode MatCreate_LMVMDFP(Mat B)
400: {
401:   Mat_LMVM          *lmvm;
402:   Mat_SymBrdn       *ldfp;

404:   MatCreate_LMVMSymBrdn(B);
405:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMDFP);
406:   B->ops->setup = MatSetUp_LMVMDFP;
407:   B->ops->destroy = MatDestroy_LMVMDFP;
408:   B->ops->setfromoptions = MatSetFromOptions_LMVMDFP;
409:   B->ops->solve = MatSolve_LMVMDFP;

411:   lmvm = (Mat_LMVM*)B->data;
412:   lmvm->ops->allocate = MatAllocate_LMVMDFP;
413:   lmvm->ops->reset = MatReset_LMVMDFP;
414:   lmvm->ops->update = MatUpdate_LMVMDFP;
415:   lmvm->ops->mult = MatMult_LMVMDFP;
416:   lmvm->ops->copy = MatCopy_LMVMDFP;

418:   ldfp = (Mat_SymBrdn*)lmvm->ctx;
419:   ldfp->needP           = PETSC_FALSE;
420:   ldfp->phi             = 1.0;
421:   return 0;
422: }

424: /*------------------------------------------------------------*/

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

431:    The provided local and global sizes must match the solution and function vectors
432:    used with MatLMVMUpdate() and MatSolve(). The resulting L-DFP matrix will have
433:    storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
434:    parallel. To use the L-DFP matrix with other vector types, the matrix must be
435:    created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
436:    This ensures that the internal storage and work vectors are duplicated from the
437:    correct type of vector.

439:    Collective

441:    Input Parameters:
442: +  comm - MPI communicator, set to PETSC_COMM_SELF
443: .  n - number of local rows for storage vectors
444: -  N - global size of the storage vectors

446:    Output Parameter:
447: .  B - the matrix

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

452:    Options Database Keys:
453: +   -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
454: .   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
455: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
456: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
457: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
458: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
459: -   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

461:    Level: intermediate

463: .seealso: MatCreate(), MATLMVM, MATLMVMDFP, MatCreateLMVMBFGS(), MatCreateLMVMSR1(),
464:            MatCreateLMVMBrdn(), MatCreateLMVMBadBrdn(), MatCreateLMVMSymBrdn()
465: @*/
466: PetscErrorCode MatCreateLMVMDFP(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
467: {
468:   MatCreate(comm, B);
469:   MatSetSizes(*B, n, n, N, N);
470:   MatSetType(*B, MATLMVMDFP);
471:   MatSetUp(*B);
472:   return 0;
473: }