Actual source code: owarmijo.c


  2: #include <petsc/private/taolinesearchimpl.h>
  3: #include <../src/tao/linesearch/impls/owarmijo/owarmijo.h>

  5: #define REPLACE_FIFO 1
  6: #define REPLACE_MRU  2

  8: #define REFERENCE_MAX  1
  9: #define REFERENCE_AVE  2
 10: #define REFERENCE_MEAN 3

 12: static PetscErrorCode ProjWork_OWLQN(Vec w, Vec x, Vec gv, PetscReal *gdx)
 13: {
 14:   const PetscReal *xptr, *gptr;
 15:   PetscReal       *wptr;
 16:   PetscInt         low, high, low1, high1, low2, high2, i;

 18:   VecGetOwnershipRange(w, &low, &high);
 19:   VecGetOwnershipRange(x, &low1, &high1);
 20:   VecGetOwnershipRange(gv, &low2, &high2);

 22:   *gdx = 0.0;
 23:   VecGetArray(w, &wptr);
 24:   VecGetArrayRead(x, &xptr);
 25:   VecGetArrayRead(gv, &gptr);

 27:   for (i = 0; i < high - low; i++) {
 28:     if (xptr[i] * wptr[i] < 0.0) wptr[i] = 0.0;
 29:     *gdx = *gdx + gptr[i] * (wptr[i] - xptr[i]);
 30:   }
 31:   VecRestoreArray(w, &wptr);
 32:   VecRestoreArrayRead(x, &xptr);
 33:   VecRestoreArrayRead(gv, &gptr);
 34:   return 0;
 35: }

 37: static PetscErrorCode TaoLineSearchDestroy_OWArmijo(TaoLineSearch ls)
 38: {
 39:   TaoLineSearch_OWARMIJO *armP = (TaoLineSearch_OWARMIJO *)ls->data;

 41:   PetscFree(armP->memory);
 42:   if (armP->x) PetscObjectDereference((PetscObject)armP->x);
 43:   VecDestroy(&armP->work);
 44:   PetscFree(ls->data);
 45:   return 0;
 46: }

 48: static PetscErrorCode TaoLineSearchSetFromOptions_OWArmijo(TaoLineSearch ls, PetscOptionItems *PetscOptionsObject)
 49: {
 50:   TaoLineSearch_OWARMIJO *armP = (TaoLineSearch_OWARMIJO *)ls->data;

 52:   PetscOptionsHeadBegin(PetscOptionsObject, "OWArmijo linesearch options");
 53:   PetscOptionsReal("-tao_ls_OWArmijo_alpha", "initial reference constant", "", armP->alpha, &armP->alpha, NULL);
 54:   PetscOptionsReal("-tao_ls_OWArmijo_beta_inf", "decrease constant one", "", armP->beta_inf, &armP->beta_inf, NULL);
 55:   PetscOptionsReal("-tao_ls_OWArmijo_beta", "decrease constant", "", armP->beta, &armP->beta, NULL);
 56:   PetscOptionsReal("-tao_ls_OWArmijo_sigma", "acceptance constant", "", armP->sigma, &armP->sigma, NULL);
 57:   PetscOptionsInt("-tao_ls_OWArmijo_memory_size", "number of historical elements", "", armP->memorySize, &armP->memorySize, NULL);
 58:   PetscOptionsInt("-tao_ls_OWArmijo_reference_policy", "policy for updating reference value", "", armP->referencePolicy, &armP->referencePolicy, NULL);
 59:   PetscOptionsInt("-tao_ls_OWArmijo_replacement_policy", "policy for updating memory", "", armP->replacementPolicy, &armP->replacementPolicy, NULL);
 60:   PetscOptionsBool("-tao_ls_OWArmijo_nondescending", "Use nondescending OWArmijo algorithm", "", armP->nondescending, &armP->nondescending, NULL);
 61:   PetscOptionsHeadEnd();
 62:   return 0;
 63: }

 65: static PetscErrorCode TaoLineSearchView_OWArmijo(TaoLineSearch ls, PetscViewer pv)
 66: {
 67:   TaoLineSearch_OWARMIJO *armP = (TaoLineSearch_OWARMIJO *)ls->data;
 68:   PetscBool               isascii;

 70:   PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii);
 71:   if (isascii) {
 72:     PetscViewerASCIIPrintf(pv, "  OWArmijo linesearch");
 73:     if (armP->nondescending) PetscViewerASCIIPrintf(pv, " (nondescending)");
 74:     PetscViewerASCIIPrintf(pv, ": alpha=%g beta=%g ", (double)armP->alpha, (double)armP->beta);
 75:     PetscViewerASCIIPrintf(pv, "sigma=%g ", (double)armP->sigma);
 76:     PetscViewerASCIIPrintf(pv, "memsize=%" PetscInt_FMT "\n", armP->memorySize);
 77:   }
 78:   return 0;
 79: }

 81: /* @ TaoApply_OWArmijo - This routine performs a linesearch. It
 82:    backtracks until the (nonmonotone) OWArmijo conditions are satisfied.

 84:    Input Parameters:
 85: +  tao - TAO_SOLVER context
 86: .  X - current iterate (on output X contains new iterate, X + step*S)
 87: .  S - search direction
 88: .  f - merit function evaluated at X
 89: .  G - gradient of merit function evaluated at X
 90: .  W - work vector
 91: -  step - initial estimate of step length

 93:    Output parameters:
 94: +  f - merit function evaluated at new iterate, X + step*S
 95: .  G - gradient of merit function evaluated at new iterate, X + step*S
 96: .  X - new iterate
 97: -  step - final step length

 99:    Info is set to one of:
100: .   0 - the line search succeeds; the sufficient decrease
101:    condition and the directional derivative condition hold

103:    negative number if an input parameter is invalid
104: -   -1 -  step < 0

106:    positive number > 1 if the line search otherwise terminates
107: +    1 -  Step is at the lower bound, stepmin.
108: @ */
109: static PetscErrorCode TaoLineSearchApply_OWArmijo(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s)
110: {
111:   TaoLineSearch_OWARMIJO *armP = (TaoLineSearch_OWARMIJO *)ls->data;
112:   PetscInt                i, its = 0;
113:   PetscReal               fact, ref, gdx;
114:   PetscInt                idx;
115:   PetscBool               g_computed = PETSC_FALSE; /* to prevent extra gradient computation */
116:   Vec                     g_old;
117:   PetscReal               owlqn_minstep = 0.005;
118:   PetscReal               partgdx;
119:   MPI_Comm                comm;

121:   PetscObjectGetComm((PetscObject)ls, &comm);
122:   fact       = 0.0;
123:   ls->nfeval = 0;
124:   ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
125:   if (!armP->work) {
126:     VecDuplicate(x, &armP->work);
127:     armP->x = x;
128:     PetscObjectReference((PetscObject)armP->x);
129:   } else if (x != armP->x) {
130:     VecDestroy(&armP->work);
131:     VecDuplicate(x, &armP->work);
132:     PetscObjectDereference((PetscObject)armP->x);
133:     armP->x = x;
134:     PetscObjectReference((PetscObject)armP->x);
135:   }

137:   TaoLineSearchMonitor(ls, 0, *f, 0.0);

139:   /* Check linesearch parameters */
140:   if (armP->alpha < 1) {
141:     PetscInfo(ls, "OWArmijo line search error: alpha (%g) < 1\n", (double)armP->alpha);
142:     ls->reason = TAOLINESEARCH_FAILED_BADPARAMETER;
143:   } else if ((armP->beta <= 0) || (armP->beta >= 1)) {
144:     PetscInfo(ls, "OWArmijo line search error: beta (%g) invalid\n", (double)armP->beta);
145:     ls->reason = TAOLINESEARCH_FAILED_BADPARAMETER;
146:   } else if ((armP->beta_inf <= 0) || (armP->beta_inf >= 1)) {
147:     PetscInfo(ls, "OWArmijo line search error: beta_inf (%g) invalid\n", (double)armP->beta_inf);
148:     ls->reason = TAOLINESEARCH_FAILED_BADPARAMETER;
149:   } else if ((armP->sigma <= 0) || (armP->sigma >= 0.5)) {
150:     PetscInfo(ls, "OWArmijo line search error: sigma (%g) invalid\n", (double)armP->sigma);
151:     ls->reason = TAOLINESEARCH_FAILED_BADPARAMETER;
152:   } else if (armP->memorySize < 1) {
153:     PetscInfo(ls, "OWArmijo line search error: memory_size (%" PetscInt_FMT ") < 1\n", armP->memorySize);
154:     ls->reason = TAOLINESEARCH_FAILED_BADPARAMETER;
155:   } else if ((armP->referencePolicy != REFERENCE_MAX) && (armP->referencePolicy != REFERENCE_AVE) && (armP->referencePolicy != REFERENCE_MEAN)) {
156:     PetscInfo(ls, "OWArmijo line search error: reference_policy invalid\n");
157:     ls->reason = TAOLINESEARCH_FAILED_BADPARAMETER;
158:   } else if ((armP->replacementPolicy != REPLACE_FIFO) && (armP->replacementPolicy != REPLACE_MRU)) {
159:     PetscInfo(ls, "OWArmijo line search error: replacement_policy invalid\n");
160:     ls->reason = TAOLINESEARCH_FAILED_BADPARAMETER;
161:   } else if (PetscIsInfOrNanReal(*f)) {
162:     PetscInfo(ls, "OWArmijo line search error: initial function inf or nan\n");
163:     ls->reason = TAOLINESEARCH_FAILED_BADPARAMETER;
164:   }

166:   if (ls->reason != TAOLINESEARCH_CONTINUE_ITERATING) return 0;

168:   /* Check to see of the memory has been allocated.  If not, allocate
169:      the historical array and populate it with the initial function
170:      values. */
171:   if (!armP->memory) PetscMalloc1(armP->memorySize, &armP->memory);

173:   if (!armP->memorySetup) {
174:     for (i = 0; i < armP->memorySize; i++) armP->memory[i] = armP->alpha * (*f);
175:     armP->current       = 0;
176:     armP->lastReference = armP->memory[0];
177:     armP->memorySetup   = PETSC_TRUE;
178:   }

180:   /* Calculate reference value (MAX) */
181:   ref = armP->memory[0];
182:   idx = 0;

184:   for (i = 1; i < armP->memorySize; i++) {
185:     if (armP->memory[i] > ref) {
186:       ref = armP->memory[i];
187:       idx = i;
188:     }
189:   }

191:   if (armP->referencePolicy == REFERENCE_AVE) {
192:     ref = 0;
193:     for (i = 0; i < armP->memorySize; i++) ref += armP->memory[i];
194:     ref = ref / armP->memorySize;
195:     ref = PetscMax(ref, armP->memory[armP->current]);
196:   } else if (armP->referencePolicy == REFERENCE_MEAN) {
197:     ref = PetscMin(ref, 0.5 * (armP->lastReference + armP->memory[armP->current]));
198:   }

200:   if (armP->nondescending) fact = armP->sigma;

202:   VecDuplicate(g, &g_old);
203:   VecCopy(g, g_old);

205:   ls->step = ls->initstep;
206:   while (ls->step >= owlqn_minstep && ls->nfeval < ls->max_funcs) {
207:     /* Calculate iterate */
208:     ++its;
209:     VecWAXPY(armP->work, ls->step, s, x);

211:     partgdx = 0.0;
212:     ProjWork_OWLQN(armP->work, x, g_old, &partgdx);
213:     MPIU_Allreduce(&partgdx, &gdx, 1, MPIU_REAL, MPIU_SUM, comm);

215:     /* Check the condition of gdx */
216:     if (PetscIsInfOrNanReal(gdx)) {
217:       PetscInfo(ls, "Initial Line Search step * g is Inf or Nan (%g)\n", (double)gdx);
218:       ls->reason = TAOLINESEARCH_FAILED_INFORNAN;
219:       return 0;
220:     }
221:     if (gdx >= 0.0) {
222:       PetscInfo(ls, "Initial Line Search step is not descent direction (g's=%g)\n", (double)gdx);
223:       ls->reason = TAOLINESEARCH_FAILED_ASCENT;
224:       return 0;
225:     }

227:     /* Calculate function at new iterate */
228:     TaoLineSearchComputeObjectiveAndGradient(ls, armP->work, f, g);
229:     g_computed = PETSC_TRUE;

231:     TaoLineSearchMonitor(ls, its, *f, ls->step);

233:     if (ls->step == ls->initstep) ls->f_fullstep = *f;

235:     if (PetscIsInfOrNanReal(*f)) {
236:       ls->step *= armP->beta_inf;
237:     } else {
238:       /* Check descent condition */
239:       if (armP->nondescending && *f <= ref - ls->step * fact * ref) break;
240:       if (!armP->nondescending && *f <= ref + armP->sigma * gdx) break;
241:       ls->step *= armP->beta;
242:     }
243:   }
244:   VecDestroy(&g_old);

246:   /* Check termination */
247:   if (PetscIsInfOrNanReal(*f)) {
248:     PetscInfo(ls, "Function is inf or nan.\n");
249:     ls->reason = TAOLINESEARCH_FAILED_BADPARAMETER;
250:   } else if (ls->step < owlqn_minstep) {
251:     PetscInfo(ls, "Step length is below tolerance.\n");
252:     ls->reason = TAOLINESEARCH_HALTED_RTOL;
253:   } else if (ls->nfeval >= ls->max_funcs) {
254:     PetscInfo(ls, "Number of line search function evals (%" PetscInt_FMT ") > maximum allowed (%" PetscInt_FMT ")\n", ls->nfeval, ls->max_funcs);
255:     ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
256:   }
257:   if (ls->reason) return 0;

259:   /* Successful termination, update memory */
260:   ls->reason          = TAOLINESEARCH_SUCCESS;
261:   armP->lastReference = ref;
262:   if (armP->replacementPolicy == REPLACE_FIFO) {
263:     armP->memory[armP->current++] = *f;
264:     if (armP->current >= armP->memorySize) armP->current = 0;
265:   } else {
266:     armP->current     = idx;
267:     armP->memory[idx] = *f;
268:   }

270:   /* Update iterate and compute gradient */
271:   VecCopy(armP->work, x);
272:   if (!g_computed) TaoLineSearchComputeGradient(ls, x, g);
273:   PetscInfo(ls, "%" PetscInt_FMT " function evals in line search, step = %10.4f\n", ls->nfeval, (double)ls->step);
274:   return 0;
275: }

