diff --git a/goffice/math/go-regression.c b/goffice/math/go-regression.c
index f207fe7..c70ef1e 100644
--- a/goffice/math/go-regression.c
+++ b/goffice/math/go-regression.c
@@ -577,12 +577,14 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xssT, int n,
 	QUAD *qresult;
 	QUAD *QTy;
 	QUAD *inv;
+	QUAD *qys;
 	int i, j, k;
 	gboolean has_result;
 	void *state;
 	gboolean has_stat;
 	int err;
 	QUAD N2;
+	DOUBLE threshold = DEFAULT_THRESHOLD;
 
 	ZERO_VECTOR (result, n);
 
@@ -597,12 +599,15 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xssT, int n,
 
 	state = SUFFIX(go_quad_start) ();
 
+	qys = g_new (QUAD, m);
+	for (i = 0; i < m; i++)
+		SUFFIX(go_quad_init)(&qys[i], ys[i]);
+
 	xss = SUFFIX(go_quad_matrix_new) (m, n);
 	for (i = 0; i < m; i++)
 		for (j = 0; j < n; j++)
 			SUFFIX(go_quad_init) (&xss->data[i][j], xssT[j][i]);
 	qr = SUFFIX(go_quad_qr_new) (xss);
-	SUFFIX(go_quad_matrix_free) (xss);
 
 	qresult = g_new0 (QUAD, n);
 	QTy = g_new (QUAD, m);
@@ -611,21 +616,42 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xssT, int n,
 	has_result = (qr != NULL);
 
 	if (has_result) {
+		DOUBLE emin, emax;
 		const SUFFIX(GOQuadMatrix) *R = SUFFIX(go_quad_qr_r) (qr);
 		regerr = GO_REG_ok;
 
+		SUFFIX(go_quad_matrix_eigen_range) (R, &emin, &emax);
+
 		/* Compute Q^T ys.  */
-		for (i = 0; i < m; i++)
-			SUFFIX(go_quad_init)(&QTy[i], ys[i]);
+		COPY_VECTOR (QTy, qys, m);
 		SUFFIX(go_quad_qr_multiply_qt)(qr, QTy);
 
 		/* Solve R res = Q^T ys */
-		if (SUFFIX(go_quad_matrix_back_solve) (R, qresult, QTy))
-			has_result = FALSE;
+		if (emin <= emax * threshold ||
+		    SUFFIX(go_quad_matrix_back_solve) (R, qresult, QTy)) {
+			/* That failed.  Retry via pseudo-inverse.  */
+			SUFFIX(GOQuadMatrix) *Z =
+				SUFFIX(go_quad_matrix_pseudo_inverse) (xss, threshold);
+			for (i = 0; i < n; i++) {
+				for (j = 0; j < m; j++)
+					g_printerr ("%2d %2d  %" FORMAT_g "\n",
+						    i, j, SUFFIX(go_quad_value)(&Z->data[i][j]));
+			}
+
+
+			for (i = 0; i < n; i++)
+				SUFFIX(go_quad_dot_product) (&qresult[i],
+							     Z->data[i],
+							     qys,
+							     m);
+
+			SUFFIX(go_quad_matrix_free) (Z);
+		}
 
 		/* Round to plain precision.  */
 		for (i = 0; i < n; i++) {
 			result[i] = SUFFIX(go_quad_value) (&qresult[i]);
+			g_printerr ("%2d: %" FORMAT_g "\n", i, result[i]);
 			SUFFIX(go_quad_init) (&qresult[i], result[i]);
 		}
 
@@ -658,6 +684,7 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xssT, int n,
 			 * This is an indication that something has gone wrong
 			 * numerically.
 			 */
+			g_printerr ("Hi: %.20" FORMAT_g "\n", stat_->ss_total);
 			regerr = GO_REG_near_singular_bad;
 		}
 
@@ -683,12 +710,14 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xssT, int n,
 
 			/* Solve R^T inv = e_k */
 			if (SUFFIX(go_quad_matrix_fwd_solve) (R, inv, inv)) {
+				g_printerr ("Hi1\n");
 				regerr = GO_REG_singular;
 				break;
 			}
 
 			/* Solve R newinv = inv */
 			if (SUFFIX(go_quad_matrix_back_solve) (R, inv, inv)) {
+				g_printerr ("Hi2\n");
 				regerr = GO_REG_singular;
 				break;
 			}
@@ -742,8 +771,11 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xssT, int n,
 	g_free (QTy);
 	g_free (qresult);
 	g_free (inv);
+	g_free (qys);
 	if (qr) SUFFIX(go_quad_qr_free) (qr);
+	SUFFIX(go_quad_matrix_free) (xss);
 out:
+
 	if (!has_stat)
 		SUFFIX(go_regression_stat_destroy) (stat_);
 #endif
