2010-11-20 59 views
74

Me preguntaba ¿cuál es la forma recomendada de calcular el inverso de una matriz?inversa de la matriz en R

Las formas que encontré parecen no ser satisfactorias. Por ejemplo,

> c=rbind(c(1, -1/4), c(-1/4, 1)) 
> c 
     [,1] [,2] 
[1,] 1.00 -0.25 
[2,] -0.25 1.00 
> inv(c) 
Error: could not find function "inv" 
> solve(c)  
      [,1]  [,2] 
[1,] 1.0666667 0.2666667 
[2,] 0.2666667 1.0666667 
> solve(c)*c 
      [,1]  [,2] 
[1,] 1.06666667 -0.06666667 
[2,] -0.06666667 1.06666667 
> qr.solve(c)*c 
      [,1]  [,2] 
[1,] 1.06666667 -0.06666667 
[2,] -0.06666667 1.06666667 

¡Gracias!

+1

Un consejo general: evitar dar objetos (como las matrices) un nombre que ya se utiliza (en este caso 'C'). – Qaswed

Respuesta

118

solve(c) da la inversa correcta. El problema con su código es que está utilizando el operador incorrecto para la multiplicación de la matriz. Debe usar solve(c) %*% c para invocar la multiplicación de matrices en R.

R realiza la multiplicación elemento por elemento cuando invoca solve(c) * c.

19

Puede utilizar la función ginv() (Moore-Penrose inversa generalizada) en el paquete MASA

+0

@xeon no estoy seguro de cómo se puede perder - ver p. 60 del Manual para el paquete mencionado en mi respuesta anterior – doug

+0

Gracias por su respuesta. Obtuve este error al ejecutar la función fem() desde el paquete FisherEM. Ejecutando Mavericks Mac OS X. –

4

Tenga en cuenta que si se preocupan por la velocidad y no es necesario preocuparse por las singularidades, solve() se debe preferir a ginv() porque es mucho más rápido, como se puede comprobar:

require(MASS) 
mat <- matrix(rnorm(1e6),nrow=1e3,ncol=1e3) 

t0 <- proc.time() 
inv0 <- ginv(mat) 
proc.time() - t0 

t1 <- proc.time() 
inv1 <- solve(mat) 
proc.time() - t1 
0

En notación matricial que hace una gran diferencia al operador "* "y el operador" %*% ". El primero tiene un elemento de multiplicación por elemento, el segundo es la fórmula correcta para la multiplicación de la matriz. Lo hou debería haber hecho es:

c = rbind(c(1, -1/4), c(-1/4, 1)) 

solve(c) %*% c 
+5

¿En qué se diferencia su respuesta de la aceptada? –