GNOME.org
 

[gnumeric] nlsolve: use modified Rosenbrock descent.

commit 26c363abb5d2c6eebecafac5ce6ed60573018e57
Author: Morten Welinder <terra gnome org>
Date:   Sat May 29 17:30:01 2010 -0400

    nlsolve: use modified Rosenbrock descent.

 NEWS                          |    2 +-
 plugins/nlsolve/gnm-nlsolve.c |  472 +++++++++++++++++++++++++++++++++--------
 2 files changed, 380 insertions(+), 94 deletions(-)
---
diff --git a/NEWS b/NEWS
index 8c1c1f7..207151d 100644
--- a/NEWS
+++ b/NEWS
@@ -6,7 +6,7 @@ Andreas:
 Morten:
 	* Fix stf crash.  [#619283]
 	* Persist solver model type.
-	* Add non-linear solver.
+	* Add non-linear solver.  [#614865]
 	* Feed the gtk+ deprecation monster.  [#572240]
 	* Fix win32/paradox problem.  [#619942]
 
diff --git a/plugins/nlsolve/gnm-nlsolve.c b/plugins/nlsolve/gnm-nlsolve.c
index 66db640..eaecef0 100644
--- a/plugins/nlsolve/gnm-nlsolve.c
+++ b/plugins/nlsolve/gnm-nlsolve.c
@@ -4,8 +4,8 @@
 #include "cell.h"
 #include "sheet.h"
 #include "value.h"
-#include "ranges.h"
 #include "regression.h"
+#include "rangefunc.h"
 #include "workbook.h"
 #include <gsf/gsf-impl-utils.h>
 #include <glib/gi18n-lib.h>
@@ -19,11 +19,21 @@ typedef struct {
 	/* Input/output cells.  */
 	GPtrArray *vars;
 	GnmCell *target;
+	GnmCellPos origin;
+	int input_width, input_height;
 
-	int iterations;
+	/* Initial point.  */
+	gnm_float *x0;
 
-	/* Next axis direction to try.  */
-	int dim;
+	/* Current point.  */
+	gnm_float *xk, yk;
+	int k;
+
+	/* Rosenbrock state */
+	gnm_float **xi;
+	int smallsteps;
+	int tentative;
+	gnm_float *tentative_xk, tentative_yk;
 
 	/* Parameters: */
 	gboolean debug;
@@ -40,17 +50,16 @@ gnm_nlsolve_cleanup (GnmNlsolve *nl)
 		g_source_remove (nl->idle_tag);
 		nl->idle_tag = 0;
 	}
-
-	if (nl->vars) {
-		g_ptr_array_free (nl->vars, TRUE);
-		nl->vars = NULL;
-	}
 }
 
 static void
 gnm_nlsolve_final (GnmNlsolve *nl)
 {
 	gnm_nlsolve_cleanup (nl);
+	if (nl->vars)
+		g_ptr_array_free (nl->vars, TRUE);
+	g_free (nl->xk);
+	g_free (nl->x0);
 	g_free (nl);
 }
 
@@ -139,10 +148,20 @@ gnm_nlsolve_set_solution (GnmNlsolve *nl)
 {
 	GnmSolver *sol = nl->parent;
 	GnmSolverResult *result = g_object_new (GNM_SOLVER_RESULT_TYPE, NULL);
+	const int n = nl->vars->len;
+	int i;
 
-	result->solution = gnm_solver_get_current_values (sol);
 	result->quality = GNM_SOLVER_RESULT_FEASIBLE;
-	result->value = get_value (nl);
+	result->value = nl->yk;
+	result->solution = value_new_array_empty (nl->input_width,
+						  nl->input_height);
+	for (i = 0; i < n; i++) {
+		GnmCell *cell = g_ptr_array_index (nl->vars, i);
+		value_array_set (result->solution,
+				 cell->pos.col - nl->origin.col,
+				 cell->pos.row - nl->origin.row,
+				 value_new_float (nl->xk[i]));
+	}
 
 	g_object_set (sol, "result", result, NULL);
 	g_object_unref (result);
@@ -152,6 +171,8 @@ static gboolean
 gnm_nlsolve_get_initial_solution (GnmNlsolve *nl, GError **err)
 {
 	GnmSolver *sol = nl->parent;
+	const int n = nl->vars->len;
+	int i;
 
 	if (gnm_solver_check_constraints (sol))
 		goto got_it;
@@ -165,7 +186,13 @@ gnm_nlsolve_get_initial_solution (GnmNlsolve *nl, GError **err)
 	return FALSE;
 
 got_it:
+	for (i = 0; i < n; i++) {
+		GnmCell *cell = g_ptr_array_index (nl->vars, i);
+		nl->xk[i] = nl->x0[i] = value_get_as_float (cell->value);
+	}
+	nl->yk = get_value (nl);
 	gnm_nlsolve_set_solution (nl);
+
 	return TRUE;
 }
 
@@ -184,9 +211,6 @@ gnm_nlsolve_prepare (GnmSolver *sol, WorkbookControl *wbc, GError **err,
 		ok = gnm_nlsolve_get_initial_solution (nl, err);
 
 	if (ok) {
-		if (nl->debug)
-			g_printerr ("Initial value:\n%15.8" GNM_FORMAT_f "\n",
-				    sol->result->value);
 		gnm_solver_set_status (sol, GNM_SOLVER_STATUS_PREPARED);
 	} else {
 		gnm_nlsolve_cleanup (nl);
@@ -233,9 +257,6 @@ compute_gradient (GnmNlsolve *nl, const gnm_float *xs)
 
 		set_value (nl, i, x0 + dx);
 		y1 = get_value (nl);
-#if 0
-		g_printerr ("  y0=%g   y1=%g\n", y0, y1);
-#endif
 		g[i] = (y1 - y0) / dx;
 
 		set_value (nl, i, x0);
@@ -288,129 +309,377 @@ compute_hessian (GnmNlsolve *nl, const gnm_float *xs, const gnm_float *g0)
 }
 
 static gboolean
-try_direction (GnmNlsolve *nl, const gnm_float *x0, const gnm_float *d)
+newton_improve (GnmNlsolve *nl, gnm_float *xs, gnm_float *y, gnm_float ymax)
 {
-	GnmSolver *sol = nl->parent;
 	const int n = nl->vars->len;
-	gnm_float factor;
+	gnm_float *g, **H, *d;
+	gboolean ok;
 
-	for (factor = 1; factor >= nl->min_factor; factor /= 2) {
-		gnm_float y;
+	g = compute_gradient (nl, xs);
+	H = compute_hessian (nl, xs, g);
+	d = g_new (gnm_float, n);
+	ok = (gnm_linear_solve (H, g, n, d) == 0);
+
+	if (ok) {
 		int i;
+		gnm_float y2, *xs2 = g_new (gnm_float, n);
+		gnm_float f, best_f = -1;
 
-		for (i = 0; i < n; i++) {
-			gnm_float x = x0[i] - factor * d[i];
-			set_value (nl, i, x);
-		}
+		ok = FALSE;
+		for (f = 1; f > 1e-4; f /= 2) {
+			int i;
+			for (i = 0; i < n; i++)
+				xs2[i] = xs[i] - f * d[i];
+			set_vector (nl, xs2);
+			y2 = get_value (nl);
+			if (nl->debug) {
+				print_vector ("xs2", xs2, n);
+				g_printerr ("YYY %.15g\n", y2);
+			}
 
-		y = get_value (nl);
-		if (nl->debug)
-			g_printerr ("Factor=%-10" GNM_FORMAT_g
-				    "  value=%.8" GNM_FORMAT_f "\n",
-				    factor, y);
-
-		if (y < sol->result->value) {
-			if (gnm_solver_check_constraints (sol)) {
-				gnm_nlsolve_set_solution (nl);
-				return TRUE;
+			if (y2 < ymax) {
+				best_f = f;
+				ymax = y2;
+				break;
 			}
-			if (nl->debug)
-				g_printerr ("Would-be solution does not satisfy constraints.\n");
 		}
+
+		if (best_f > 0) {
+			for (i = 0; i < n; i++)
+				xs[i] = xs[i] - best_f * d[i];
+			*y = ymax;
+			ok = TRUE;
+		}
+
+		g_free (xs2);
+	} else {
+		if (nl->debug)
+			g_printerr ("Failed to solve Newton step.\n");
 	}
 
-	return FALSE;
+	g_free (d);
+	g_free (g);
+	free_matrix (H, n);
+
+	return ok;
 }
 
-static gint
-gnm_nlsolve_idle (gpointer data)
+static void
+rosenbrock_init (GnmNlsolve *nl)
 {
-	GnmNlsolve *nl = data;
-	GnmSolver *sol = nl->parent;
-	gnm_float **H, *g, *d, *x0;
-	int i;
 	const int n = nl->vars->len;
-	gboolean ok;
-	gnm_float y0;
-	gboolean call_again = TRUE;
+	int i, j;
 
-	nl->iterations++;
+	nl->xi = g_new (gnm_float *, n);
+	for (i = 0; i < n; i++) {
+		nl->xi[i] = g_new (gnm_float, n);
+		for (j = 0; j < n; j++)
+			nl->xi[i][j] = (i == j);
+	}
 
-	y0 = get_value (nl);
+	nl->smallsteps = 0;
 
-	x0 = g_new (gnm_float, n);
-	for (i = 0; i < n; i++) {
-		GnmCell *cell = g_ptr_array_index (nl->vars, i);
-		x0[i] = value_get_as_float (cell->value);
+	nl->tentative = 0;
+	nl->tentative_xk = NULL;
+}
+
+static void
+rosenbrock_tentative_end (GnmNlsolve *nl, gboolean accept)
+{
+	const int n = nl->vars->len;
+
+	if (!accept && nl->tentative_xk) {
+		nl->yk = nl->tentative_yk;
+		memcpy (nl->xk, nl->tentative_xk, n * sizeof (gnm_float));
 	}
-	g = compute_gradient (nl, x0);
-	H = compute_hessian (nl, x0, g);
 
-	if (nl->debug) {
-		int i;
+	nl->tentative = 0;
+	g_free (nl->tentative_xk);
+	nl->tentative_xk = NULL;
 
-		print_vector ("x0", x0, n);
-		print_vector ("Gradient", g, n);
+	nl->smallsteps = 0;
+}
 
-		g_printerr ("Hessian:\n");
-		for (i = 0; i < n; i++) {
-			print_vector (NULL, H[i], n);
+static gboolean
+rosenbrock_iter (GnmNlsolve *nl)
+{
+	const int n = nl->vars->len;
+	int i, j;
+	const gnm_float alpha = 3;
+	const gnm_float beta = 0.5;
+	gboolean any_at_all = FALSE;
+	gnm_float *d, **A, *x, *dx, *t;
+	char *state;
+	int dones = 0;
+	gnm_float ykm1 = nl->yk, *xkm1;
+	gnm_float eps = gnm_pow2 (-16);
+
+	if (nl->tentative) {
+		nl->tentative--;
+		if (nl->tentative == 0) {
+			if (nl->debug)
+				g_printerr ("Tentative move rejected\n");
+			rosenbrock_tentative_end (nl, FALSE);
 		}
 	}
 
+	if (nl->k % 20 == 0) {
+		for (i = 0; i < n; i++)
+			for (j = 0; j < n; j++)
+				nl->xi[i][j] = (i == j);
+	}
+
+	A = g_new (gnm_float *, n);
+	for (i = 0; i < n; i++)
+		A[i] = g_new (gnm_float, n);
+
+	dx = g_new (gnm_float, n);
+	for (i = 0; i < n; i++)
+		dx[i] = 0;
+
+	x = g_new (gnm_float, n);
+	t = g_new (gnm_float, n);
+
 	d = g_new (gnm_float, n);
-	ok = (gnm_linear_solve (H, g, n, d) == 0);
-	if (ok) {
-		if (nl->debug)
-			print_vector ("Delta", d, n);
-		ok = try_direction (nl, x0, d);
+	for (i = 0; i < n; i++) {
+		d[i] = (nl->xk[i] == 0)
+			? eps
+			: gnm_abs (nl->xk[i]) * eps;
 	}
 
-	if (!ok) {
-		int i, j;
-		gnm_float *d1 = g_new (gnm_float, n);
-
-		/*
-		 * The Newton step failed.  This is likely because the
-		 * surface is seriously non-linear.  Try one axis
-		 * direction at a time.
-		 */
-
-		for (j = 0; !ok && j < n; j++) {
-			int c = nl->dim++ % n;
-			if (g[c] == 0)
+	xkm1 = g_memdup (nl->xk, n * sizeof (gnm_float));
+
+	state = g_new0 (char, n);
+
+	while (dones < n) {
+		for (i = 0; i < n; i++) {
+			gnm_float y;
+
+			if (state[i] == 2)
 				continue;
 
+			/* x = xk + (d[i] * xi[i])  */
+			for (j = 0; j < n; j++)
+				x[j] = nl->xk[j] + d[i] * nl->xi[i][j];
+
+			set_vector (nl, x);
+			y = get_value (nl);
+
+			if (y <= nl->yk) {
+				if (y < nl->yk) {
+					nl->yk = y;
+					memcpy (nl->xk, x, n * sizeof (gnm_float));
+					dx[i] += d[i];
+					any_at_all = TRUE;
+				}
+				switch (state[i]) {
+				case 0:
+					state[i] = 1;
+					/* Fall through */
+				case 1:
+					d[i] *= alpha;
+					break;
+				default:
+				case 2:
+					break;
+				}
+			} else {
+				switch (state[i]) {
+				case 1:
+					state[i] = 2;
+					dones++;
+					/* Fall through */
+				case 0:
+					d[i] *= -beta;
+					break;
+				default:
+				case 2:
+					/* No sign change. */
+					d[i] *= 0.5;
+					break;
+				}
+			}
+		}
+	}
+
+	if (any_at_all) {
+		gnm_float div, sum;
+
+                for (j = n - 1; j >= 0; j--)
 			for (i = 0; i < n; i++)
-				d1[i] = (i == c)
-					? y0 / g[c]
-					: 0;
+				A[j][i] = (j == n - 1 ? 0 : A[j + 1][i]) + dx[j] * nl->xi[j][i];
+
+		sum = 0;
+                for (i = n - 1; i >= 0; i--) {
+			sum += dx[i] * dx[i];
+			t[i] = sum;
+		}
 
-			ok = try_direction (nl, x0, d1);
+                for (i = n - 1; i > 0; i--) {
+			div = gnm_sqrt (t[i - 1] * t[i]);
+			if (div != 0)
+				for (j = 0; j < n; j++) {
+					nl->xi[i][j] = (dx[i - 1] * A[i][j] -
+							nl->xi[i - 1][j] * t[i]) / div;
+					g_assert (gnm_finite (nl->xi[i][j]));
+				}
+                }
+
+                div = sqrt (t[0]);
+                for (i = 0; i < n; i++) {
+			nl->xi[0][i] = A[0][i] / div;
+			g_assert (gnm_finite (nl->xi[0][i]));
 		}
 
-		g_free (d1);
+		/* ---------------------------------------- */
+
+		if (!nl->tentative)
+			gnm_nlsolve_set_solution (nl);
+
+		if (nl->tentative) {
+			if (nl->yk < nl->tentative_yk) {
+				if (nl->debug)
+					g_printerr ("Tentative move accepted!\n");
+				rosenbrock_tentative_end (nl, TRUE);
+			}
+		} else if (gnm_abs (nl->yk - ykm1) > gnm_abs (ykm1) * 0.01) {
+			/* A big step.  */
+			nl->smallsteps = 0;
+		} else {
+			nl->smallsteps++;
+		}
+
+		if (!nl->tentative && nl->smallsteps > 50) {
+			gnm_float yk = nl->yk;
+
+			nl->tentative = 10;
+			nl->tentative_xk = g_memdup (nl->xk, n * sizeof (gnm_float));
+			nl->tentative_yk = yk;
+
+			for (i = 0; i < 4; i++) {
+				gnm_float ymax = yk +
+					gnm_abs (yk) * (0.10 / (i + 1));
+				if (!newton_improve (nl, nl->xk, &nl->yk, ymax))
+					break;
+			}
+
+			if (nl->debug)
+				print_vector ("Tentative move to", nl->xk, n);
+		}
+	}
+
+	g_free (x);
+	g_free (xkm1);
+	g_free (dx);
+	g_free (t);
+	g_free (d);
+	free_matrix (A, n);
+	g_free (state);
+
+	return any_at_all;
+}
+
+static void
+rosenbrock_shutdown (GnmNlsolve *nl)
+{
+	const int n = nl->vars->len;
+
+	rosenbrock_tentative_end (nl, FALSE);
+
+	free_matrix (nl->xi, n);
+	nl->xi = NULL;
+}
+
+static gboolean
+polish_iter (GnmNlsolve *nl)
+{
+	const int n = nl->vars->len;
+	gnm_float *x;
+	gnm_float step;
+	gboolean any_at_all = FALSE;
+
+	x = g_new (gnm_float, n);
+	for (step = gnm_pow2 (-10); step > GNM_EPSILON; step *= 0.75) {
+		gboolean any = FALSE;
+		int c, s;
+
+		for (c = 0; c < n; c++) {
+			for (s = 0; s <= 1; s++) {
+				gnm_float y;
+				gnm_float dx = step * gnm_abs (nl->xk[c]);
+
+				if (dx == 0) dx = step;
+				if (s) dx = -dx;
+
+				memcpy (x, nl->xk, n * sizeof (gnm_float));
+				x[c] += dx;
+				set_vector (nl, x);
+				y = get_value (nl);
+
+				if (y < nl->yk)  {
+					nl->yk = y;
+					memcpy (nl->xk, x, n * sizeof (gnm_float));
+					any = TRUE;
+					any_at_all = TRUE;
+					if (nl->debug)
+						g_printerr ("Polish step %.15" GNM_FORMAT_g
+							    " in direction %d\n",
+							    dx, c);
+					break;
+				}
+			}
+		}
+	}
+
+	g_free (x);
+
+	if (any_at_all)
+		gnm_nlsolve_set_solution (nl);
+
+	return any_at_all;
+}
+
+static gint
+gnm_nlsolve_idle (gpointer data)
+{
+	GnmNlsolve *nl = data;
+	GnmSolver *sol = nl->parent;
+	const int n = nl->vars->len;
+	gboolean ok;
+	gboolean call_again = TRUE;
+
+	if (nl->k == 0)
+		rosenbrock_init (nl);
+
+	if (nl->debug) {
+		g_printerr ("Iteration %d at %.15" GNM_FORMAT_g "\n",
+			    nl->k, nl->yk);
+		print_vector ("Current point", nl->xk, n);
+	}
+
+	nl->k++;
+	ok = rosenbrock_iter (nl);
+
+	if (!ok && !nl->tentative) {
+		ok = polish_iter (nl);
 	}
 
 	if (!ok) {
-		set_vector (nl, x0);
 		gnm_solver_set_status (sol, GNM_SOLVER_STATUS_DONE);
 		call_again = FALSE;
 	}
 
-	if (call_again && nl->iterations >= nl->max_iter) {
+	if (call_again && nl->k >= nl->max_iter) {
 		gnm_solver_set_status (sol, GNM_SOLVER_STATUS_DONE);
 		call_again = FALSE;
 	}
 
 	if (!call_again) {
+		set_vector (nl, nl->x0);
 		workbook_recalc (sol->params->sheet->workbook);
-	}
 
-	g_free (d);
-	g_free (x0);
-	g_free (g);
-	free_matrix (H, n);
+		rosenbrock_shutdown (nl);
+	}
 
 	return call_again;
 }
@@ -462,9 +731,22 @@ nlsolve_solver_factory (GnmSolverFactory *factory, GnmSolverParameters *params)
 				       NULL);
 	GnmNlsolve *nl = g_new0 (GnmNlsolve, 1);
 	GSList *input_cells, *l;
+	int n;
+	GnmValue const *vinput = gnm_solver_param_get_input (params);
+	GnmEvalPos ep;
+	GnmCellRef origin;
 
 	nl->parent = GNM_SOLVER (res);
 
+	eval_pos_init_sheet (&ep, params->sheet);
+	if (vinput) {
+		gnm_cellref_make_abs (&origin, &vinput->v_range.cell.a, &ep);
+		nl->origin.col = origin.col;
+		nl->origin.row = origin.row;
+		nl->input_width = value_area_get_width (vinput, &ep);
+		nl->input_height = value_area_get_height (vinput, &ep);
+	}
+
 	nl->debug = gnm_solver_debug ();
 	nl->max_iter = params->options.max_iter;
 	nl->min_factor = 1e-10;
@@ -476,6 +758,10 @@ nlsolve_solver_factory (GnmSolverFactory *factory, GnmSolverParameters *params)
 	for (l = input_cells; l; l = l->next)
 		g_ptr_array_add (nl->vars, l->data);
 	g_slist_free (input_cells);
+	n = nl->vars->len;
+
+	nl->x0 = g_new (gnm_float, n);
+	nl->xk = g_new (gnm_float, n);
 
 	g_signal_connect (res, "prepare", G_CALLBACK (gnm_nlsolve_prepare), nl);
 	g_signal_connect (res, "start", G_CALLBACK (gnm_nlsolve_start), nl);
[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]