atom feed2 messages in org.r-project.r-help[R] row echelon form
FromSent OnAttachments
Scott HydeAug 31, 2007 8:56 pm 
Peng, CSep 6, 2010 10:31 pm 
Subject:[R] row echelon form
From:Scott Hyde (hyd@byuh.edu)
Date:Aug 31, 2007 8:56:39 pm
List:org.r-project.r-help

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