[R] Solving Matrices

[John Fox]

Dear Elizabeth,

If this is the rowEchelonForm function that I posted to r-help some
time ago, as I mentioned in a subsequent post, the tolerance was set
too low; if you set the tol argument (better, the default) to a larger
number than the machine precision, then the function should work more
reliably.

I hope this helps,
 John

On 16 Apr 2004 11:09:33 -0500
 etb <lizzy at noradd.org> wrote:
> "Richard" writes:
> 
> [ ... ]
> >
> > So here are three different functions all saying much the same
> thing
> > about x: it is numerically close to a matrix which does NOT span
> the
> > whole of R**4.
> 
> Theorem 4 (in my text):
> 
> 	Let A be an m x n matrix. Then the following statements are
> 	logically equivalent. That is, for a particular A, either they
> 	are all true statements or they are all false:
> 
> 	  a) For each b in R**m, the equation Ax = b has a solution.
> 	  b) The columns of A span R**m.
> 	  c) A has a pivot position in every row.
> 
> The back of the book is in agreement with you but I can manipulate
> the
> matrix to produce (c) above which should mean that (b) is correct:
> 
> > x <- matrix(data=c(7,-5,6,-7,2,-3,10,9,-5,4,-2,2,8,-9,7,15),nrow=4)
> > rowEchelonForm(x)
>      [,1] [,2] [,3] [,4]
> [1,]    1    0    0    0
> [2,]    0    1    0    0
> [3,]    0    0    1    0
> [4,]    0    0    0    1
> 
> Likewise, I can produce the above using the swap(), scale(), gauss()
> and bgauss() functions that the author of the text originally wrote
> for maple, Mathematica, etc. but from all I can tell, the row echelon
> form of matrix x meets the criteria of theorem 4 above, yet the
> columns
> of x do not span R**4.
> 
> Using the gauss()[1] function:
> 
>      [,1]      [,2]       [,3]          [,4]
> [1,]    7  2.000000 -5.0000000  8.000000e+00
> [2,]    0 -1.571429  0.4285714 -3.285714e+00
> [3,]    0  0.000000  4.5454545 -1.718182e+01
> [4,]    0  0.000000  0.0000000 -5.035972e-15
> 
> Again, this has a pivot column in each row although [4,4] seems a bit
> suspicious.
> 
> Elizabeth
> 
> [1]:
> 
> gauss <- function(A, r, v = c()) {
>   m <- dim(A)[1]
>   n <- dim(A)[2]
>   if(length(r) > 1 | r < 1 | r > m)
>     stop("The second entry in gauss(...) must be the pivot row
> index.")
>   if(length(v) == 0) {
>     if(r+1 > m)
>       v <- m
>     else
>       v <- (r+1):m
>   }
>   for(j in v) {
>     if(j == r)
>       stop("Pivot ow cannot change itself.")
>   }
> 
>   col <- 1
>   while(A[r, col] == 0 && col < n)
>     col <- col + 1
>   if(col == n && A[r, col] == 0)
>     stop("Row r has only zeroes and cannot be a pivot row.")
>   for(j in v) {
>     multiplier <- -A[j, col]/A[r, col]
>     A[j,] <- A[j,] + multiplier * A[r,]
>   }
>   A
> }
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html


--------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/