xref: /petsc/src/tao/linesearch/impls/owarmijo/owarmijo.c (revision bae903cb6e5fe998dd743ea97ba6d0e4e52494ec)
1 
2 #include <petsc/private/taolinesearchimpl.h>
3 #include <../src/tao/linesearch/impls/owarmijo/owarmijo.h>
4 
5 #define REPLACE_FIFO 1
6 #define REPLACE_MRU  2
7 
8 #define REFERENCE_MAX  1
9 #define REFERENCE_AVE  2
10 #define REFERENCE_MEAN 3
11 
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;
17 
18   PetscFunctionBegin;
19   PetscCall(VecGetOwnershipRange(w,&low,&high));
20   PetscCall(VecGetOwnershipRange(x,&low1,&high1));
21   PetscCall(VecGetOwnershipRange(gv,&low2,&high2));
22 
23   *gdx=0.0;
24   PetscCall(VecGetArray(w,&wptr));
25   PetscCall(VecGetArrayRead(x,&xptr));
26   PetscCall(VecGetArrayRead(gv,&gptr));
27 
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(0);
36 }
37 
38 static PetscErrorCode TaoLineSearchDestroy_OWArmijo(TaoLineSearch ls)
39 {
40   TaoLineSearch_OWARMIJO *armP = (TaoLineSearch_OWARMIJO *)ls->data;
41 
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(0);
48 }
49 
50 static PetscErrorCode TaoLineSearchSetFromOptions_OWArmijo(PetscOptionItems *PetscOptionsObject,TaoLineSearch ls)
51 {
52   TaoLineSearch_OWARMIJO *armP = (TaoLineSearch_OWARMIJO *)ls->data;
53 
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(0);
66 }
67 
68 static PetscErrorCode TaoLineSearchView_OWArmijo(TaoLineSearch ls, PetscViewer pv)
69 {
70   TaoLineSearch_OWARMIJO *armP = (TaoLineSearch_OWARMIJO *)ls->data;
71   PetscBool              isascii;
72 
73   PetscFunctionBegin;
74   PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
75   if (isascii) {
76     PetscCall(PetscViewerASCIIPrintf(pv,"  OWArmijo linesearch"));
77     if (armP->nondescending) {
78       PetscCall(PetscViewerASCIIPrintf(pv, " (nondescending)"));
79     }
80     PetscCall(PetscViewerASCIIPrintf(pv,": alpha=%g beta=%g ",(double)armP->alpha,(double)armP->beta));
81     PetscCall(PetscViewerASCIIPrintf(pv,"sigma=%g ",(double)armP->sigma));
82     PetscCall(PetscViewerASCIIPrintf(pv,"memsize=%" PetscInt_FMT "\n",armP->memorySize));
83   }
84   PetscFunctionReturn(0);
85 }
86 
87 /* @ TaoApply_OWArmijo - This routine performs a linesearch. It
88    backtracks until the (nonmonotone) OWArmijo conditions are satisfied.
89 
90    Input Parameters:
91 +  tao - TAO_SOLVER context
92 .  X - current iterate (on output X contains new iterate, X + step*S)
93 .  S - search direction
94 .  f - merit function evaluated at X
95 .  G - gradient of merit function evaluated at X
96 .  W - work vector
97 -  step - initial estimate of step length
98 
99    Output parameters:
100 +  f - merit function evaluated at new iterate, X + step*S
101 .  G - gradient of merit function evaluated at new iterate, X + step*S
102 .  X - new iterate
103 -  step - final step length
104 
105    Info is set to one of:
106 .   0 - the line search succeeds; the sufficient decrease
107    condition and the directional derivative condition hold
108 
109    negative number if an input parameter is invalid
110 -   -1 -  step < 0
111 
112    positive number > 1 if the line search otherwise terminates
113 +    1 -  Step is at the lower bound, stepmin.
114 @ */
115 static PetscErrorCode TaoLineSearchApply_OWArmijo(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s)
116 {
117   TaoLineSearch_OWARMIJO *armP = (TaoLineSearch_OWARMIJO *)ls->data;
118   PetscInt               i, its=0;
119   PetscReal              fact, ref, gdx;
120   PetscInt               idx;
121   PetscBool              g_computed=PETSC_FALSE; /* to prevent extra gradient computation */
122   Vec                    g_old;
123   PetscReal              owlqn_minstep=0.005;
124   PetscReal              partgdx;
125   MPI_Comm               comm;
126 
127   PetscFunctionBegin;
128   PetscCall(PetscObjectGetComm((PetscObject)ls,&comm));
129   fact = 0.0;
130   ls->nfeval=0;
131   ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
132   if (!armP->work) {
133     PetscCall(VecDuplicate(x,&armP->work));
134     armP->x = x;
135     PetscCall(PetscObjectReference((PetscObject)armP->x));
136   } else if (x != armP->x) {
137     PetscCall(VecDestroy(&armP->work));
138     PetscCall(VecDuplicate(x,&armP->work));
139     PetscCall(PetscObjectDereference((PetscObject)armP->x));
140     armP->x = x;
141     PetscCall(PetscObjectReference((PetscObject)armP->x));
142   }
143 
144   PetscCall(TaoLineSearchMonitor(ls, 0, *f, 0.0));
145 
146   /* Check linesearch parameters */
147   if (armP->alpha < 1) {
148     PetscCall(PetscInfo(ls,"OWArmijo line search error: alpha (%g) < 1\n", (double)armP->alpha));
149     ls->reason=TAOLINESEARCH_FAILED_BADPARAMETER;
150   } else if ((armP->beta <= 0) || (armP->beta >= 1)) {
151     PetscCall(PetscInfo(ls,"OWArmijo line search error: beta (%g) invalid\n", (double)armP->beta));
152     ls->reason=TAOLINESEARCH_FAILED_BADPARAMETER;
153   } else if ((armP->beta_inf <= 0) || (armP->beta_inf >= 1)) {
154     PetscCall(PetscInfo(ls,"OWArmijo line search error: beta_inf (%g) invalid\n", (double)armP->beta_inf));
155     ls->reason=TAOLINESEARCH_FAILED_BADPARAMETER;
156   } else if ((armP->sigma <= 0) || (armP->sigma >= 0.5)) {
157     PetscCall(PetscInfo(ls,"OWArmijo line search error: sigma (%g) invalid\n", (double)armP->sigma));
158     ls->reason=TAOLINESEARCH_FAILED_BADPARAMETER;
159   } else if (armP->memorySize < 1) {
160     PetscCall(PetscInfo(ls,"OWArmijo line search error: memory_size (%" PetscInt_FMT ") < 1\n", armP->memorySize));
161     ls->reason=TAOLINESEARCH_FAILED_BADPARAMETER;
162   }  else if ((armP->referencePolicy != REFERENCE_MAX) && (armP->referencePolicy != REFERENCE_AVE) && (armP->referencePolicy != REFERENCE_MEAN)) {
163     PetscCall(PetscInfo(ls,"OWArmijo line search error: reference_policy invalid\n"));
164     ls->reason=TAOLINESEARCH_FAILED_BADPARAMETER;
165   } else if ((armP->replacementPolicy != REPLACE_FIFO) && (armP->replacementPolicy != REPLACE_MRU)) {
166     PetscCall(PetscInfo(ls,"OWArmijo line search error: replacement_policy invalid\n"));
167     ls->reason=TAOLINESEARCH_FAILED_BADPARAMETER;
168   } else if (PetscIsInfOrNanReal(*f)) {
169     PetscCall(PetscInfo(ls,"OWArmijo line search error: initial function inf or nan\n"));
170     ls->reason=TAOLINESEARCH_FAILED_BADPARAMETER;
171   }
172 
173   if (ls->reason != TAOLINESEARCH_CONTINUE_ITERATING) PetscFunctionReturn(0);
174 
175   /* Check to see of the memory has been allocated.  If not, allocate
176      the historical array and populate it with the initial function
177      values. */
178   if (!armP->memory) {
179     PetscCall(PetscMalloc1(armP->memorySize, &armP->memory));
180   }
181 
182   if (!armP->memorySetup) {
183     for (i = 0; i < armP->memorySize; i++) {
184       armP->memory[i] = armP->alpha*(*f);
185     }
186     armP->current = 0;
187     armP->lastReference = armP->memory[0];
188     armP->memorySetup=PETSC_TRUE;
189   }
190 
191   /* Calculate reference value (MAX) */
192   ref = armP->memory[0];
193   idx = 0;
194 
195   for (i = 1; i < armP->memorySize; i++) {
196     if (armP->memory[i] > ref) {
197       ref = armP->memory[i];
198       idx = i;
199     }
200   }
201 
202   if (armP->referencePolicy == REFERENCE_AVE) {
203     ref = 0;
204     for (i = 0; i < armP->memorySize; i++) {
205       ref += armP->memory[i];
206     }
207     ref = ref / armP->memorySize;
208     ref = PetscMax(ref, armP->memory[armP->current]);
209   } else if (armP->referencePolicy == REFERENCE_MEAN) {
210     ref = PetscMin(ref, 0.5*(armP->lastReference + armP->memory[armP->current]));
211   }
212 
213   if (armP->nondescending) {
214     fact = armP->sigma;
215   }
216 
217   PetscCall(VecDuplicate(g,&g_old));
218   PetscCall(VecCopy(g,g_old));
219 
220   ls->step = ls->initstep;
221   while (ls->step >= owlqn_minstep && ls->nfeval < ls->max_funcs) {
222     /* Calculate iterate */
223     ++its;
224     PetscCall(VecWAXPY(armP->work,ls->step,s,x));
225 
226     partgdx=0.0;
227     PetscCall(ProjWork_OWLQN(armP->work,x,g_old,&partgdx));
228     PetscCall(MPIU_Allreduce(&partgdx,&gdx,1,MPIU_REAL,MPIU_SUM,comm));
229 
230     /* Check the condition of gdx */
231     if (PetscIsInfOrNanReal(gdx)) {
232       PetscCall(PetscInfo(ls,"Initial Line Search step * g is Inf or Nan (%g)\n",(double)gdx));
233       ls->reason=TAOLINESEARCH_FAILED_INFORNAN;
234       PetscFunctionReturn(0);
235     }
236     if (gdx >= 0.0) {
237       PetscCall(PetscInfo(ls,"Initial Line Search step is not descent direction (g's=%g)\n",(double)gdx));
238       ls->reason = TAOLINESEARCH_FAILED_ASCENT;
239       PetscFunctionReturn(0);
240     }
241 
242     /* Calculate function at new iterate */
243     PetscCall(TaoLineSearchComputeObjectiveAndGradient(ls,armP->work,f,g));
244     g_computed=PETSC_TRUE;
245 
246     PetscCall(TaoLineSearchMonitor(ls, its, *f, ls->step));
247 
248     if (ls->step == ls->initstep) {
249       ls->f_fullstep = *f;
250     }
251 
252     if (PetscIsInfOrNanReal(*f)) {
253       ls->step *= armP->beta_inf;
254     } else {
255       /* Check descent condition */
256       if (armP->nondescending && *f <= ref - ls->step*fact*ref) break;
257       if (!armP->nondescending && *f <= ref + armP->sigma * gdx) break;
258       ls->step *= armP->beta;
259     }
260   }
261   PetscCall(VecDestroy(&g_old));
262 
263   /* Check termination */
264   if (PetscIsInfOrNanReal(*f)) {
265     PetscCall(PetscInfo(ls, "Function is inf or nan.\n"));
266     ls->reason = TAOLINESEARCH_FAILED_BADPARAMETER;
267   } else if (ls->step < owlqn_minstep) {
268     PetscCall(PetscInfo(ls, "Step length is below tolerance.\n"));
269     ls->reason = TAOLINESEARCH_HALTED_RTOL;
270   } else if (ls->nfeval >= ls->max_funcs) {
271     PetscCall(PetscInfo(ls, "Number of line search function evals (%" PetscInt_FMT ") > maximum allowed (%" PetscInt_FMT ")\n",ls->nfeval, ls->max_funcs));
272     ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
273   }
274   if (ls->reason) PetscFunctionReturn(0);
275 
276   /* Successful termination, update memory */
277   ls->reason = TAOLINESEARCH_SUCCESS;
278   armP->lastReference = ref;
279   if (armP->replacementPolicy == REPLACE_FIFO) {
280     armP->memory[armP->current++] = *f;
281     if (armP->current >= armP->memorySize) {
282       armP->current = 0;
283     }
284   } else {
285     armP->current = idx;
286     armP->memory[idx] = *f;
287   }
288 
289   /* Update iterate and compute gradient */
290   PetscCall(VecCopy(armP->work,x));
291   if (!g_computed) {
292     PetscCall(TaoLineSearchComputeGradient(ls, x, g));
293   }
294   PetscCall(PetscInfo(ls, "%" PetscInt_FMT " function evals in line search, step = %10.4f\n",ls->nfeval, (double)ls->step));
295   PetscFunctionReturn(0);
296 }
297 
298 /*MC
299    TAOLINESEARCHOWARMIJO - Special line-search type for the Orthant-Wise Limited Quasi-Newton (TAOOWLQN) algorithm.
300    Should not be used with any other algorithm.
301 
302    Level: developer
303 
304 .keywords: Tao, linesearch
305 M*/
306 PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_OWArmijo(TaoLineSearch ls)
307 {
308   TaoLineSearch_OWARMIJO *armP;
309 
310   PetscFunctionBegin;
311   PetscValidHeaderSpecific(ls,TAOLINESEARCH_CLASSID,1);
312   PetscCall(PetscNewLog(ls,&armP));
313 
314   armP->memory = NULL;
315   armP->alpha = 1.0;
316   armP->beta = 0.25;
317   armP->beta_inf = 0.25;
318   armP->sigma = 1e-4;
319   armP->memorySize = 1;
320   armP->referencePolicy = REFERENCE_MAX;
321   armP->replacementPolicy = REPLACE_MRU;
322   armP->nondescending=PETSC_FALSE;
323   ls->data = (void*)armP;
324   ls->initstep = 0.1;
325   ls->ops->monitor = NULL;
326   ls->ops->setup = NULL;
327   ls->ops->reset = NULL;
328   ls->ops->apply = TaoLineSearchApply_OWArmijo;
329   ls->ops->view = TaoLineSearchView_OWArmijo;
330   ls->ops->destroy = TaoLineSearchDestroy_OWArmijo;
331   ls->ops->setfromoptions = TaoLineSearchSetFromOptions_OWArmijo;
332   PetscFunctionReturn(0);
333 }
334