[R] row echelon form

Scott Hyde hydes at byuh.edu
Sat Sep 1 05:57:02 CEST 2007

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

rref <- function(A, tol=sqrt(.Machine$double.eps),verbose=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.

Any suggestions of problems that anyone sees?


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

More information about the R-help mailing list