xref: /petsc/src/tao/unconstrained/impls/neldermead/neldermead.c (revision a7fac7c2b47dbcd8ece0cc2dfdfe0e63be1bb7b5) 
1 #include <../src/tao/unconstrained/impls/neldermead/neldermead.h>
2 #include <petscvec.h>
3 
4 static PetscErrorCode NelderMeadSort(TAO_NelderMead *nm);
5 static PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f);
6 /* ---------------------------------------------------------- */
7 #undef __FUNCT__
8 #define __FUNCT__ "TaoSetUp_NM"
9 static PetscErrorCode TaoSetUp_NM(Tao tao)
10 {
11   PetscErrorCode ierr;
12   TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
13   PetscInt       n;
14 
15   PetscFunctionBegin;
16   ierr = VecGetSize(tao->solution,&n);CHKERRQ(ierr);
17   nm->N = n;
18   nm->oneOverN = 1.0/n;
19   ierr = VecDuplicateVecs(tao->solution,nm->N+1,&nm->simplex);CHKERRQ(ierr);
20   ierr = PetscMalloc1(nm->N+1,&nm->f_values);CHKERRQ(ierr);
21   ierr = PetscMalloc1(nm->N+1,&nm->indices);CHKERRQ(ierr);
22   ierr = VecDuplicate(tao->solution,&nm->Xbar);CHKERRQ(ierr);
23   ierr = VecDuplicate(tao->solution,&nm->Xmur);CHKERRQ(ierr);
24   ierr = VecDuplicate(tao->solution,&nm->Xmue);CHKERRQ(ierr);
25   ierr = VecDuplicate(tao->solution,&nm->Xmuc);CHKERRQ(ierr);
26 
27   tao->gradient=0;
28   tao->step=0;
29   PetscFunctionReturn(0);
30 }
31 
32 /* ---------------------------------------------------------- */
33 #undef __FUNCT__
34 #define __FUNCT__ "TaoDestroy_NM"
35 PetscErrorCode TaoDestroy_NM(Tao tao)
36 {
37   TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
38   PetscErrorCode ierr;
39 
40   PetscFunctionBegin;
41   if (tao->setupcalled) {
42     ierr = VecDestroyVecs(nm->N+1,&nm->simplex);CHKERRQ(ierr);
43     ierr = VecDestroy(&nm->Xmuc);CHKERRQ(ierr);
44     ierr = VecDestroy(&nm->Xmue);CHKERRQ(ierr);
45     ierr = VecDestroy(&nm->Xmur);CHKERRQ(ierr);
46     ierr = VecDestroy(&nm->Xbar);CHKERRQ(ierr);
47   }
48   ierr = PetscFree(nm->indices);CHKERRQ(ierr);
49   ierr = PetscFree(nm->f_values);CHKERRQ(ierr);
50   ierr = PetscFree(tao->data);CHKERRQ(ierr);
51   tao->data = 0;
52   PetscFunctionReturn(0);
53 }
54 
55 /*------------------------------------------------------------*/
56 #undef __FUNCT__
57 #define __FUNCT__ "TaoSetFromOptions_NM"
58 PetscErrorCode TaoSetFromOptions_NM(PetscOptions *PetscOptionsObject,Tao tao)
59 {
60   TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
61   PetscErrorCode ierr;
62 
63   PetscFunctionBegin;
64   ierr = PetscOptionsHead(PetscOptionsObject,"Nelder-Mead options");CHKERRQ(ierr);
65   ierr = PetscOptionsReal("-tao_nm_lamda","initial step length","",nm->lamda,&nm->lamda,NULL); CHKERRQ(ierr);
66   ierr = PetscOptionsReal("-tao_nm_mu","mu","",nm->mu_oc,&nm->mu_oc,NULL);CHKERRQ(ierr);
67   nm->mu_ic = -nm->mu_oc;
68   nm->mu_r = nm->mu_oc*2.0;
69   nm->mu_e = nm->mu_oc*4.0;
70   ierr = PetscOptionsTail();CHKERRQ(ierr);
71   PetscFunctionReturn(0);
72 }
73 
74 /*------------------------------------------------------------*/
75 #undef __FUNCT__
76 #define __FUNCT__ "TaoView_NM"
77 PetscErrorCode TaoView_NM(Tao tao,PetscViewer viewer)
78 {
79   TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
80   PetscBool      isascii;
81   PetscErrorCode ierr;
82 
83   PetscFunctionBegin;
84   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
85   if (isascii) {
86     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
87     ierr = PetscViewerASCIIPrintf(viewer,"expansions: %D\n",nm->nexpand);CHKERRQ(ierr);
88     ierr = PetscViewerASCIIPrintf(viewer,"reflections: %D\n",nm->nreflect);CHKERRQ(ierr);
89     ierr = PetscViewerASCIIPrintf(viewer,"inside contractions: %D\n",nm->nincontract);CHKERRQ(ierr);
90     ierr = PetscViewerASCIIPrintf(viewer,"outside contractionss: %D\n",nm->noutcontract);CHKERRQ(ierr);
91     ierr = PetscViewerASCIIPrintf(viewer,"Shrink steps: %D\n",nm->nshrink);CHKERRQ(ierr);
92     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
93   }
94   PetscFunctionReturn(0);
95 }
96 
97 /*------------------------------------------------------------*/
98 #undef __FUNCT__
99 #define __FUNCT__ "TaoSolve_NM"
100 PetscErrorCode TaoSolve_NM(Tao tao)
101 {
102   PetscErrorCode     ierr;
103   TAO_NelderMead     *nm = (TAO_NelderMead*)tao->data;
104   TaoConvergedReason reason;
105   PetscReal          *x;
106   PetscInt           i;
107   Vec                Xmur=nm->Xmur, Xmue=nm->Xmue, Xmuc=nm->Xmuc, Xbar=nm->Xbar;
108   PetscReal          fr,fe,fc;
109   PetscInt           shrink;
110   PetscInt           low,high;
111 
112   PetscFunctionBegin;
113   nm->nshrink =      0;
114   nm->nreflect =     0;
115   nm->nincontract =  0;
116   nm->noutcontract = 0;
117   nm->nexpand =      0;
118 
119   if (tao->XL || tao->XU || tao->ops->computebounds) {
120     ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n");CHKERRQ(ierr);
121   }
122 
123   ierr = VecCopy(tao->solution,nm->simplex[0]);CHKERRQ(ierr);
124   ierr = TaoComputeObjective(tao,nm->simplex[0],&nm->f_values[0]);CHKERRQ(ierr);
125   nm->indices[0]=0;
126   for (i=1;i<nm->N+1;i++){
127     ierr = VecCopy(tao->solution,nm->simplex[i]);CHKERRQ(ierr);
128     ierr = VecGetOwnershipRange(nm->simplex[i],&low,&high);CHKERRQ(ierr);
129     if (i-1 >= low && i-1 < high) {
130       ierr = VecGetArray(nm->simplex[i],&x);CHKERRQ(ierr);
131       x[i-1-low] += nm->lamda;
132       ierr = VecRestoreArray(nm->simplex[i],&x);CHKERRQ(ierr);
133     }
134 
135     ierr = TaoComputeObjective(tao,nm->simplex[i],&nm->f_values[i]);CHKERRQ(ierr);
136     nm->indices[i] = i;
137   }
138 
139   /*  Xbar  = (Sum of all simplex vectors - worst vector)/N */
140   ierr = NelderMeadSort(nm);CHKERRQ(ierr);
141   ierr = VecSet(Xbar,0.0);CHKERRQ(ierr);
142   for (i=0;i<nm->N;i++) {
143     ierr = VecAXPY(Xbar,1.0,nm->simplex[nm->indices[i]]);CHKERRQ(ierr);
144   }
145   ierr = VecScale(Xbar,nm->oneOverN);CHKERRQ(ierr);
146   reason = TAO_CONTINUE_ITERATING;
147   while (1) {
148     shrink = 0;
149     ierr = VecCopy(nm->simplex[nm->indices[0]],tao->solution);CHKERRQ(ierr);
150     ierr = TaoMonitor(tao,tao->niter++,nm->f_values[nm->indices[0]],nm->f_values[nm->indices[nm->N]]-nm->f_values[nm->indices[0]],0.0,1.0,&reason);CHKERRQ(ierr);
151     if (reason != TAO_CONTINUE_ITERATING) break;
152 
153     /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */
154     ierr = VecAXPBYPCZ(Xmur,1+nm->mu_r,-nm->mu_r,0,Xbar,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr);
155     ierr = TaoComputeObjective(tao,Xmur,&fr);CHKERRQ(ierr);
156 
157     if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N-1]]) {
158       /*  reflect */
159       nm->nreflect++;
160       ierr = PetscInfo(0,"Reflect\n");CHKERRQ(ierr);
161       ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);CHKERRQ(ierr);
162     } else if (fr < nm->f_values[nm->indices[0]]) {
163       /*  expand */
164       nm->nexpand++;
165       ierr = PetscInfo(0,"Expand\n");CHKERRQ(ierr);
166       ierr = VecAXPBYPCZ(Xmue,1+nm->mu_e,-nm->mu_e,0,Xbar,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr);
167       ierr = TaoComputeObjective(tao,Xmue,&fe);CHKERRQ(ierr);
168       if (fe < fr) {
169         ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe);CHKERRQ(ierr);
170       } else {
171         ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);CHKERRQ(ierr);
172       }
173     } else if (nm->f_values[nm->indices[nm->N-1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) {
174       /* outside contraction */
175       nm->noutcontract++;
176       ierr = PetscInfo(0,"Outside Contraction\n");CHKERRQ(ierr);
177       ierr = VecAXPBYPCZ(Xmuc,1+nm->mu_oc,-nm->mu_oc,0,Xbar,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr);
178 
179       ierr = TaoComputeObjective(tao,Xmuc,&fc);CHKERRQ(ierr);
180       if (fc <= fr) {
181         ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);CHKERRQ(ierr);
182       } else shrink=1;
183     } else {
184       /* inside contraction */
185       nm->nincontract++;
186       ierr = PetscInfo(0,"Inside Contraction\n");CHKERRQ(ierr);
187       ierr = VecAXPBYPCZ(Xmuc,1+nm->mu_ic,-nm->mu_ic,0,Xbar,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr);
188       ierr = TaoComputeObjective(tao,Xmuc,&fc);CHKERRQ(ierr);
189       if (fc < nm->f_values[nm->indices[nm->N]]) {
190         ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);CHKERRQ(ierr);
191       } else shrink = 1;
192     }
193 
194     if (shrink) {
195       nm->nshrink++;
196       ierr = PetscInfo(0,"Shrink\n");CHKERRQ(ierr);
197 
198       for (i=1;i<nm->N+1;i++) {
199         ierr = VecAXPBY(nm->simplex[nm->indices[i]],1.5,-0.5,nm->simplex[nm->indices[0]]);CHKERRQ(ierr);
200         ierr = TaoComputeObjective(tao,nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]]);CHKERRQ(ierr);
201       }
202       ierr = VecAXPBY(Xbar,1.5*nm->oneOverN,-0.5,nm->simplex[nm->indices[0]]);CHKERRQ(ierr);
203 
204       /*  Add last vector's fraction of average */
205       ierr = VecAXPY(Xbar,nm->oneOverN,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr);
206       ierr = NelderMeadSort(nm);CHKERRQ(ierr);
207       /*  Subtract new last vector from average */
208       ierr = VecAXPY(Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr);
209     }
210   }
211   PetscFunctionReturn(0);
212 }
213 
214 /* ---------------------------------------------------------- */
215 /*MC
216  TAONM - Nelder-Mead solver for derivative free, unconstrained minimization
217 
218  Options Database Keys:
219 + -tao_nm_lamda - initial step length
220 . -tao_nm_mu - expansion/contraction factor
221 
222  Level: beginner
223 M*/
224 
225 #undef __FUNCT__
226 #define __FUNCT__ "TaoCreate_NM"
227 PETSC_EXTERN PetscErrorCode TaoCreate_NM(Tao tao)
228 {
229   TAO_NelderMead *nm;
230   PetscErrorCode ierr;
231 
232   PetscFunctionBegin;
233   ierr = PetscNewLog(tao,&nm);CHKERRQ(ierr);
234   tao->data = (void*)nm;
235 
236   tao->ops->setup = TaoSetUp_NM;
237   tao->ops->solve = TaoSolve_NM;
238   tao->ops->view = TaoView_NM;
239   tao->ops->setfromoptions = TaoSetFromOptions_NM;
240   tao->ops->destroy = TaoDestroy_NM;
241 
242   /* Override default settings (unless already changed) */
243   if (!tao->max_it_changed) tao->max_it = 2000;
244   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
245 #if defined(PETSC_USE_REAL_SINGLE)
246   if (!tao->fatol_changed) tao->fatol = 1.0e-5;
247   if (!tao->frtol_changed) tao->frtol = 1.0e-5;
248 #else
249   if (!tao->fatol_changed) tao->fatol = 1.0e-8;
250   if (!tao->frtol_changed) tao->frtol = 1.0e-8;
251 #endif
252 
253   nm->simplex = 0;
254   nm->lamda = 1;
255 
256   nm->mu_ic = -0.5;
257   nm->mu_oc = 0.5;
258   nm->mu_r = 1.0;
259   nm->mu_e = 2.0;
260 
261   PetscFunctionReturn(0);
262 }
263 
264 
265 /*------------------------------------------------------------*/
266 #undef __FUNCT__
267 #define __FUNCT__ "NelderMeadSort"
268 PetscErrorCode NelderMeadSort(TAO_NelderMead *nm)
269 {
270   PetscReal *values = nm->f_values;
271   PetscInt  *indices = nm->indices;
272   PetscInt  dim = nm->N+1;
273   PetscInt  i,j,index;
274   PetscReal val;
275 
276   PetscFunctionBegin;
277   for (i=1;i<dim;i++) {
278     index = indices[i];
279     val = values[index];
280     for (j=i-1; j>=0 && values[indices[j]] > val; j--) {
281       indices[j+1] = indices[j];
282     }
283     indices[j+1] = index;
284   }
285   PetscFunctionReturn(0);
286 }
287 
288 
289 /*------------------------------------------------------------*/
290 #undef __FUNCT__
291 #define __FUNCT__ "NelderMeadReplace"
292 PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f)
293 {
294   PetscErrorCode ierr;
295 
296   PetscFunctionBegin;
297   /*  Add new vector's fraction of average */
298   ierr = VecAXPY(nm->Xbar,nm->oneOverN,Xmu);CHKERRQ(ierr);
299   ierr = VecCopy(Xmu,nm->simplex[index]);CHKERRQ(ierr);
300   nm->f_values[index] = f;
301 
302   ierr = NelderMeadSort(nm);CHKERRQ(ierr);
303 
304   /*  Subtract last vector from average */
305   ierr = VecAXPY(nm->Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr);
306   PetscFunctionReturn(0);
307 }
308