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
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
