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:   PetscFunctionBegin;
 19:   PetscCall(VecGetOwnershipRange(w, &low, &high));
 20:   PetscCall(VecGetOwnershipRange(x, &low1, &high1));
 21:   PetscCall(VecGetOwnershipRange(gv, &low2, &high2));

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

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

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

 42:   PetscFunctionBegin;
 43:   PetscCall(PetscFree(armP->memory));
 44:   if (armP->x) PetscCall(PetscObjectDereference((PetscObject)armP->x));
 45:   PetscCall(VecDestroy(&armP->work));
 46:   PetscCall(PetscFree(ls->data));
 47:   PetscFunctionReturn(PETSC_SUCCESS);
 48: }

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

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

 68: static PetscErrorCode TaoLineSearchView_OWArmijo(TaoLineSearch ls, PetscViewer pv)
 69: {
 70:   TaoLineSearch_OWARMIJO *armP = (TaoLineSearch_OWARMIJO *)ls->data;
 71:   PetscBool               isascii;

 73:   PetscFunctionBegin;
 74:   PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
 75:   if (isascii) {
 76:     PetscCall(PetscViewerASCIIPrintf(pv, "  OWArmijo linesearch"));
 77:     if (armP->nondescending) PetscCall(PetscViewerASCIIPrintf(pv, " (nondescending)"));
 78:     PetscCall(PetscViewerASCIIPrintf(pv, ": alpha=%g beta=%g ", (double)armP->alpha, (double)armP->beta));
 79:     PetscCall(PetscViewerASCIIPrintf(pv, "sigma=%g ", (double)armP->sigma));
 80:     PetscCall(PetscViewerASCIIPrintf(pv, "memsize=%" PetscInt_FMT "\n", armP->memorySize));
 81:   }
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }

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

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

 97:    Output parameters:
 98: +  f - merit function evaluated at new iterate, X + step*S
 99: .  G - gradient of merit function evaluated at new iterate, X + step*S
100: .  X - new iterate
101: -  step - final step length

103:    Info is set to one of:
104: .   0 - the line search succeeds; the sufficient decrease
105:    condition and the directional derivative condition hold

107:    negative number if an input parameter is invalid
108: -   -1 -  step < 0

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

125:   PetscFunctionBegin;
126:   PetscCall(PetscObjectGetComm((PetscObject)ls, &comm));
127:   fact       = 0.0;
128:   ls->nfeval = 0;
129:   ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
130:   if (!armP->work) {
131:     PetscCall(VecDuplicate(x, &armP->work));
132:     armP->x = x;
133:     PetscCall(PetscObjectReference((PetscObject)armP->x));
134:   } else if (x != armP->x) {
135:     PetscCall(VecDestroy(&armP->work));
136:     PetscCall(VecDuplicate(x, &armP->work));
137:     PetscCall(PetscObjectDereference((PetscObject)armP->x));
138:     armP->x = x;
139:     PetscCall(PetscObjectReference((PetscObject)armP->x));
140:   }

142:   PetscCall(TaoLineSearchMonitor(ls, 0, *f, 0.0));

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

171:   if (ls->reason != TAOLINESEARCH_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);

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

178:   if (!armP->memorySetup) {
179:     for (i = 0; i < armP->memorySize; i++) armP->memory[i] = armP->alpha * (*f);
180:     armP->current       = 0;
181:     armP->lastReference = armP->memory[0];
182:     armP->memorySetup   = PETSC_TRUE;
183:   }

185:   /* Calculate reference value (MAX) */
186:   ref = armP->memory[0];
187:   idx = 0;

189:   for (i = 1; i < armP->memorySize; i++) {
190:     if (armP->memory[i] > ref) {
191:       ref = armP->memory[i];
192:       idx = i;
193:     }
194:   }

196:   if (armP->referencePolicy == REFERENCE_AVE) {
197:     ref = 0;
198:     for (i = 0; i < armP->memorySize; i++) ref += armP->memory[i];
199:     ref = ref / armP->memorySize;
200:     ref = PetscMax(ref, armP->memory[armP->current]);
201:   } else if (armP->referencePolicy == REFERENCE_MEAN) {
202:     ref = PetscMin(ref, 0.5 * (armP->lastReference + armP->memory[armP->current]));
203:   }

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

207:   PetscCall(VecDuplicate(g, &g_old));
208:   PetscCall(VecCopy(g, g_old));

