lib/linreg/sweep.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2005, 2009 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 /*
  48   The matrix A will be overwritten. In ordinary uses of the sweep
  49   operator, A will be the matrix
  50
  51    __       __
  52   |X'X    X'Y|
  53   |          |
  54   |Y'X    Y'Y|
  55    --        --
  56
  57    X refers to the design matrix and Y to the vector of dependent
  58    observations. reg_sweep sweeps on the diagonal elements of
  59    X'X.
  60
  61    The matrix A is assumed to be symmetric, so the sweep operation is
  62    performed only for the upper triangle of A.
  63
  64    LAST_COL is considered to be the final column in the augmented matrix,
  65    that is, the column to the right of the '=' sign of the system.
  66  */
  67
  68 int
  69 reg_sweep (gsl_matrix * A, int last_col)
  70 {
  71   double sweep_element;
  72   double tmp;
  73   int i;
  74   int j;
  75   int k;
  76   int row_i;
  77   int row_k;
  78   int row_j;
  79   int *ordered_cols;
  80   gsl_matrix *B;
  81
  82   if (A != NULL)
  83     {
  84       if (A->size1 == A->size2)
  85         {
  86           ordered_cols = malloc (A->size1 * sizeof (*ordered_cols));
  87           for (i = 0; i < last_col; i++)
  88             {
  89               ordered_cols[i] = i;
  90             }
  91           for (i = last_col + 1; i < A->size1; i++)
  92             {
  93               ordered_cols[i - 1] = i;
  94             }
  95           ordered_cols[A->size1 - 1] = last_col;
  96           B = gsl_matrix_alloc (A->size1, A->size2);
  97           for (k = 0; k < (A->size1 - 1); k++)
  98             {
  99               row_k = ordered_cols[k];
 100               sweep_element = gsl_matrix_get (A, row_k, row_k);
 101               if (fabs (sweep_element) > GSL_DBL_MIN)
 102                 {
 103                   tmp = -1.0 / sweep_element;
 104                   gsl_matrix_set (B, row_k, row_k, tmp);
 105                   /*
 106                      Rows before current row k.
 107                    */
 108                   for (i = 0; i < k; i++)
 109                     {
 110                       row_i = ordered_cols[i];
 111                       for (j = i; j < A->size2; j++)
 112                         {
 113                           row_j = ordered_cols[j];
 114                           /*
 115                              Use only the upper triangle of A.
 116                            */
 117                           if (row_j < row_k)
 118                             {
 119                               tmp = gsl_matrix_get (A, row_i, row_j) -
 120                                 gsl_matrix_get (A, row_i, row_k)
 121                                 * gsl_matrix_get (A, row_j, row_k) / sweep_element;
 122                               gsl_matrix_set (B, row_i, row_j, tmp);
 123                             }
 124                           else if (row_j > row_k)
 125                             {
 126                               tmp = gsl_matrix_get (A, row_i, row_j) -
 127                                 gsl_matrix_get (A, row_i, row_k)
 128                                 * gsl_matrix_get (A, row_k, row_j) / sweep_element;
 129                               gsl_matrix_set (B, row_i, row_j, tmp);
 130                             }
 131                           else
 132                             {
 133                               tmp = gsl_matrix_get (A, row_i, row_k) / sweep_element;
 134                               gsl_matrix_set (B, row_i, row_j, tmp);
 135                             }
 136                         }
 137                     }
 138                   /*
 139                      Current row k.
 140                    */
 141                   for (j = k + 1; j < A->size1; j++)
 142                     {
 143                       row_j = ordered_cols[j];
 144                       tmp = gsl_matrix_get (A, row_k, row_j) / sweep_element;
 145                       gsl_matrix_set (B, row_k, row_j, tmp);
 146                     }
 147                   /*
 148                      Rows after the current row k.
 149                    */
 150                   for (i = k + 1; i < A->size1; i++)
 151                     {
 152                       row_i = ordered_cols[i];
 153                       for (j = i; j < A->size2; j++)
 154                         {
 155                           row_j = ordered_cols[j];
 156                           tmp = gsl_matrix_get (A, row_i, row_j) -
 157                             gsl_matrix_get (A, row_k, row_i)
 158                             * gsl_matrix_get (A, row_k, row_j) / sweep_element;
 159                           gsl_matrix_set (B, row_i, row_j, tmp);
 160                         }
 161                     }
 162                 }
 163               for (i = 0; i < A->size1; i++)
 164                 for (j = i; j < A->size2; j++)
 165                   {
 166                     row_i = ordered_cols[i];
 167                     row_j = ordered_cols[j];
 168                     gsl_matrix_set (A, row_i, row_j, gsl_matrix_get (B, row_i, row_j));
 169                   }
 170             }
 171           gsl_matrix_free (B);
 172           free (ordered_cols);
 173           return GSL_SUCCESS;
 174         }
 175       return GSL_ENOTSQR;
 176     }
 177   return GSL_EFAULT;
 178 }