[goffice] Regression: first cut at pseudo-inverse calculation.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [goffice] Regression: first cut at pseudo-inverse calculation.
- Date: Tue, 30 Apr 2013 20:26:48 +0000 (UTC)
commit a7d66ce39d1722f975306409bcd1bd8474cbd8d5
Author: Morten Welinder <terra gnome org>
Date: Tue Apr 30 16:26:16 2013 -0400
Regression: first cut at pseudo-inverse calculation.
ChangeLog | 5 +
goffice/math/go-regression.c | 179 +++++++++++++++++++++++++++++++++++++++++-
goffice/math/go-regression.h | 8 ++
3 files changed, 190 insertions(+), 2 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 6b1e47e..c7077fe 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2013-04-30 Morten Welinder <terra gnome org>
+
+ * goffice/math/go-regression.c (go_matrix_pseudo_inverse): New
+ function.
+
2013-04-26 Morten Welinder <terra gnome org>
* configure.ac: Post-release bump.
diff --git a/goffice/math/go-regression.c b/goffice/math/go-regression.c
index 7a72429..ed681d7 100644
--- a/goffice/math/go-regression.c
+++ b/goffice/math/go-regression.c
@@ -307,9 +307,11 @@ copy_stats (go_regression_stat_t *s1,
*
* V is a pre-allocated output m-times-n matrix. V will contrain
* n vectors of different lengths: n, n-1, ..., 1. These are the
- * Householder vectors (or null for the degenerate case).
+ * Householder vectors (or null for the degenerate case). The
+ * matrix Q of size m-times-n is implied from V.
*
- * Rfinal is a pre-allocated output matrix of size n-times-n.
+ * Rfinal is a pre-allocated output matrix of size n-times-n. (To
+ * get the m-times-n version of R, simply add m-n null rows.)
*
* qdet is an optional output location for det(Q).
*
@@ -2134,3 +2136,176 @@ SUFFIX(go_linear_regression_leverage) (MATRIX A, DOUBLE *d, int m, int n)
return regres;
}
+
+/*
+ * Compute the pseudo-inverse of a matrix, see
+ * http://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse
+ *
+ * Think of this as A+ = (AT A)^-1 AT, except that the pseudo-inverse
+ * is defined even if AT*A isn't invertible.
+ *
+ * Properties:
+ *
+ * A A+ A = A
+ * A+ A A+ = A+
+ * (A A+)^T = A A+
+ * (A+ A)^T = A+ A
+ *
+ * assuming threshold==0 and no rounding errors.
+ */
+
+void
+SUFFIX(go_matrix_pseudo_inverse) (CONSTMATRIX A, int m, int n,
+ DOUBLE threshold,
+ MATRIX B)
+{
+ QMATRIX V;
+ QMATRIX R;
+ QMATRIX R2 = NULL;
+ QMATRIX R2inv = NULL;
+ gboolean *null_dim;
+ QUAD *x = NULL;
+ gboolean has_result;
+ void *state = SUFFIX(go_quad_start) ();
+ int i, j, i2, j2, rank;
+ DOUBLE emax;
+ gboolean debug = FALSE;
+
+ ALLOC_MATRIX2 (V, m, n, QUAD);
+ ALLOC_MATRIX2 (R, n, n, QUAD);
+ null_dim = g_new0 (gboolean, n);
+
+ has_result = SUFFIX(QRH)(A, FALSE, V, R, m, n, NULL);
+ if (!has_result)
+ goto null_result;
+
+ if (debug) {
+ g_printerr ("R:\n");
+ PRINT_QMATRIX (R, n, n);
+ }
+
+ emax = 0;
+ for (i = 0; i < n; i++) {
+ DOUBLE abs_e = SUFFIX(fabs) (SUFFIX(go_quad_value)(&R[i][i]));
+ emax = MAX (emax, abs_e);
+ }
+
+ rank = 0;
+ for (i = 0; i < n; i++) {
+ DOUBLE abs_e = SUFFIX(fabs) (SUFFIX(go_quad_value)(&R[i][i]));
+ null_dim[i] = abs_e <= emax * threshold;
+ if (!null_dim[i])
+ rank++;
+ }
+ if (rank == 0)
+ goto null_result;
+
+ /*
+ * R2 is R after removing columns and rows that were deemed
+ * null dimensions.
+ */
+ ALLOC_MATRIX2 (R2, rank, rank, QUAD);
+ for (i = i2 = 0; i < n; i++) {
+ if (null_dim[i])
+ continue;
+ for (j = j2 = 0; j < n; j++) {
+ if (null_dim[j])
+ continue;
+ R2[i2][j2] = R[i][j];
+ j2++;
+ }
+ i2++;
+ }
+ if (debug) {
+ g_printerr ("R2:\n");
+ PRINT_QMATRIX (R2, rank, rank);
+ }
+
+ /* R2inv := R2^-1 */
+ ALLOC_MATRIX2 (R2inv, rank, rank, QUAD);
+ x = g_new (QUAD, rank);
+ for (j = 0; j < rank; j++) {
+ for (i = 0; i < rank; i++)
+ SUFFIX(go_quad_init)(&x[i], i == j ? 1 : 0);
+
+ if (SUFFIX(back_solve) (R2, x, x, rank))
+ goto cannot_happen;
+
+ for (i = 0; i < rank; i++)
+ R2inv[i][j] = x[i];
+ }
+ g_free (x);
+ x = NULL;
+ if (debug) {
+ g_printerr ("R2inv:\n");
+ PRINT_QMATRIX (R2inv, rank, rank);
+ }
+
+ /* R := R2inv with null columns/rows inserted back in */
+ for (i = i2 = 0; i < n; i++) {
+ if (null_dim[i]) {
+ for (j = 0; j < n; j++)
+ SUFFIX(go_quad_init)(&R[i][j], 0.0);
+ } else {
+ for (j = j2 = 0; j < n; j++) {
+ if (null_dim[j])
+ SUFFIX(go_quad_init)(&R[i][j], 0.0);
+ else {
+ R[i][j] = R2inv[i2][j2];
+ j2++;
+ }
+ }
+ i2++;
+ }
+ }
+ if (debug) {
+ g_printerr ("R2-inv-expand:\n");
+ PRINT_QMATRIX (R, n, n);
+ }
+
+ /* B := R Q^T */
+ x = g_new (QUAD, m);
+ for (j = 0; j < m; j++) {
+ QUAD acc;
+
+ /* Compute Q^T e_j. */
+ for (i = 0; i < m; i++)
+ SUFFIX(go_quad_init)(&x[i], i == j ? 1 : 0);
+ SUFFIX(multiply_Qt) (x, V, m, n);
+
+ SUFFIX(go_quad_init) (&acc, 0);
+ for (i = 0; i < n; i++) {
+ QUAD acc;
+ int k;
+
+ SUFFIX(go_quad_init) (&acc, 0);
+ for (k = 0; k < n /* Only n */; k++) {
+ QUAD p;
+ SUFFIX(go_quad_mul) (&p, &R[i][k], &x[k]);
+ SUFFIX(go_quad_add) (&acc, &acc, &p);
+ }
+
+ B[i][j] = SUFFIX(go_quad_value)(&acc);
+ }
+ }
+
+out:
+ FREE_MATRIX (V, m, n);
+ FREE_MATRIX (R, n, n);
+ if (R2) FREE_MATRIX (R2, rank, rank);
+ if (R2inv) FREE_MATRIX (R2inv, rank, rank);
+ g_free (null_dim);
+ g_free (x);
+
+ SUFFIX(go_quad_end) (state);
+ return;
+
+cannot_happen:
+ g_printerr ("This cannot happen\n");
+
+null_result:
+ for (i = 0; i < n; i++)
+ for (j = 0; j < m; j++)
+ B[i][j] = 0;
+ goto out;
+}
diff --git a/goffice/math/go-regression.h b/goffice/math/go-regression.h
index 2192efe..36811b2 100644
--- a/goffice/math/go-regression.h
+++ b/goffice/math/go-regression.h
@@ -106,6 +106,10 @@ GORegressionResult go_non_linear_regression (GORegressionFunction f,
gboolean go_matrix_invert (double **A, int n);
double go_matrix_determinant (double *const *const A, int n);
+void go_matrix_pseudo_inverse (double *const * const A, int m, int n,
+ double threshold,
+ double **B);
+
GORegressionResult go_linear_solve (double *const *const A,
const double *b,
int n,
@@ -190,6 +194,10 @@ GORegressionResult go_non_linear_regressionl (GORegressionFunctionl f,
gboolean go_matrix_invertl (long double **A, int n);
long double go_matrix_determinantl (long double *const * const A, int n);
+void go_matrix_pseudo_inversel (long double *const * const A, int m, int n,
+ long double threshold,
+ long double **B);
+
GORegressionResult go_linear_solvel (long double *const *const A,
const long double *b,
int n,
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]