[R] row echelon form

[Scott Hyde]

Hi everyone,

I am looking to use R as a MATLAB replacement for linear algebra.
I've done a fairly good job for finding replacements for most of the
functions I'm interested in, I
John Fox wrote a program for implementing the reduced row echelon form
of a matrix (by doing the Gauss-Jordan elimination).  I modified it a
bit:

rref <- function(A, tol=sqrt(.Machine$double.eps),verbose=FALSE,
                 fractions=FALSE){
  ## A: coefficient matrix
  ## tol: tolerance for checking for 0 pivot
  ## verbose: if TRUE, print intermediate steps
  ## fractions: try to express nonintegers as rational numbers
  ## Written by John Fox
  if (fractions) {
    mass <- require(MASS)
    if (!mass) stop("fractions=TRUE needs MASS package")
  }
  if ((!is.matrix(A)) || (!is.numeric(A)))
    stop("argument must be a numeric matrix")
  n <- nrow(A)
  m <- ncol(A)
  for (i in 1:min(c(m, n))){
    col <- A[,i]
    col[1:n < i] <- 0
    # find maximum pivot in current column at or below current row
    which <- which.max(abs(col))
    pivot <- A[which, i]
    if (abs(pivot) <= tol) next     # check for 0 pivot
    if (which > i) A[c(i, which),] <- A[c(which, i),]  # exchange rows
    A[i,] <- A[i,]/pivot            # pivot
    row <- A[i,]
    A <- A - outer(A[,i], row)      # sweep
    A[i,] <- row                    # restore current row
    if (verbose)
      if (fractions) print(fractions(A))
      else print(round(A,round(abs(log(tol,10)))))
  }
  for (i in 1:n)
    if (max(abs(A[i,1:m])) <= tol)
      A[c(i,n),] <- A[c(n,i),] # 0 rows to bottom
  if (fractions) fractions (A)
  else round(A, round(abs(log(tol,10))))
}

I've found its useful for my students.

I wrote a program for implementing row echelon form, which is only
concerned with the getting the lower triangular part to be zero, along
with the first non-zero entry in each non-zero row being a one.  Here
is what I wrote:

ref <- function(A) {
  ## Implements gaussian elimination using the QR factorization.
  ## Note: ref form is not unique.
  ## written by S. K. Hyde
  ##

  qrA <- qr(A)
  r <- qrA$rank  ##computed rank of the matrix
  R <- qr.R(qrA)

  if (r < nrow(R))
    R[-(1:r),] <- 0 #zero out rows that should be zero (according to the rank).

  #Get ones as the first nonzero entry in each row and return result.
  diag(c(1/diag(R)[1:r],rep(1,nrow(R)-r)))%*%R
}

Any suggestions of problems that anyone sees?

-Scott



-- 
*****************************************************************
Scott K. Hyde
Assistant Professor of Statistics and Mathematics
Brigham Young University -- Hawaii