[gnumeric] Solver: add generic line search.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnumeric] Solver: add generic line search.
- Date: Sat, 2 May 2015 21:22:31 +0000 (UTC)
commit 6913d655617a06f8b4b00800fa8336c01f8d1d0a
Author: Morten Welinder <terra gnome org>
Date: Sat May 2 17:22:06 2015 -0400
Solver: add generic line search.
src/tools/ChangeLog | 6 +
src/tools/gnm-solver.c | 256 +++++++++++++++++++++++++++++++++++++-----------
src/tools/gnm-solver.h | 6 +
3 files changed, 211 insertions(+), 57 deletions(-)
---
diff --git a/src/tools/ChangeLog b/src/tools/ChangeLog
index f4a0047..0b40274 100644
--- a/src/tools/ChangeLog
+++ b/src/tools/ChangeLog
@@ -1,3 +1,9 @@
+2015-05-02 Morten Welinder <terra gnome org>
+
+ * gnm-solver.c (cb_polish_iter): Implement in terms of line search.
+ (gnm_solver_line_search): Generic line search with Fibonacci
+ interval reduction.
+
2015-04-28 Morten Welinder <terra gnome org>
* gnm-solver.c (gnm_solver_iterator_new_polish): New function with
diff --git a/src/tools/gnm-solver.c b/src/tools/gnm-solver.c
index 26bd3d0..5ad3cf8 100644
--- a/src/tools/gnm-solver.c
+++ b/src/tools/gnm-solver.c
@@ -1707,8 +1707,8 @@ gnm_solver_set_vars (GnmSolver *sol, gnm_float const *xs)
* @xs: Point to compute gradient at
*
* Returns: (transfer full): A vector containing the gradient. This
- * function takes the flip_sign into account. Note, that this is a
- * numerical approximation.
+ * function takes the flip-sign property into account. Note, that this
+ * is a numerical approximation.
*/
gnm_float *
gnm_solver_compute_gradient (GnmSolver *sol, gnm_float const *xs)
@@ -1743,6 +1743,178 @@ gnm_solver_compute_gradient (GnmSolver *sol, gnm_float const *xs)
return g;
}
+static gnm_float
+try_step (GnmSolver *sol, gnm_float const *x0, gnm_float const *dir, gnm_float step)
+{
+ int const n = sol->input_cells->len;
+ gnm_float *x = g_new (gnm_float, n);
+ int i;
+ gnm_float y;
+
+ for (i = 0; i < n; i++)
+ x[i] = x0[i] + step * dir[i];
+ gnm_solver_set_vars (sol, x);
+ y = gnm_solver_get_target_value (sol);
+ g_free (x);
+ return y;
+}
+
+/**
+ * gnm_solver_line_search:
+ * @sol: Solver
+ * @x0: Starting point
+ * @dir: direction
+ * @try_reverse: whether to try reverse direction at all
+ * @step: initial step size
+ * @max_step: largest allowed step
+ * @eps: tolerance for optimal step
+ * @py: (out): location to store resulting objective function value
+ *
+ * Returns: optimal step size.
+ */
+gnm_float
+gnm_solver_line_search (GnmSolver *sol,
+ gnm_float const *x0, gnm_float const *dir,
+ gboolean try_reverse,
+ gnm_float step, gnm_float max_step, gnm_float eps,
+ gnm_float *py)
+{
+ /*
+ * 0: Initial
+ * 1: Found improvement, but not far side of it
+ * 2: Have (s0,s1,s2) with s1 lowest
+ */
+ int phase = 0;
+ gnm_float s, s0, s1, s2;
+ gnm_float y, y0, y1, y2;
+ gnm_float const phi = (gnm_sqrt (5) + 1) / 2;
+ gboolean rbig;
+ gboolean debug = gnm_solver_debug ();
+
+ g_return_val_if_fail (eps >= 0, gnm_nan);
+ g_return_val_if_fail (step > 0, gnm_nan);
+ g_return_val_if_fail (max_step >= step, gnm_nan);
+
+ if (debug)
+ g_printerr ("LS: step=%" GNM_FORMAT_g ", max=%" GNM_FORMAT_g ", eps=%" GNM_FORMAT_g "\n",
step, max_step, eps);
+
+ gnm_solver_set_vars (sol, x0);
+ y0 = gnm_solver_get_target_value (sol);
+ s0 = 0;
+
+ s = step;
+ while (phase == 0) {
+ gboolean flat = TRUE;
+
+ y = try_step (sol, x0, dir, s);
+ if (debug)
+ g_printerr ("LS0: s:%.6" GNM_FORMAT_g " y:%.6" GNM_FORMAT_g "\n",
+ s, y);
+ if (y < y0 && gnm_solver_check_constraints (sol)) {
+ y1 = y;
+ s1 = s;
+ phase = 1;
+ break;
+ } else if (y != y0)
+ flat = FALSE;
+
+ if (try_reverse) {
+ y = try_step (sol, x0, dir, -s);
+ if (debug)
+ g_printerr ("LS0: s:%.6" GNM_FORMAT_g " y:%.6" GNM_FORMAT_g "\n",
+ -s, y);
+ if (y < y0 && gnm_solver_check_constraints (sol)) {
+ y1 = y;
+ s1 = -s;
+ phase = 1;
+ break;
+ } else if (y != y0)
+ flat = FALSE;
+ }
+
+ s /= 32;
+
+ if (s <= 0 || flat)
+ return gnm_nan;
+ }
+
+ while (phase == 1) {
+ s = s1 * (phi + 1);
+
+ if (gnm_abs (s) >= max_step)
+ return gnm_nan;
+
+ y = try_step (sol, x0, dir, s);
+ if (!gnm_finite (y) || !gnm_solver_check_constraints (sol))
+ return gnm_nan;
+
+ if (y < y1) {
+ y1 = y;
+ s1 = s;
+ continue;
+ }
+
+ y2 = y;
+ s2 = s;
+ phase = 2;
+ }
+
+ /*
+ * Phase 2: we have three steps, s0/s1/s2, in order (descending or ascending) such
+ * that
+ * 1. y1<=y0 (equality unlikely)
+ * 2. y1<=y2 (equality unlikely)
+ * 3a. (s2-s1)=phi*(s1-s0) or
+ * 3b. (s2-s1)*phi=(s1-s0)
+ */
+ rbig = TRUE; /* Initially 3a holds. */
+ while (phase == 2) {
+ if (debug) {
+ gnm_float s01 = s1 - s0;
+ gnm_float s12 = s2 - s1;
+ g_printerr ("LS2: s0:%.6" GNM_FORMAT_g " s01=%.6" GNM_FORMAT_g " s12=%.6"
GNM_FORMAT_g
+ " r=%" GNM_FORMAT_g
+ " y:%.10" GNM_FORMAT_g "/%.10" GNM_FORMAT_g "/%.10" GNM_FORMAT_g "\n",
+ s0, s01, s12, s12 / s01, y0, y1, y2);
+ }
+
+ s = rbig ? s1 + (s1 - s0) * (phi - 1) : s1 - (s2 - s1) * (phi - 1);
+ if (s <= s0 || s >= s2 || gnm_abs (s - s1) <= eps)
+ break;
+
+ y = try_step (sol, x0, dir, s);
+ if (!gnm_finite (y) || !gnm_solver_check_constraints (sol))
+ return gnm_nan;
+
+ if (y < y1) {
+ if (rbig) {
+ y0 = y1;
+ s0 = s1;
+ } else {
+ y2 = y1;
+ s2 = s1;
+ }
+ y1 = y;
+ s1 = s;
+ } else {
+ if (rbig) {
+ y2 = y;
+ s2 = s;
+ } else {
+ y0 = y;
+ s0 = s;
+ }
+ rbig = !rbig;
+
+ if (y0 == y1 && y1 == y2)
+ break;
+ }
+ }
+
+ *py = y1;
+ return s1;
+}
+
static void
gnm_solver_class_init (GObjectClass *object_class)
@@ -2280,67 +2452,32 @@ cb_polish_iter (GnmSolverIterator *iter, GnmIterSolver *isol)
{
GnmSolver *sol = GNM_SOLVER (isol);
const int n = sol->input_cells->len;
- gnm_float *x;
+ gnm_float *dir;
gboolean progress = FALSE;
- gboolean debug = gnm_solver_debug ();
int c;
- x = g_new (gnm_float, n);
-
+ dir = g_new0 (gnm_float, n);
for (c = 0; c < n; c++) {
- gnm_float xc = isol->xk[c], dx;
+ gnm_float s, y, s0, xc = isol->xk[c];
int e;
- gboolean try_reverse = TRUE;
- gboolean try_bigger = TRUE;
(void)gnm_frexp (xc, &e);
- dx = gnm_ldexp (1, e + 10);
- if (dx == 0) dx = GNM_MIN;
-
- memcpy (x, isol->xk, n * sizeof (gnm_float));
- while (1) {
- gnm_float y;
-
- x[c] = xc + dx;
- if (x[c] == xc)
- break;
- gnm_solver_set_vars (sol, x);
- y = gnm_solver_get_target_value (sol);
-
- if (y < isol->yk && gnm_solver_check_constraints (sol)) {
- /* Success! */
-
- isol->yk = y;
- isol->xk[c] = xc = x[c];
- progress = TRUE;
- try_reverse = FALSE;
- if (debug)
- g_printerr ("Polish step %.15" GNM_FORMAT_g
- " in direction %d\n",
- dx, c);
- if (try_bigger) {
- dx *= 2;
- if (gnm_finite (dx))
- continue;
- }
- break;
- } else {
- /* Step didn't work. Restore and find new step to try. */
- x[c] = xc;
-
- if (try_reverse) {
- dx = -dx;
- if (dx < 0)
- continue;
- }
-
- dx /= 2;
- try_bigger = FALSE;
- }
+ s0 = gnm_ldexp (1, e - 10);
+ if (s0 == 0) s0 = GNM_MIN;
+
+ dir[c] = 1;
+ s = gnm_solver_line_search (sol, isol->xk, dir, TRUE,
+ s0, gnm_abs (xc), 0.0,
+ &y);
+ dir[c] = 0;
+
+ if (gnm_finite (s)) {
+ isol->xk[c] += s;
+ isol->yk = y;
+ progress = TRUE;
}
}
-
- g_free (x);
+ g_free (dir);
if (progress)
gnm_iter_solver_set_solution (isol);
@@ -2632,9 +2769,14 @@ gnm_iter_solver_idle (gpointer data)
progress = isol->iterator && gnm_solver_iterator_iterate (isol->iterator);
isol->iterations++;
- if (!gnm_solver_finished (sol) &&
- (!progress || isol->iterations >= params->options.max_iter))
- gnm_solver_set_status (sol, GNM_SOLVER_STATUS_DONE);
+ if (!gnm_solver_finished (sol)) {
+ if (!progress) {
+ gnm_solver_set_status (sol, GNM_SOLVER_STATUS_DONE);
+ } else if (isol->iterations >= params->options.max_iter) {
+ gnm_solver_stop (sol, NULL);
+ gnm_solver_set_reason (sol, _("Iteration limit exceeded"));
+ }
+ }
if (gnm_solver_finished (sol)) {
isol->idle_tag = 0;
diff --git a/src/tools/gnm-solver.h b/src/tools/gnm-solver.h
index 50be250..bd454e4 100644
--- a/src/tools/gnm-solver.h
+++ b/src/tools/gnm-solver.h
@@ -259,6 +259,12 @@ void gnm_solver_set_vars (GnmSolver *sol, gnm_float const *xs);
gnm_float *gnm_solver_compute_gradient (GnmSolver *sol, gnm_float const *xs);
+gnm_float gnm_solver_line_search (GnmSolver *sol,
+ gnm_float const *x0, gnm_float const *dir,
+ gboolean try_reverse,
+ gnm_float step, gnm_float max_step, gnm_float eps,
+ gnm_float *py);
+
/* ------------------------------------------------------------------------- */
/* Solver subclass for subprocesses. */
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
