[R] Solving Matrices

[etb]

"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
}