210:   ls->step = ls->initstep;
211:   while (ls->step >= owlqn_minstep && ls->nfeval < ls->max_funcs) {
212:     /* Calculate iterate */
213:     ++its;
214:     PetscCall(VecWAXPY(armP->work, ls->step, s, x));

216:     partgdx = 0.0;
217:     PetscCall(ProjWork_OWLQN(armP->work, x, g_old, &partgdx));
218:     PetscCall(MPIU_Allreduce(&partgdx, &gdx, 1, MPIU_REAL, MPIU_SUM, comm));

220:     /* Check the condition of gdx */
221:     if (PetscIsInfOrNanReal(gdx)) {
222:       PetscCall(PetscInfo(ls, "Initial Line Search step * g is Inf or Nan (%g)\n", (double)gdx));
223:       ls->reason = TAOLINESEARCH_FAILED_INFORNAN;
224:       PetscFunctionReturn(PETSC_SUCCESS);
225:     }
226:     if (gdx >= 0.0) {
227:       PetscCall(PetscInfo(ls, "Initial Line Search step is not descent direction (g's=%g)\n", (double)gdx));
228:       ls->reason = TAOLINESEARCH_FAILED_ASCENT;
229:       PetscFunctionReturn(PETSC_SUCCESS);
230:     }

232:     /* Calculate function at new iterate */
233:     PetscCall(TaoLineSearchComputeObjectiveAndGradient(ls, armP->work, f, g));
234:     g_computed = PETSC_TRUE;

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

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

240:     if (PetscIsInfOrNanReal(*f)) {
241:       ls->step *= armP->beta_inf;
242:     } else {
243:       /* Check descent condition */
244:       if (armP->nondescending && *f <= ref - ls->step * fact * ref) break;
245:       if (!armP->nondescending && *f <= ref + armP->sigma * gdx) break;
246:       ls->step *= armP->beta;
247:     }
248:   }
249:   PetscCall(VecDestroy(&g_old));

251:   /* Check termination */
252:   if (PetscIsInfOrNanReal(*f)) {
253:     PetscCall(PetscInfo(ls, "Function is inf or nan.\n"));
254:     ls->reason = TAOLINESEARCH_FAILED_BADPARAMETER;
255:   } else if (ls->step < owlqn_minstep) {
256:     PetscCall(PetscInfo(ls, "Step length is below tolerance.\n"));
257:     ls->reason = TAOLINESEARCH_HALTED_RTOL;
258:   } else if (ls->nfeval >= ls->max_funcs) {
259:     PetscCall(PetscInfo(ls, "Number of line search function evals (%" PetscInt_FMT ") > maximum allowed (%" PetscInt_FMT ")\n", ls->nfeval, ls->max_funcs));
260:     ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
261:   }
262:   if (ls->reason) PetscFunctionReturn(PETSC_SUCCESS);

264:   /* Successful termination, update memory */
265:   ls->reason          = TAOLINESEARCH_SUCCESS;
266:   armP->lastReference = ref;
267:   if (armP->replacementPolicy == REPLACE_FIFO) {
268:     armP->memory[armP->current++] = *f;
269:     if (armP->current >= armP->memorySize) armP->current = 0;
270:   } else {
271:     armP->current     = idx;
272:     armP->memory[idx] = *f;
273:   }

275:   /* Update iterate and compute gradient */
276:   PetscCall(VecCopy(armP->work, x));
277:   if (!g_computed) PetscCall(TaoLineSearchComputeGradient(ls, x, g));
278:   PetscCall(PetscInfo(ls, "%" PetscInt_FMT " function evals in line search, step = %10.4f\n", ls->nfeval, (double)ls->step));
279:   PetscFunctionReturn(PETSC_SUCCESS);
280: }

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

286:    Level: developer

288: .keywords: Tao, linesearch
289: M*/
290: PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_OWArmijo(TaoLineSearch ls)
291: {
292:   TaoLineSearch_OWARMIJO *armP;

294:   PetscFunctionBegin;
296:   PetscCall(PetscNew(&armP));

298:   armP->memory            = NULL;
299:   armP->alpha             = 1.0;
300:   armP->beta              = 0.25;
301:   armP->beta_inf          = 0.25;
302:   armP->sigma             = 1e-4;
303:   armP->memorySize        = 1;
304:   armP->referencePolicy   = REFERENCE_MAX;
305:   armP->replacementPolicy = REPLACE_MRU;
306:   armP->nondescending     = PETSC_FALSE;
307:   ls->data                = (void *)armP;
308:   ls->initstep            = 0.1;
309:   ls->ops->monitor        = NULL;
310:   ls->ops->setup          = NULL;
311:   ls->ops->reset          = NULL;
312:   ls->ops->apply          = TaoLineSearchApply_OWArmijo;
313:   ls->ops->view           = TaoLineSearchView_OWArmijo;
314:   ls->ops->destroy        = TaoLineSearchDestroy_OWArmijo;
315:   ls->ops->setfromoptions = TaoLineSearchSetFromOptions_OWArmijo;
316:   PetscFunctionReturn(PETSC_SUCCESS);
317: }