1 #include <gnumeric-config.h>
3 #include <tools/gnm-solver.h>
7 #include "regression.h"
11 #include <gsf/gsf-impl-utils.h>
12 #include <glib/gi18n-lib.h>
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.
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.
24 * There are a couple of attrictive features of the Rosenbrock method:
25 * 1. It's monotonic. Unlike Newton-style methods it cannot suddenly
27 * 2. We don't need the Hessian.
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.)
34 #define PRIVATE_KEY "::nlsolve::"
37 * Note: the solver code assumes the problem is a minimization problem.
38 * When used for a maximization problem, we flip the objective function
44 /* The solver object in two forms. */
51 /* Rosenbrock state */
55 gnm_float
*tentative_xk
, tentative_yk
;
58 // (nothing right now)
66 check_program (GnmSolver
*sol
, GError
**err
)
69 const GnmSolverParameters
*params
= sol
->params
;
72 for (l
= params
->constraints
; l
; l
= l
->next
) {
73 GnmSolverConstraint
*c
= l
->data
;
77 * This catches also equalities where the sides are not
86 for (ui
= 0; ui
< sol
->input_cells
->len
; ui
++) {
87 if (sol
->discrete
[ui
])
91 * This also catches using two inequality constraints used
92 * to emulate equality.
94 if (sol
->min
[ui
] == sol
->max
[ui
])
104 _("This solver does not handle discrete variables."));
111 _("This solver does not handle equality constraints."));
116 print_vector (const char *name
, const gnm_float
*v
, int n
)
121 g_printerr ("%s:\n", name
);
122 for (i
= 0; i
< n
; i
++)
123 g_printerr ("%15.8" GNM_FORMAT_f
" ", v
[i
]);
129 set_value (GnmNlsolve
*nl
, int i
, gnm_float x
)
131 gnm_solver_set_var (nl
->sol
, i
, x
);
136 set_vector (GnmNlsolve
*nl
, const gnm_float
*xs
)
138 gnm_solver_set_vars (nl
->sol
, xs
);
141 /* Get the target value as-if we were minimizing. */
143 get_value (GnmNlsolve
*nl
)
145 /* nl->sol has been taught to flip sign if needed. */
146 return gnm_solver_get_target_value (nl
->sol
);
150 free_matrix (gnm_float
**m
, int n
)
153 for (i
= 0; i
< n
; i
++)
159 set_solution (GnmNlsolve
*nl
)
161 /* nl->isol has been taught to flip sign if needed. */
162 gnm_iter_solver_set_solution (nl
->isol
);
166 gnm_nlsolve_prepare (GnmSolver
*sol
, WorkbookControl
*wbc
, GError
**err
,
171 g_return_val_if_fail (sol
->status
== GNM_SOLVER_STATUS_READY
, FALSE
);
173 gnm_solver_set_status (sol
, GNM_SOLVER_STATUS_PREPARING
);
175 ok
= check_program (sol
, err
);
177 ok
= gnm_iter_solver_get_initial_solution (nl
->isol
, err
);
180 gnm_solver_set_status (sol
, GNM_SOLVER_STATUS_PREPARED
);
182 gnm_solver_set_status (sol
, GNM_SOLVER_STATUS_ERROR
);
189 compute_gradient (GnmNlsolve
*nl
, const gnm_float
*xs
)
191 return gnm_solver_compute_gradient (nl
->sol
, xs
);
195 newton_improve (GnmNlsolve
*nl
, gnm_float
*xs
)
197 GnmSolver
*sol
= nl
->sol
;
198 GnmIterSolver
*isol
= nl
->isol
;
201 gnm_float
*g
, *d
, *xs2
;
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
);
211 for (i
= 0; i
< n
; 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
);
222 print_vector ("d", d
, n
);
224 g_printerr ("Failed to solve Newton step.\n");
230 for (i
= 0; i
< n
; i
++)
231 xs2
[i
] = xs
[i
] + d
[i
];
232 set_vector (nl
, xs2
);
236 if (y2
< isol
->yk
&& gnm_solver_check_constraints (sol
)) {
238 g_printerr ("Accepting newton step\n");
239 memcpy (isol
->xk
, xs2
, n
* sizeof (gnm_float
));
245 g_printerr ("Full newton step would go to %g\n", y2
);
247 f
= gnm_solver_line_search
248 (sol
, xs
, d
, FALSE
, 0.75, 1, 0.01, &y2
);
250 if (f
> 0 && f
< 1 && y2
< isol
->yk
) {
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
];
271 nlsolve_init (GnmNlsolve
*nl
)
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
);
286 nl
->tentative_xk
= NULL
;
290 rosenbrock_tentative_end (GnmNlsolve
*nl
, gboolean accept
)
293 GnmIterSolver
*isol
= nl
->isol
;
295 if (!accept
&& nl
->tentative_xk
) {
296 nl
->isol
->yk
= nl
->tentative_yk
;
297 memcpy (isol
->xk
, nl
->tentative_xk
, n
* sizeof (gnm_float
));
301 g_free (nl
->tentative_xk
);
302 nl
->tentative_xk
= NULL
;
308 rosenbrock_iter (GnmNlsolve
*nl
)
310 GnmSolver
*sol
= nl
->sol
;
311 GnmIterSolver
*isol
= nl
->isol
;
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
;
320 gnm_float ykm1
= isol
->yk
, *xkm1
;
321 gnm_float eps
= gnm_pow2 (-16);
326 if (nl
->tentative
== 0) {
328 g_printerr ("Tentative move rejected\n");
329 rosenbrock_tentative_end (nl
, FALSE
);
333 if ((isol
->iterations
< 20 || isol
->iterations
% 100 == 0) &&
334 gnm_solver_has_analytic_hessian (sol
)) {
335 if (newton_improve (nl
, isol
->xk
))
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
);
345 A
= g_new (gnm_float
*, n
);
346 for (i
= 0; i
< n
; i
++)
347 A
[i
] = g_new (gnm_float
, n
);
349 dx
= g_new (gnm_float
, n
);
350 for (i
= 0; i
< n
; i
++)
353 x
= g_new (gnm_float
, n
);
354 t
= g_new (gnm_float
, n
);
356 d
= g_new (gnm_float
, n
);
357 for (i
= 0; i
< n
; i
++) {
358 d
[i
] = (isol
->xk
[i
] == 0)
360 : gnm_abs (isol
->xk
[i
]) * eps
;
363 xkm1
= g_memdup (isol
->xk
, n
* sizeof (gnm_float
));
365 state
= g_new0 (char, n
);
369 * A safety that shouldn't get hit, but might if the function
370 * being optimized is non-deterministic.
372 if (safety
++ > n
* GNM_MANT_DIG
)
375 for (i
= 0; i
< n
; i
++) {
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
];
388 if (y
<= isol
->yk
&& gnm_solver_check_constraints (sol
)) {
391 memcpy (isol
->xk
, x
, n
* sizeof (gnm_float
));
417 /* No sign change. */
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
];
433 for (i
= n
- 1; i
>= 0; i
--) {
434 sum
+= dx
[i
] * dx
[i
];
438 for (i
= n
- 1; i
> 0; i
--) {
439 div
= gnm_sqrt (t
[i
- 1] * t
[i
]);
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
]));
448 gnm_range_hypot (dx
, n
, &div
);
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
]));
460 /* ---------------------------------------- */
462 if (!nl
->tentative
) {
463 set_vector (nl
, isol
->xk
);
468 if (isol
->yk
< nl
->tentative_yk
) {
470 g_printerr ("Tentative move accepted!\n");
471 rosenbrock_tentative_end (nl
, TRUE
);
473 } else if (gnm_abs (isol
->yk
- ykm1
) > gnm_abs (ykm1
) * 0.01) {
493 gnm_nlsolve_iterate (GnmSolverIterator
*iter
, GnmNlsolve
*nl
)
495 GnmIterSolver
*isol
= nl
->isol
;
498 if (isol
->iterations
== 0)
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
);
507 return rosenbrock_iter (nl
);
511 gnm_nlsolve_final (GnmNlsolve
*nl
)
515 /* Accept, i.e., don't try to restore. */
516 rosenbrock_tentative_end (nl
, TRUE
);
519 free_matrix (nl
->xi
, n
);
526 /* ------------------------------------------------------------------------- */
527 /* Plug-in interface. */
529 gboolean
nlsolve_solver_factory_functional (GnmSolverFactory
*factory
);
530 GnmSolver
*nlsolve_solver_factory (GnmSolverFactory
*factory
, GnmSolverParameters
*params
);
533 nlsolve_solver_factory_functional (GnmSolverFactory
*factory
)
539 nlsolve_solver_factory (GnmSolverFactory
*factory
, GnmSolverParameters
*params
)
541 GnmIterSolver
*isol
= g_object_new
542 (GNM_ITER_SOLVER_TYPE
,
544 "flip-sign", (params
->problem_type
== GNM_SOLVER_MAXIMIZE
),
546 GnmSolver
*sol
= GNM_SOLVER (isol
);
547 GnmNlsolve
*nl
= g_new0 (GnmNlsolve
, 1);
548 GnmSolverIteratorCompound
*citer
;
549 GnmSolverIterator
*iter
;
551 citer
= g_object_new (GNM_SOLVER_ITERATOR_COMPOUND_TYPE
, NULL
);
553 iter
= gnm_solver_iterator_new_func (G_CALLBACK (gnm_nlsolve_iterate
), nl
);
554 gnm_solver_iterator_compound_add (citer
, iter
, 1);
556 gnm_solver_iterator_compound_add (citer
, gnm_solver_iterator_new_polish (isol
), 0);
558 gnm_iter_solver_set_iterator (isol
, GNM_SOLVER_ITERATOR (citer
));
559 g_object_unref (citer
);
563 nl
->debug
= gnm_solver_debug ();
564 nl
->min_factor
= 1e-10;
565 nl
->n
= nl
->sol
->input_cells
->len
;
567 g_signal_connect (isol
, "prepare", G_CALLBACK (gnm_nlsolve_prepare
), nl
);
569 g_object_set_data_full (G_OBJECT (isol
), PRIVATE_KEY
, nl
,
570 (GDestroyNotify
)gnm_nlsolve_final
);
575 /* ------------------------------------------------------------------------- */