Reduced row echelon form

Question

Is there a function in R that produces the reduced row echelon form of a matrix?. This reference says there isn't. Do you agree?

Solution 1

The pracma package also contains an implementation. See pracma::rref.

Solution 2

I don't have enough rep to comment, but the function given above by soldier.moth in the accepted answer [edit 2018: no longer the accepted answer] is buggy - it doesn't handle matrices where the RREF solution has zeroes on its main diagonal. Try e.g.

m<-matrix(c(1,0,1,0,0,2),byrow=TRUE,nrow=2) rref(m)

and note that the output is not in RREF.

I think I have it working, but you may want to check outputs for yourself:

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
  # Modified by Geoffrey Brent 2014-12-17 to fix a bug
  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)
  x.position<-1
  y.position<-1
  # change loop:
  while((x.position<=m) & (y.position<=n)){
    col <- A[,x.position]
    col[1:n < y.position] <- 0
    # find maximum pivot in current column at or below current row
    which <- which.max(abs(col))
    pivot <- col[which]
    if (abs(pivot) <= tol) x.position<-x.position+1     # check for 0 pivot
    else{
      if (which > y.position) { A[c(y.position,which),]<-A[c(which,y.position),] } # exchange rows
      A[y.position,]<-A[y.position,]/pivot # pivot
      row <-A[y.position,]
      A <- A - outer(A[,x.position],row) # sweep
      A[y.position,]<-row # restore current row
      if (verbose)
        if (fractions) print(fractions(A))
        else print(round(A,round(abs(log(tol,10)))))
      x.position<-x.position+1
      y.position<-y.position+1
    }
  }
  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))))
}

Solution 3

There is also a recent package developed for teaching Linear Algebra (matlib) which both computes the echelon form of a matrix, and shows the steps used along the way.

Example from the reference docs:

library('matlib')
A <- matrix(c(2, 1, -1,-3, -1, 2,-2, 1, 2), 3, 3, byrow=TRUE)
b <- c(8, -11, -3)
echelon(A, b, verbose=TRUE, fractions=TRUE)

Initial matrix:
     [,1] [,2] [,3] [,4]
[1,]   2    1   -1    8 
[2,]  -3   -1    2  -11 
[3,]  -2    1    2   -3 

row: 1 

 exchange rows 1 and 2 
     [,1] [,2] [,3] [,4]
[1,]  -3   -1    2  -11 
[2,]   2    1   -1    8 
[3,]  -2    1    2   -3 

 multiply row 1 by -1/3 
     [,1] [,2] [,3] [,4]
[1,]    1  1/3 -2/3 11/3
[2,]    2    1   -1    8
[3,]   -2    1    2   -3

 multiply row 1 by 2 and subtract from row 2 
     [,1] [,2] [,3] [,4]
[1,]    1  1/3 -2/3 11/3
[2,]    0  1/3  1/3  2/3
[3,]   -2    1    2   -3

 multiply row 1 by 2 and add to row 3 
     [,1] [,2] [,3] [,4]
[1,]    1  1/3 -2/3 11/3
[2,]    0  1/3  1/3  2/3
[3,]    0  5/3  2/3 13/3

row: 2 

 exchange rows 2 and 3 
     [,1] [,2] [,3] [,4]
[1,]    1  1/3 -2/3 11/3
[2,]    0  5/3  2/3 13/3
[3,]    0  1/3  1/3  2/3

 multiply row 2 by 3/5 
     [,1] [,2] [,3] [,4]
[1,]    1  1/3 -2/3 11/3
[2,]    0    1  2/5 13/5
[3,]    0  1/3  1/3  2/3

 multiply row 2 by 1/3 and subtract from row 1 
     [,1] [,2] [,3] [,4]
[1,]    1    0 -4/5 14/5
[2,]    0    1  2/5 13/5
[3,]    0  1/3  1/3  2/3

 multiply row 2 by 1/3 and subtract from row 3 
     [,1] [,2] [,3] [,4]
[1,]    1    0 -4/5 14/5
[2,]    0    1  2/5 13/5
[3,]    0    0  1/5 -1/5

row: 3 

 multiply row 3 by 5 
     [,1] [,2] [,3] [,4]
[1,]    1    0 -4/5 14/5
[2,]    0    1  2/5 13/5
[3,]    0    0    1   -1

 multiply row 3 by 4/5 and add to row 1 
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    2
[2,]    0    1  2/5 13/5
[3,]    0    0    1   -1

 multiply row 3 by 2/5 and subtract from row 2 
     [,1] [,2] [,3] [,4]
[1,]  1    0    0    2  
[2,]  0    1    0    3  
[3,]  0    0    1   -1  

Solution 4

Doesn't look like there is one built in but I found this rref function on this page.

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

View more solutions

Comments

[László Papp]

This does not provide an answer to the question. To critique or request clarification from an author, leave a comment below their post - you can always comment on your own posts, and once you have sufficient reputation you will be able to comment on any post.

[Geoffrey Brent]

Sorry, I'm new here and may have missed something, but: the "accepted answer" provided by soldier.moth above is buggy (as I discovered the hard way when I tried to use it myself!) so I thought it was important to flag that. I don't have enough rep to comment on soldier.moth's answer directly, so I created a new answer - what should I have done here?

[László Papp]

Gain enough reputation for commenting first?

[Geoffrey Brent]

shrug I really didn't think anybody would be bothered by me pointing out a non-obvious bug that had gone four years without correction, and given that the buggy "solution" posted above was accepted as an answer, I'm perplexed as to why a fixed version of the same code is considered less so.

[Kevin]

@GeoffreyBrent, Welcome to Stackoverflow. Where rules supersede common sense.