lib/linreg/sweep.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2005 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 "sweep.h"
  44
  45 /*
  46   The matrix A will be overwritten. In ordinary uses of the sweep
  47   operator, A will be the matrix
  48
  49    __       __
  50   |X'X    X'Y|
  51   |          |
  52   |Y'X    Y'Y|
  53    --        --
  54
  55    X refers to the design matrix and Y to the vector of dependent
  56    observations. reg_sweep sweeps on the diagonal elements of
  57    X'X.
  58
  59    The matrix A is assumed to be symmetric, so the sweep operation is
  60    performed only for the upper triangle of A.
  61  */
  62
  63 int
  64 reg_sweep (gsl_matrix * A)
  65 {
  66   double sweep_element;
  67   double tmp;
  68   int i;
  69   int j;
  70   int k;
  71   gsl_matrix *B;
  72
  73   if (A != NULL)
  74     {
  75       if (A->size1 == A->size2)
  76         {
  77           B = gsl_matrix_alloc (A->size1, A->size2);
  78           for (k = 0; k < (A->size1 - 1); k++)
  79             {
  80               sweep_element = gsl_matrix_get (A, k, k);
  81               if (fabs (sweep_element) > GSL_DBL_MIN)
  82                 {
  83                   tmp = -1.0 / sweep_element;
  84                   gsl_matrix_set (B, k, k, tmp);
  85                   /*
  86                      Rows before current row k.
  87                    */
  88                   for (i = 0; i < k; i++)
  89                     {
  90                       for (j = i; j < A->size2; j++)
  91                         {
  92                           /*
  93                              Use only the upper triangle of A.
  94                            */
  95                           if (j < k)
  96                             {
  97                               tmp = gsl_matrix_get (A, i, j) -
  98                                 gsl_matrix_get (A, i, k)
  99                                 * gsl_matrix_get (A, j, k) / sweep_element;
 100                               gsl_matrix_set (B, i, j, tmp);
 101                             }
 102                           else if (j > k)
 103                             {
 104                               tmp = gsl_matrix_get (A, i, j) -
 105                                 gsl_matrix_get (A, i, k)
 106                                 * gsl_matrix_get (A, k, j) / sweep_element;
 107                               gsl_matrix_set (B, i, j, tmp);
 108                             }
 109                           else
 110                             {
 111                               tmp = gsl_matrix_get (A, i, k) / sweep_element;
 112                               gsl_matrix_set (B, i, j, tmp);
 113                             }
 114                         }
 115                     }
 116                   /*
 117                      Current row k.
 118                    */
 119                   for (j = k + 1; j < A->size1; j++)
 120                     {
 121                       tmp = gsl_matrix_get (A, k, j) / sweep_element;
 122                       gsl_matrix_set (B, k, j, tmp);
 123                     }
 124                   /*
 125                      Rows after the current row k.
 126                    */
 127                   for (i = k + 1; i < A->size1; i++)
 128                     {
 129                       for (j = i; j < A->size2; j++)
 130                         {
 131                           tmp = gsl_matrix_get (A, i, j) -
 132                             gsl_matrix_get (A, k, i)
 133                             * gsl_matrix_get (A, k, j) / sweep_element;
 134                           gsl_matrix_set (B, i, j, tmp);
 135                         }
 136                     }
 137                 }
 138               for (i = 0; i < A->size1; i++)
 139                 for (j = i; j < A->size2; j++)
 140                   {
 141                     gsl_matrix_set (A, i, j, gsl_matrix_get (B, i, j));
 142                   }
 143             }
 144           gsl_matrix_free (B);
 145           return GSL_SUCCESS;
 146         }
 147       return GSL_ENOTSQR;
 148     }
 149   return GSL_EFAULT;
 150 }