2010-06-27 18 views

Respuesta

8

El paquete pracma también contiene una implementación. Ver pracma :: rref.

4

Parece que no hay uno incorporado, pero encontré esta función rref en la página this.

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

No tengo suficientes representantes para comentar, pero la función dada anteriormente en la respuesta aceptada es incorrecta: no maneja matrices donde la solución RREF tiene ceros en su diagonal principal. Pruebe, por ejemplo,

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

y tenga en cuenta que la salida no está en RREF .

creo que tengo que trabajar, pero es posible que desee comprobar las salidas por sí mismo:

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

Esto no proporciona una respuesta a la pregunta. Para criticar o solicitar aclaraciones de un autor, deje un comentario debajo de su publicación; siempre puede comentar sus propias publicaciones, y una vez que tenga suficiente [reputación] (http://stackoverflow.com/help/whats-reputation) lo hará poder [comentar cualquier publicación] (http://stackoverflow.com/help/privileges/comment). – lpapp

+1

Lo siento, soy nuevo aquí y es posible que haya pasado algo por alto, pero: la "respuesta aceptada" proporcionada por soldier.moth es incorrecta (como descubrí de la peor manera cuando traté de usarla yo mismo) así que pensé que era importante para marcar eso. No tengo suficientes representantes para comentar la respuesta de soldier.moth directamente, así que creé una nueva respuesta. ¿Qué debería haber hecho aquí? –

+0

¿Tiene suficiente reputación para comentar primero? – lpapp

5

También hay un paquete de reciente desarrollado para enseñanza Álgebra Lineal (matlib), que tanto calcula la forma escalonada de una matriz, y muestra los pasos utilizados en el camino.

Ejemplo de la 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