lib/linreg/sweep.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2005, 2009, 2011 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 /*
  18   Find the least-squares estimate of b for the linear model:
  19
  20   Y = Xb + Z
  21
  22   where Y is an n-by-1 column vector, X is an n-by-p matrix of
  23   independent variables, b is a p-by-1 vector of regression coefficients,
  24   and Z is an n-by-1 normally-distributed random vector with independent
  25   identically distributed components with mean 0.
  26
  27   This estimate is found via the sweep operator, which is a modification
  28   of Gauss-Jordan pivoting.
  29
  30
  31   References:
  32
  33   Matrix Computations, third edition. GH Golub and CF Van Loan.
  34   The Johns Hopkins University Press. 1996. ISBN 0-8018-5414-8.
  35
  36   Numerical Analysis for Statisticians. K Lange. Springer. 1999.
  37   ISBN 0-387-94979-8.
  38
  39   Numerical Linear Algebra for Applications in Statistics. JE Gentle.
  40   Springer. 1998. ISBN 0-387-98542-5.
  41 */
  42
  43 #include <config.h>
  44
  45 #include "sweep.h"
  46
  47 #include <assert.h>
  48 /*
  49   The matrix A will be overwritten. In ordinary uses of the sweep
  50   operator, A will be the matrix
  51
  52   __       __
  53   |X'X    X'Y|
  54   |          |
  55   |Y'X    Y'Y|
  56   --        --
  57
  58   X refers to the design matrix and Y to the vector of dependent
  59   observations. reg_sweep sweeps on the diagonal elements of
  60   X'X.
  61
  62   The matrix A is assumed to be symmetric, so the sweep operation is
  63   performed only for the upper triangle of A.
  64
  65   LAST_COL is considered to be the final column in the augmented matrix,
  66   that is, the column to the right of the '=' sign of the system.
  67 */
  68
  69 int
  70 reg_sweep (gsl_matrix * A, int last_col)
  71 {
  72   int i;
  73   int j;
  74   int k;
  75   gsl_matrix *B;
  76
  77   if (A == NULL)
  78     return GSL_EFAULT;
  79
  80   if (A->size1 != A->size2)
  81     return GSL_ENOTSQR;
  82
  83   assert (last_col < A->size1);
  84   gsl_matrix_swap_rows (A, A->size1 - 1, last_col);
  85   gsl_matrix_swap_columns (A, A->size1 - 1 , last_col);
  86
  87   B = gsl_matrix_alloc (A->size1, A->size2);
  88   for (k = 0; k < (A->size1 - 1); k++)
  89     {
  90       const double sweep_element = gsl_matrix_get (A, k, k);
  91       if (fabs (sweep_element) > GSL_DBL_MIN)
  92         {
  93           gsl_matrix_set (B, k, k, -1.0 / sweep_element);
  94           /*
  95             Rows before current row k.
  96           */
  97           for (i = 0; i < k; i++)
  98             {
  99               for (j = i; j < A->size2; j++)
 100                 {
 101                   /* Use only the upper triangle of A. */
 102                   double tmp;
 103                   if (j < k)
 104                     {
 105                       tmp = gsl_matrix_get (A, i, j) -
 106                         gsl_matrix_get (A, i, k)
 107                         * gsl_matrix_get (A, j, k) / sweep_element;
 108                     }
 109                   else if (j > k)
 110                     {
 111                       tmp = gsl_matrix_get (A, i, j) -
 112                         gsl_matrix_get (A, i, k)
 113                         * gsl_matrix_get (A, k, j) / sweep_element;
 114                     }
 115                   else
 116                     {
 117                       tmp = gsl_matrix_get (A, i, k) / sweep_element;
 118                     }
 119                   gsl_matrix_set (B, i, j, tmp);
 120                 }
 121             }
 122           /*
 123             Current row k.
 124           */
 125           for (j = k + 1; j < A->size1; j++)
 126             {
 127               double tmp = gsl_matrix_get (A, k, j) / sweep_element;
 128               gsl_matrix_set (B, k, j, tmp);
 129             }
 130           /*
 131             Rows after the current row k.
 132           */
 133           for (i = k + 1; i < A->size1; i++)
 134             {
 135               for (j = i; j < A->size2; j++)
 136                 {
 137                   double tmp = gsl_matrix_get (A, i, j) -
 138                     gsl_matrix_get (A, k, i)
 139                     * gsl_matrix_get (A, k, j) / sweep_element;
 140                   gsl_matrix_set (B, i, j, tmp);
 141                 }
 142             }
 143         }
 144       gsl_matrix_memcpy (A, B);
 145     }
 146   gsl_matrix_free (B);
 147
 148   gsl_matrix_swap_columns (A, A->size1 - 1 , last_col);
 149   gsl_matrix_swap_rows (A, A->size1 - 1, last_col);
 150
 151   return GSL_SUCCESS;
 152 }