Actual source code: neldermead.c

  1: #include <../src/tao/unconstrained/impls/neldermead/neldermead.h>
  2: #include <petscvec.h>

  4: /*------------------------------------------------------------*/
  5: static PetscErrorCode NelderMeadSort(TAO_NelderMead *nm)
  6: {
  7:   PetscReal *values  = nm->f_values;
  8:   PetscInt  *indices = nm->indices;
  9:   PetscInt   dim     = nm->N + 1;
 10:   PetscInt   i, j, index;
 11:   PetscReal  val;

 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:   return 0;
 20: }

 22: /*------------------------------------------------------------*/
 23: static PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f)
 24: {
 25:   /*  Add new vector's fraction of average */
 26:   VecAXPY(nm->Xbar, nm->oneOverN, Xmu);
 27:   VecCopy(Xmu, nm->simplex[index]);
 28:   nm->f_values[index] = f;

 30:   NelderMeadSort(nm);

 32:   /*  Subtract last vector from average */
 33:   VecAXPY(nm->Xbar, -nm->oneOverN, nm->simplex[nm->indices[nm->N]]);
 34:   return 0;
 35: }

 37: /* ---------------------------------------------------------- */
 38: static PetscErrorCode TaoSetUp_NM(Tao tao)
 39: {
 40:   TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
 41:   PetscInt        n;

 43:   VecGetSize(tao->solution, &n);
 44:   nm->N        = n;
 45:   nm->oneOverN = 1.0 / n;
 46:   VecDuplicateVecs(tao->solution, nm->N + 1, &nm->simplex);
 47:   PetscMalloc1(nm->N + 1, &nm->f_values);
 48:   PetscMalloc1(nm->N + 1, &nm->indices);
 49:   VecDuplicate(tao->solution, &nm->Xbar);
 50:   VecDuplicate(tao->solution, &nm->Xmur);
 51:   VecDuplicate(tao->solution, &nm->Xmue);
 52:   VecDuplicate(tao->solution, &nm->Xmuc);

 54:   tao->gradient = NULL;
 55:   tao->step     = 0;
 56:   return 0;
 57: }

 59: /* ---------------------------------------------------------- */
 60: static PetscErrorCode TaoDestroy_NM(Tao tao)
 61: {
 62:   TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;

 64:   if (tao->setupcalled) {
 65:     VecDestroyVecs(nm->N + 1, &nm->simplex);
 66:     VecDestroy(&nm->Xmuc);
 67:     VecDestroy(&nm->Xmue);
 68:     VecDestroy(&nm->Xmur);
 69:     VecDestroy(&nm->Xbar);
 70:   }
 71:   PetscFree(nm->indices);
 72:   PetscFree(nm->f_values);
 73:   PetscFree(tao->data);
 74:   return 0;
 75: }

 77: /*------------------------------------------------------------*/
 78: static PetscErrorCode TaoSetFromOptions_NM(Tao tao, PetscOptionItems *PetscOptionsObject)
 79: {
 80:   TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;

 82:   PetscOptionsHeadBegin(PetscOptionsObject, "Nelder-Mead options");
 83:   PetscOptionsDeprecated("-tao_nm_lamda", "-tao_nm_lambda", "3.18.4", NULL);
 84:   PetscOptionsReal("-tao_nm_lambda", "initial step length", "", nm->lambda, &nm->lambda, NULL);
 85:   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:   return 0;
 91: }

 93: /*------------------------------------------------------------*/
 94: static PetscErrorCode TaoView_NM(Tao tao, PetscViewer viewer)
 95: {
 96:   TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
 97:   PetscBool       isascii;

 99:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
100:   if (isascii) {
101:     PetscViewerASCIIPushTab(viewer);
102:     PetscViewerASCIIPrintf(viewer, "expansions: %" PetscInt_FMT "\n", nm->nexpand);
103:     PetscViewerASCIIPrintf(viewer, "reflections: %" PetscInt_FMT "\n", nm->nreflect);
104:     PetscViewerASCIIPrintf(viewer, "inside contractions: %" PetscInt_FMT "\n", nm->nincontract);
105:     PetscViewerASCIIPrintf(viewer, "outside contractionss: %" PetscInt_FMT "\n", nm->noutcontract);
106:     PetscViewerASCIIPrintf(viewer, "Shrink steps: %" PetscInt_FMT "\n", nm->nshrink);
107:     PetscViewerASCIIPopTab(viewer);
108:   }
109:   return 0;
110: }