277: /*MC
278:    TAOLINESEARCHOWARMIJO - Special line-search type for the Orthant-Wise Limited Quasi-Newton (TAOOWLQN) algorithm.
279:    Should not be used with any other algorithm.

281:    Level: developer

283: .keywords: Tao, linesearch
284: M*/
285: PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_OWArmijo(TaoLineSearch ls)
286: {
287:   TaoLineSearch_OWARMIJO *armP;

290:   PetscNew(&armP);

292:   armP->memory            = NULL;
293:   armP->alpha             = 1.0;
294:   armP->beta              = 0.25;
295:   armP->beta_inf          = 0.25;
296:   armP->sigma             = 1e-4;
297:   armP->memorySize        = 1;
298:   armP->referencePolicy   = REFERENCE_MAX;
299:   armP->replacementPolicy = REPLACE_MRU;
300:   armP->nondescending     = PETSC_FALSE;
301:   ls->data                = (void *)armP;
302:   ls->initstep            = 0.1;
303:   ls->ops->monitor        = NULL;
304:   ls->ops->setup          = NULL;
305:   ls->ops->reset          = NULL;
306:   ls->ops->apply          = TaoLineSearchApply_OWArmijo;
307:   ls->ops->view           = TaoLineSearchView_OWArmijo;
308:   ls->ops->destroy        = TaoLineSearchDestroy_OWArmijo;
309:   ls->ops->setfromoptions = TaoLineSearchSetFromOptions_OWArmijo;
310:   return 0;
311: }