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