112: /*------------------------------------------------------------*/
113: static PetscErrorCode TaoSolve_NM(Tao tao)
114: {
115:   TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
116:   PetscReal      *x;
117:   PetscInt        i;
118:   Vec             Xmur = nm->Xmur, Xmue = nm->Xmue, Xmuc = nm->Xmuc, Xbar = nm->Xbar;
119:   PetscReal       fr, fe, fc;
120:   PetscInt        shrink;
121:   PetscInt        low, high;

123:   nm->nshrink      = 0;
124:   nm->nreflect     = 0;
125:   nm->nincontract  = 0;
126:   nm->noutcontract = 0;
127:   nm->nexpand      = 0;

129:   if (tao->XL || tao->XU || tao->ops->computebounds) PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n");

131:   VecCopy(tao->solution, nm->simplex[0]);
132:   TaoComputeObjective(tao, nm->simplex[0], &nm->f_values[0]);
133:   nm->indices[0] = 0;
134:   for (i = 1; i < nm->N + 1; i++) {
135:     VecCopy(tao->solution, nm->simplex[i]);
136:     VecGetOwnershipRange(nm->simplex[i], &low, &high);
137:     if (i - 1 >= low && i - 1 < high) {
138:       VecGetArray(nm->simplex[i], &x);
139:       x[i - 1 - low] += nm->lambda;
140:       VecRestoreArray(nm->simplex[i], &x);
141:     }

143:     TaoComputeObjective(tao, nm->simplex[i], &nm->f_values[i]);
144:     nm->indices[i] = i;
145:   }

147:   /*  Xbar  = (Sum of all simplex vectors - worst vector)/N */
148:   NelderMeadSort(nm);
149:   VecSet(Xbar, 0.0);
150:   for (i = 0; i < nm->N; i++) VecAXPY(Xbar, 1.0, nm->simplex[nm->indices[i]]);
151:   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:     VecCopy(nm->simplex[nm->indices[0]], tao->solution);
159:     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:     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;

164:     /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */
165:     VecAXPBYPCZ(Xmur, 1 + nm->mu_r, -nm->mu_r, 0, Xbar, nm->simplex[nm->indices[nm->N]]);
166:     TaoComputeObjective(tao, Xmur, &fr);

168:     if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N - 1]]) {
169:       /*  reflect */
170:       nm->nreflect++;
171:       PetscInfo(0, "Reflect\n");
172:       NelderMeadReplace(nm, nm->indices[nm->N], Xmur, fr);
173:     } else if (fr < nm->f_values[nm->indices[0]]) {
174:       /*  expand */
175:       nm->nexpand++;
176:       PetscInfo(0, "Expand\n");
177:       VecAXPBYPCZ(Xmue, 1 + nm->mu_e, -nm->mu_e, 0, Xbar, nm->simplex[nm->indices[nm->N]]);
178:       TaoComputeObjective(tao, Xmue, &fe);
179:       if (fe < fr) {
180:         NelderMeadReplace(nm, nm->indices[nm->N], Xmue, fe);
181:       } else {
182:         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:       PetscInfo(0, "Outside Contraction\n");
188:       VecAXPBYPCZ(Xmuc, 1 + nm->mu_oc, -nm->mu_oc, 0, Xbar, nm->simplex[nm->indices[nm->N]]);

190:       TaoComputeObjective(tao, Xmuc, &fc);
191:       if (fc <= fr) {
192:         NelderMeadReplace(nm, nm->indices[nm->N], Xmuc, fc);
193:       } else shrink = 1;
194:     } else {
195:       /* inside contraction */
196:       nm->nincontract++;
197:       PetscInfo(0, "Inside Contraction\n");
198:       VecAXPBYPCZ(Xmuc, 1 + nm->mu_ic, -nm->mu_ic, 0, Xbar, nm->simplex[nm->indices[nm->N]]);
199:       TaoComputeObjective(tao, Xmuc, &fc);
200:       if (fc < nm->f_values[nm->indices[nm->N]]) {
201:         NelderMeadReplace(nm, nm->indices[nm->N], Xmuc, fc);
202:       } else shrink = 1;
203:     }

205:     if (shrink) {
206:       nm->nshrink++;
207:       PetscInfo(0, "Shrink\n");

209:       for (i = 1; i < nm->N + 1; i++) {
210:         VecAXPBY(nm->simplex[nm->indices[i]], 1.5, -0.5, nm->simplex[nm->indices[0]]);
211:         TaoComputeObjective(tao, nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]]);
212:       }
213:       VecAXPBY(Xbar, 1.5 * nm->oneOverN, -0.5, nm->simplex[nm->indices[0]]);

215:       /*  Add last vector's fraction of average */
216:       VecAXPY(Xbar, nm->oneOverN, nm->simplex[nm->indices[nm->N]]);
217:       NelderMeadSort(nm);
218:       /*  Subtract new last vector from average */
219:       VecAXPY(Xbar, -nm->oneOverN, nm->simplex[nm->indices[nm->N]]);
220:     }
221:   }
222:   return 0;
223: }

225: /* ---------------------------------------------------------- */
226: /*MC
227:  TAONM - Nelder-Mead solver for derivative free, unconstrained minimization

229:  Options Database Keys:
230: + -tao_nm_lambda - initial step length
231: - -tao_nm_mu - expansion/contraction factor

233:  Level: beginner
234: M*/

236: PETSC_EXTERN PetscErrorCode TaoCreate_NM(Tao tao)
237: {
238:   TAO_NelderMead *nm;

240:   PetscNew(&nm);
241:   tao->data = (void *)nm;

243:   tao->ops->setup          = TaoSetUp_NM;
244:   tao->ops->solve          = TaoSolve_NM;
245:   tao->ops->view           = TaoView_NM;
246:   tao->ops->setfromoptions = TaoSetFromOptions_NM;
247:   tao->ops->destroy        = TaoDestroy_NM;

249:   /* Override default settings (unless already changed) */
250:   if (!tao->max_it_changed) tao->max_it = 2000;
251:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;

253:   nm->simplex = NULL;
254:   nm->lambda  = 1;

256:   nm->mu_ic = -0.5;
257:   nm->mu_oc = 0.5;
258:   nm->mu_r  = 1.0;
259:   nm->mu_e  = 2.0;

261:   return 0;
262: }