[R] Row-Echelon Form
John Fox
jfox at mcmaster.ca
Tue Sep 4 04:21:45 CEST 2007
Dear Peter,
Some time ago, I posted a function RREF to r-help that computes the reduced
row-echelon form of a matrix. Just last week, Scott Hyde posted a revised
version of this function to r-help (see
https://stat.ethz.ch/pipermail/r-help/2007-September/139923.html). My
original function didn't work correctly in some cases, and in response to
Scott's message, I tried to post a newer version to r-help, but it doesn't
seem to have gotten through (perhaps because I attached the function in a
file). I append the function below, along with some other simple
linear-algebra functions.
I hope this helps,
John
------------- snip -------------
# Last modified 9 July 2007 by J. Fox
GaussianElimination <- function(A, B, tol=sqrt(.Machine$double.eps),
verbose=FALSE, fractions=FALSE){
# A: coefficient matrix
# B: right-hand side vector or matrix
# tol: tolerance for checking for 0 pivot
# verbose: if TRUE, print intermediate steps
# fractions: try to express nonintegers as rational numbers
# If B is absent returns the reduced row-echelon form of A.
# If B is present, reduces A to RREF carrying B along.
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)
if (!missing(B)){
B <- as.matrix(B)
if (!(nrow(B) == nrow(A)) || !is.numeric(B))
stop("argument must be numeric and must match the number of row of
A")
A <- cbind(A, B)
}
i <- j <- 1
while (i <= n && j <= m){
while (j <= m){
currentColumn <- A[,j]
currentColumn[1:n < i] <- 0
# find maximum pivot in current column at or below current row
which <- which.max(abs(currentColumn))
pivot <- currentColumn[which]
if (abs(pivot) <= tol) { # check for 0 pivot
j <- j + 1
next
}
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[,j], row) # sweep
A[i,] <- row # restore current row
if (verbose) if (fractions) print(fractions(A))
else print(round(A, round(abs(log(tol,10)))))
j <- j + 1
break
}
i <- i + 1
}
# 0 rows to bottom
zeros <- which(apply(A[,1:m], 1, function(x) max(abs(x)) <= tol))
if (length(zeros) > 0){
zeroRows <- A[zeros,]
A <- A[-zeros,]
A <- rbind(A, zeroRows)
}
rownames(A) <- NULL
if (fractions) fractions (A) else round(A, round(abs(log(tol, 10))))
}
matrixInverse <- function(X, tol=sqrt(.Machine$double.eps), ...){
# returns the inverse of nonsingular X
if ((!is.matrix(X)) || (nrow(X) != ncol(X)) || (!is.numeric(X)))
stop("X must be a square numeric matrix")
n <- nrow(X)
X <- GaussianElimination(X, diag(n), tol=tol, ...) # append identity
matrix
# check for 0 rows in the RREF of X:
if (any(apply(abs(X[,1:n]) <= sqrt(.Machine$double.eps), 1, all)))
stop ("X is numerically singular")
X[,(n + 1):(2*n)] # return inverse
}
RREF <- function(X, ...) GaussianElimination(X, ...)
# returns the reduced row-echelon form of X
Ginv <- function(A, tol=sqrt(.Machine$double.eps), verbose=FALSE,
fractions=FALSE){
# return an arbitrary generalized inverse of the matrix A
# A: a matrix
# tol: tolerance for checking for 0 pivot
# verbose: if TRUE, print intermediate steps
# fractions: try to express nonintegers as rational numbers
m <- nrow(A)
n <- ncol(A)
B <- GaussianElimination(A, diag(m), tol=tol, verbose=verbose,
fractions=fractions)
L <- B[,-(1:n)]
AR <- B[,1:n]
C <- GaussianElimination(t(AR), diag(n), tol=tol, verbose=verbose,
fractions=fractions)
R <- t(C[,-(1:m)])
AC <- t(C[,1:m])
ginv <- R %*% t(AC) %*% L
if (fractions) fractions (ginv) else round(ginv, round(abs(log(tol,
10))))
}
cholesky <- function(X, tol=sqrt(.Machine$double.eps)){
# returns the Cholesky square root of the nonsingular, symmetric matrix
X
# tol: tolerance for checking for 0 pivot
# algorithm from Kennedy & Gentle (1980)
if (!is.numeric(X)) stop("argument is not numeric")
if (!is.matrix(X)) stop("argument is not a matrix")
n <- nrow(X)
if (ncol(X) != n) stop("matrix is not square")
if (max(abs(X - t(X))) > tol) stop("matrix is not symmetric")
D <- rep(0, n)
L <- diag(n)
i <- 2:n
D[1] <- X[1, 1]
if (abs(D[1]) < tol) stop("matrix is numerically singular")
L[i, 1] <- X[i, 1]/D[1]
for (j in 2:(n - 1)){
k <- 1:(j - 1)
D[j] <- X[j, j] - sum((L[j, k]^2) * D[k])
if (abs(D[j]) < tol) stop("matrix is numerically singular")
i <- (j + 1):n
L[i, j] <- (X[i, j] -
colSums(L[j, k] * t(L[i, k, drop=FALSE]) *
D[k]))/D[j]
}
k <- 1:(n - 1)
D[n] <- X[n, n] - sum((L[n, k]^2) * D[k])
if (abs(D[n]) < tol) stop("matrix is numerically singular")
L %*% diag(sqrt(D))
}
--------------------------------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
--------------------------------
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Peter Danenberg
> Sent: Monday, September 03, 2007 2:17 PM
> To: r-help at stat.math.ethz.ch
> Subject: [R] Row-Echelon Form
>
> I was looking for an R-package that would reduce matrices to
> row-echelon form, but Google was not my friend; any leads?
>
> If not, I wonder if the problem could be expressed in terms
> of constraint satisfaction...
>
> ______________________________________________
> R-help at 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.
>
More information about the R-help
mailing list