plugins/nlsolve/gnm-nlsolve.c

   1 #include <gnumeric-config.h>
   2 #include "gnumeric.h"
   3 #include <tools/gnm-solver.h>
   4 #include "cell.h"
   5 #include "sheet.h"
   6 #include "value.h"
   7 #include "regression.h"
   8 #include "rangefunc.h"
   9 #include "workbook.h"
  10 #include "mathfunc.h"
  11 #include <gsf/gsf-impl-utils.h>
  12 #include <glib/gi18n-lib.h>
  13 #include <string.h>
  14
  15 /*
  16  * This is based on the algorithm from "An Automatic Method for finding
  17  * the Greatest or Least Value of a Function" by H. H. Rosenbrock
  18  * published in _The Computer Journal_ (1960) 3(3): 175-184.
  19  *
  20  * It is thus 50+ years old.  You would think that advances in computer
  21  * science would have produced improvements that would run circles
  22  * around this, but that is not obviously true.
  23  *
  24  * There are a couple of attrictive features of the Rosenbrock method:
  25  * 1. It's monotonic.  Unlike Newton-style methods it cannot suddenly
  26  *    warp far away.
  27  * 2. We don't need the Hessian.
  28  *
  29  * Note, that in order to speed convergence we occasionally perform
  30  * a tentative Newton iteration step.  (It's tentative because we will
  31  * discard it if it doesn't lead to an immediate improvement.)
  32  */
  33
  34 #define PRIVATE_KEY "::nlsolve::"
  35
  36 /*
  37  * Note: the solver code assumes the problem is a minimization problem.
  38  * When used for a maximization problem, we flip the objective function
  39  * sign.
  40  */
  41
  42
  43 typedef struct {
  44         /* The solver object in two forms.  */
  45         GnmSolver *sol;
  46         GnmIterSolver *isol;
  47
  48         /* Number of vars.  */
  49         int n;
  50
  51         /* Rosenbrock state */
  52         gnm_float **xi;
  53         int smallsteps;
  54         int tentative;
  55         gnm_float *tentative_xk, tentative_yk;
  56
  57         // Newton state
  58         // (nothing right now)
  59
  60         /* Parameters: */
  61         gboolean debug;
  62         gnm_float min_factor;
  63 } GnmNlsolve;
  64
  65 static gboolean
  66 check_program (GnmSolver *sol, GError **err)
  67 {
  68         unsigned ui;
  69         const GnmSolverParameters *params = sol->params;
  70         GSList *l;
  71
  72         for (l = params->constraints; l; l = l->next) {
  73                 GnmSolverConstraint *c  = l->data;
  74                 switch (c->type) {
  75                 case GNM_SOLVER_EQ:
  76                         /*
  77                          * This catches also equalities where the sides are not
  78                          * input variables.
  79                          */
  80                         goto no_equal;
  81                 default:
  82                         break;
  83                 }
  84         }
  85
  86         for (ui = 0; ui < sol->input_cells->len; ui++) {
  87                 if (sol->discrete[ui])
  88                         goto no_discrete;
  89
  90                 /*
  91                  * This also catches using two inequality constraints used
  92                  * to emulate equality.
  93                  */
  94                 if (sol->min[ui] == sol->max[ui])
  95                         goto no_equal;
  96         }
  97
  98         return TRUE;
  99
 100 no_discrete:
 101         g_set_error (err,
 102                      go_error_invalid (),
 103                      0,
 104                      _("This solver does not handle discrete variables."));
 105         return FALSE;
 106
 107 no_equal:
 108         g_set_error (err,
 109                      go_error_invalid (),
 110                      0,
 111                      _("This solver does not handle equality constraints."));
 112         return FALSE;
 113 }
 114
 115 static void
 116 print_vector (const char *name, const gnm_float *v, int n)
 117 {
 118         int i;
 119
 120         if (name)
 121                 g_printerr ("%s:\n", name);
 122         for (i = 0; i < n; i++)
 123                 g_printerr ("%15.8" GNM_FORMAT_f " ", v[i]);
 124         g_printerr ("\n");
 125 }
 126
 127 #if 0
 128 static void
 129 set_value (GnmNlsolve *nl, int i, gnm_float x)
 130 {
 131         gnm_solver_set_var (nl->sol, i, x);
 132 }
 133 #endif
 134
 135 static void
 136 set_vector (GnmNlsolve *nl, const gnm_float *xs)
 137 {
 138         gnm_solver_set_vars (nl->sol, xs);
 139 }
 140
 141 /* Get the target value as-if we were minimizing.  */
 142 static gnm_float
 143 get_value (GnmNlsolve *nl)
 144 {
 145         /* nl->sol has been taught to flip sign if needed.  */
 146         return gnm_solver_get_target_value (nl->sol);
 147 }
 148
 149 static void
 150 free_matrix (gnm_float **m, int n)
 151 {
 152         int i;
 153         for (i = 0; i < n; i++)
 154                 g_free (m[i]);
 155         g_free (m);
 156 }
 157
 158 static void
 159 set_solution (GnmNlsolve *nl)
 160 {
 161         /* nl->isol has been taught to flip sign if needed.  */
 162         gnm_iter_solver_set_solution (nl->isol);
 163 }
 164
 165 static gboolean
 166 gnm_nlsolve_prepare (GnmSolver *sol, WorkbookControl *wbc, GError **err,
 167                      GnmNlsolve *nl)
 168 {
 169         gboolean ok;
 170
 171         g_return_val_if_fail (sol->status == GNM_SOLVER_STATUS_READY, FALSE);
 172
 173         gnm_solver_set_status (sol, GNM_SOLVER_STATUS_PREPARING);
 174
 175         ok = check_program (sol, err);
 176         if (ok)
 177                 ok = gnm_iter_solver_get_initial_solution (nl->isol, err);
 178
 179         if (ok) {
 180                 gnm_solver_set_status (sol, GNM_SOLVER_STATUS_PREPARED);
 181         } else {
 182                 gnm_solver_set_status (sol, GNM_SOLVER_STATUS_ERROR);
 183         }
 184
 185         return ok;
 186 }
 187
 188 static gnm_float *
 189 compute_gradient (GnmNlsolve *nl, const gnm_float *xs)
 190 {
 191         return gnm_solver_compute_gradient (nl->sol, xs);
 192 }
 193
 194 static gboolean
 195 newton_improve (GnmNlsolve *nl, gnm_float *xs)
 196 {
 197         GnmSolver *sol = nl->sol;
 198         GnmIterSolver *isol = nl->isol;
 199         const int n = nl->n;
 200         int i;
 201         gnm_float *g, *d, *xs2;
 202         GnmMatrix *H;
 203         gboolean ok;
 204
 205         xs2 = g_new (gnm_float, n);
 206         g = compute_gradient (nl, xs);
 207         H = gnm_solver_compute_hessian (sol, xs);
 208         d = g_new (gnm_float, n);
 209         ok = (gnm_linear_solve_posdef (H, g, d) == GO_REG_ok);
 210         if (ok) {
 211                 for (i = 0; i < n; i++)
 212                         d[i] = 0 - d[i];
 213         }
 214
 215         if (nl->debug) {
 216                 int i;
 217                 g_printerr ("Hessian:\n");
 218                 for (i = 0; i < n; i++)
 219                         print_vector (NULL, H->data[i], n);
 220                 print_vector ("g", g, n);
 221                 if (ok)
 222                         print_vector ("d", d, n);
 223                 else
 224                         g_printerr ("Failed to solve Newton step.\n");
 225         }
 226
 227         if (ok) {
 228                 gnm_float y2, f;
 229
 230                 for (i = 0; i < n; i++)
 231                         xs2[i] = xs[i] + d[i];
 232                 set_vector (nl, xs2);
 233                 y2 = get_value (nl);
 234
 235                 ok = FALSE;
 236                 if (y2 < isol->yk && gnm_solver_check_constraints (sol)) {
 237                         if (nl->debug)
 238                                 g_printerr ("Accepting newton step\n");
 239                         memcpy (isol->xk, xs2, n * sizeof (gnm_float));
 240                         isol->yk = y2;
 241                         set_solution (nl);
 242                         ok = TRUE;
 243                 } else {
 244                         if (nl->debug)
 245                                 g_printerr ("Full newton step would go to %g\n", y2);
 246
 247                         f = gnm_solver_line_search
 248                                 (sol, xs, d, FALSE, 0.75, 1, 0.01, &y2);
 249
 250                         if (f > 0 && f < 1 && y2 < isol->yk) {
 251                                 if (nl->debug)
 252                                         g_printerr ("Accepting reduced newton step with f=%g\n", f);
 253                                 for (i = 0; i < n; i++)
 254                                         isol->xk[i] = xs[i] + f * d[i];
 255                                 isol->yk = y2;
 256                                 set_solution (nl);
 257                                 ok = TRUE;
 258                         }
 259                 }
 260         }
 261
 262         g_free (d);
 263         g_free (g);
 264         gnm_matrix_free (H);
 265         g_free (xs2);
 266
 267         return ok;
 268 }
 269
 270 static void
 271 nlsolve_init (GnmNlsolve *nl)
 272 {
 273         const int n = nl->n;
 274         int i, j;
 275
 276         nl->xi = g_new (gnm_float *, n);
 277         for (i = 0; i < n; i++) {
 278                 nl->xi[i] = g_new (gnm_float, n);
 279                 for (j = 0; j < n; j++)
 280                         nl->xi[i][j] = (i == j);
 281         }
 282
 283         nl->smallsteps = 0;
 284
 285         nl->tentative = 0;
 286         nl->tentative_xk = NULL;
 287 }
 288
 289 static void
 290 rosenbrock_tentative_end (GnmNlsolve *nl, gboolean accept)
 291 {
 292         const int n = nl->n;
 293         GnmIterSolver *isol = nl->isol;
 294
 295         if (!accept && nl->tentative_xk) {
 296                 nl->isol->yk = nl->tentative_yk;
 297                 memcpy (isol->xk, nl->tentative_xk, n * sizeof (gnm_float));
 298         }
 299
 300         nl->tentative = 0;
 301         g_free (nl->tentative_xk);
 302         nl->tentative_xk = NULL;
 303
 304         nl->smallsteps = 0;
 305 }
 306
 307 static gboolean
 308 rosenbrock_iter (GnmNlsolve *nl)
 309 {
 310         GnmSolver *sol = nl->sol;
 311         GnmIterSolver *isol = nl->isol;
 312         const int n = nl->n;
 313         int i, j;
 314         const gnm_float alpha = 3;
 315         const gnm_float beta = 0.5;
 316         gboolean any_at_all = FALSE;
 317         gnm_float *d, **A, *x, *dx, *t;
 318         char *state;
 319         int dones = 0;
 320         gnm_float ykm1 = isol->yk, *xkm1;
 321         gnm_float eps = gnm_pow2 (-16);
 322         int safety = 0;
 323
 324         if (nl->tentative) {
 325                 nl->tentative--;
 326                 if (nl->tentative == 0) {
 327                         if (nl->debug)
 328                                 g_printerr ("Tentative move rejected\n");
 329                         rosenbrock_tentative_end (nl, FALSE);
 330                 }
 331         }
 332
 333         if ((isol->iterations < 20 || isol->iterations % 100 == 0) &&
 334             gnm_solver_has_analytic_hessian (sol)) {
 335                 if (newton_improve (nl, isol->xk))
 336                         return TRUE;
 337         }
 338
 339         if (isol->iterations % 20 == 0) {
 340                 for (i = 0; i < n; i++)
 341                         for (j = 0; j < n; j++)
 342                                 nl->xi[i][j] = (i == j);
 343         }
 344
 345         A = g_new (gnm_float *, n);
 346         for (i = 0; i < n; i++)
 347                 A[i] = g_new (gnm_float, n);
 348
 349         dx = g_new (gnm_float, n);
 350         for (i = 0; i < n; i++)
 351                 dx[i] = 0;
 352
 353         x = g_new (gnm_float, n);
 354         t = g_new (gnm_float, n);
 355
 356         d = g_new (gnm_float, n);
 357         for (i = 0; i < n; i++) {
 358                 d[i] = (isol->xk[i] == 0)
 359                         ? eps
 360                         : gnm_abs (isol->xk[i]) * eps;
 361         }
 362
 363         xkm1 = g_memdup (isol->xk, n * sizeof (gnm_float));
 364
 365         state = g_new0 (char, n);
 366
 367         while (dones < n) {
 368                 /*
 369                  * A safety that shouldn't get hit, but might if the function
 370                  * being optimized is non-deterministic.
 371                  */
 372                 if (safety++ > n * GNM_MANT_DIG)
 373                         break;
 374
 375                 for (i = 0; i < n; i++) {
 376                         gnm_float y;
 377
 378                         if (state[i] == 2)
 379                                 continue;
 380
 381                         /* x = xk + (d[i] * xi[i])  */
 382                         for (j = 0; j < n; j++)
 383                                 x[j] = isol->xk[j] + d[i] * nl->xi[i][j];
 384
 385                         set_vector (nl, x);
 386                         y = get_value (nl);
 387
 388                         if (y <= isol->yk && gnm_solver_check_constraints (sol)) {
 389                                 if (y < isol->yk) {
 390                                         isol->yk = y;
 391                                         memcpy (isol->xk, x, n * sizeof (gnm_float));
 392                                         dx[i] += d[i];
 393                                         any_at_all = TRUE;
 394                                 }
 395                                 switch (state[i]) {
 396                                 case 0:
 397                                         state[i] = 1;
 398                                         /* Fall through */
 399                                 case 1:
 400                                         d[i] *= alpha;
 401                                         break;
 402                                 default:
 403                                 case 2:
 404                                         break;
 405                                 }
 406                         } else {
 407                                 switch (state[i]) {
 408                                 case 1:
 409                                         state[i] = 2;
 410                                         dones++;
 411                                         /* Fall through */
 412                                 case 0:
 413                                         d[i] *= -beta;
 414                                         break;
 415                                 default:
 416                                 case 2:
 417                                         /* No sign change. */
 418                                         d[i] *= 0.5;
 419                                         break;
 420                                 }
 421                         }
 422                 }
 423         }
 424
 425         if (any_at_all) {
 426                 gnm_float div, sum;
 427
 428                 for (j = n - 1; j >= 0; j--)
 429                         for (i = 0; i < n; i++)
 430                                 A[j][i] = (j == n - 1 ? 0 : A[j + 1][i]) + dx[j] * nl->xi[j][i];
 431
 432                 sum = 0;
 433                 for (i = n - 1; i >= 0; i--) {
 434                         sum += dx[i] * dx[i];
 435                         t[i] = sum;
 436                 }
 437
 438                 for (i = n - 1; i > 0; i--) {
 439                         div = gnm_sqrt (t[i - 1] * t[i]);
 440                         if (div != 0)
 441                                 for (j = 0; j < n; j++) {
 442                                         nl->xi[i][j] = (dx[i - 1] * A[i][j] -
 443                                                         nl->xi[i - 1][j] * t[i]) / div;
 444                                         g_assert (gnm_finite (nl->xi[i][j]));
 445                                 }
 446                 }
 447
 448                 gnm_range_hypot (dx, n, &div);
 449                 if (div != 0) {
 450                         for (i = 0; i < n; i++) {
 451                                 nl->xi[0][i] = A[0][i] / div;
 452                                 if (!gnm_finite (nl->xi[0][i])) {
 453                                         g_printerr ("%g %g %g\n",
 454                                                     div, A[0][i], nl->xi[0][i]);
 455                                         g_assert (gnm_finite (nl->xi[0][i]));
 456                                 }
 457                         }
 458                 }
 459
 460                 /* ---------------------------------------- */
 461
 462                 if (!nl->tentative) {
 463                         set_vector (nl, isol->xk);
 464                         set_solution (nl);
 465                 }
 466
 467                 if (nl->tentative) {
 468                         if (isol->yk < nl->tentative_yk) {
 469                                 if (nl->debug)
 470                                         g_printerr ("Tentative move accepted!\n");
 471                                 rosenbrock_tentative_end (nl, TRUE);
 472                         }
 473                 } else if (gnm_abs (isol->yk - ykm1) > gnm_abs (ykm1) * 0.01) {
 474                         /* A big step.  */
 475                         nl->smallsteps = 0;
 476                 } else {
 477                         nl->smallsteps++;
 478                 }
 479         }
 480
 481         g_free (x);
 482         g_free (xkm1);
 483         g_free (dx);
 484         g_free (t);
 485         g_free (d);
 486         free_matrix (A, n);
 487         g_free (state);
 488
 489         return any_at_all;
 490 }
 491
 492 static gboolean
 493 gnm_nlsolve_iterate (GnmSolverIterator *iter, GnmNlsolve *nl)
 494 {
 495         GnmIterSolver *isol = nl->isol;
 496         const int n = nl->n;
 497
 498         if (isol->iterations == 0)
 499                 nlsolve_init (nl);
 500
 501         if (nl->debug) {
 502                 g_printerr ("Iteration %ld at %.15" GNM_FORMAT_g "\n",
 503                             (long)(isol->iterations), isol->yk);
 504                 print_vector ("Current point", isol->xk, n);
 505         }
 506
 507         return rosenbrock_iter (nl);
 508 }
 509
 510 static void
 511 gnm_nlsolve_final (GnmNlsolve *nl)
 512 {
 513         const int n = nl->n;
 514
 515         /* Accept, i.e., don't try to restore.  */
 516         rosenbrock_tentative_end (nl, TRUE);
 517
 518         if (nl->xi) {
 519                 free_matrix (nl->xi, n);
 520                 nl->xi = NULL;
 521         }
 522
 523         g_free (nl);
 524 }
 525
 526 /* ------------------------------------------------------------------------- */
 527 /* Plug-in interface.  */
 528
 529 gboolean nlsolve_solver_factory_functional (GnmSolverFactory *factory);
 530 GnmSolver *nlsolve_solver_factory (GnmSolverFactory *factory, GnmSolverParameters *params);
 531
 532 gboolean
 533 nlsolve_solver_factory_functional (GnmSolverFactory *factory)
 534 {
 535         return TRUE;
 536 }
 537
 538 GnmSolver *
 539 nlsolve_solver_factory (GnmSolverFactory *factory, GnmSolverParameters *params)
 540 {
 541         GnmIterSolver *isol = g_object_new
 542                 (GNM_ITER_SOLVER_TYPE,
 543                  "params", params,
 544                  "flip-sign", (params->problem_type == GNM_SOLVER_MAXIMIZE),
 545                  NULL);
 546         GnmSolver *sol = GNM_SOLVER (isol);
 547         GnmNlsolve *nl = g_new0 (GnmNlsolve, 1);
 548         GnmSolverIteratorCompound *citer;
 549         GnmSolverIterator *iter;
 550
 551         citer = g_object_new (GNM_SOLVER_ITERATOR_COMPOUND_TYPE, NULL);
 552
 553         iter = gnm_solver_iterator_new_func (G_CALLBACK (gnm_nlsolve_iterate), nl);
 554         gnm_solver_iterator_compound_add (citer, iter, 1);
 555
 556         gnm_solver_iterator_compound_add (citer, gnm_solver_iterator_new_polish (isol), 0);
 557
 558         gnm_iter_solver_set_iterator (isol, GNM_SOLVER_ITERATOR (citer));
 559         g_object_unref (citer);
 560
 561         nl->sol = sol;
 562         nl->isol = isol;
 563         nl->debug = gnm_solver_debug ();
 564         nl->min_factor = 1e-10;
 565         nl->n = nl->sol->input_cells->len;
 566
 567         g_signal_connect (isol, "prepare", G_CALLBACK (gnm_nlsolve_prepare), nl);
 568
 569         g_object_set_data_full (G_OBJECT (isol), PRIVATE_KEY, nl,
 570                                 (GDestroyNotify)gnm_nlsolve_final);
 571
 572         return sol;
 573 }
 574
 575 /* ------------------------------------------------------------------------- */