lib/linreg/sweep.c

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