Matrices

Hay varias maneras de definir una matriz en R. Si es una matriz pequeña podemos utilizar la siguiente sintaxis:

A=matrix(nrow=3,ncol=3, c(1,2,3,4,5,6,7,8,9),byrow=TRUE); 
A
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
## [3,]    7    8    9
A_1=matrix(nrow=3,ncol=3, c(1,2,3,4,5,6,7,8,9),byrow=FALSE)
A_1
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9
B_1 <- matrix(c(1,2,3,4,5,6,7,8,9),byrow=TRUE,nrow=3,ncol=3); 
B_1
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6
## [3,]    7    8    9
B_2 <- matrix(c(1,2,3,4,5,6,7,8,9),byrow=FALSE,nrow=3,ncol=3)
B_2
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9

Si disponemos de varios vectores de la misma longitud que queremos utilizar como filas (o columnas) de una matriz, podemos utilizar la función rbind() para unirlos por filas, o la función cbind() para unirlos por columnas, como vemos en el siguiente ejemplo:

vector1 = c(1,2,3,4)
vector2 = c(5,6,7,8)
vector3 = c(9,10,11,12)
M1 = cbind(vector1,vector2,vector3) # Unimos por columnas
M1  
##      vector1 vector2 vector3
## [1,]       1       5       9
## [2,]       2       6      10
## [3,]       3       7      11
## [4,]       4       8      12

Se pueden seleccionar partes de una matriz utilizando los índices de posición [fila, columna] entre corchetes. El siguiente ejemplo ilustra la forma de hacerlo:

A[2,3]   # Se selecciona el valor de la fila 2, columna 3
## [1] 6
A[2,]    # Se selecciona la fila 2 completa
## [1] 4 5 6
A[,3]    # Se selecciona la columna 3 completa
## [1] 3 6 9
A[1,2:3] # Se seleccionan el segundo y tercer valor de la fila 1
## [1] 2 3

Operaciones con matrices

La función diag() extrae la diagonal principal de una matriz:

diag(A)
## [1] 1 5 9

También permite crear matrices diagonales:

diag(c(1,2,3,4))
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    0    2    0    0
## [3,]    0    0    3    0
## [4,]    0    0    0    4

La matriz identidad es muy fácil de crear en R. Por ejemplo, la matriz identidad de dimensión 4 es:

Id4=diag(1,nrow=4)
Id4
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    0    1    0    0
## [3,]    0    0    1    0
## [4,]    0    0    0    1

Hay que tener cierto cuidado con los operadores aritméticos básicos. Si se suman una matriz y una constante, el efecto es que dicha constante se suma a todos los elementos de la matriz. Lo mismo ocurre con la diferencia, la multiplicación y la división:

M = matrix(nrow=2,c(1,2,5,7),byrow=F)
M
##      [,1] [,2]
## [1,]    1    5
## [2,]    2    7
M+2
##      [,1] [,2]
## [1,]    3    7
## [2,]    4    9

Asimismo, si a una matriz se le suma un vector cuya longitud sea igual al número de filas de la matriz, se obtiene como resultado una nueva matriz cuyas columnas son la suma de las columnas de la matriz original más dicho vector.

v=c(3,4)
M+v
##      [,1] [,2]
## [1,]    4    8
## [2,]    6   11

La suma o resta de matrices de la misma dimensión se realiza con los operadores + y -; el producto de matrices (siempre que sean compatibles) se realiza con el símbolo%*%:

M+M
##      [,1] [,2]
## [1,]    2   10
## [2,]    4   14
M-M
##      [,1] [,2]
## [1,]    0    0
## [2,]    0    0
M%*%M
##      [,1] [,2]
## [1,]   11   40
## [2,]   16   59

Una fuente de posibles errores en el cálculo matricial, cuando se utilizan matrices de la misma dimensión, es utilizar los operadores * y / ya que multiplican (o dividen) las matrices término a término:

M*M
##      [,1] [,2]
## [1,]    1   25
## [2,]    4   49
M/M
##      [,1] [,2]
## [1,]    1    1
## [2,]    1    1

Preguntas

  1. Crea una matriz 4x4 que contenga los números del 5 al 20, colocando los números por columna. ¿Cómo accederías al valor en la tercera fila y cuarta columna?

  2. Dado los vectores: v1 <- c(2, 4, 6), v2 <- c(8, 10, 12) Únelos para formar una matriz por columnas. ¿Cómo sería la matriz resultante si los unes por filas?

  3. Utiliza la matriz B, definida como: B <- matrix(11:26, nrow = 4, ncol = 4) ¿Cómo seleccionarías la segunda fila completa y la tercera columna completa?

  4. Crea una matriz diagonal utilizando los números c(2, 4, 6, 8). ¿Qué comando usarías para obtener la matriz identidad de dimensión 4?

  5. Si tienes una matriz M definida como: M <- matrix(c(2, 3, 5, 7), nrow = 2) Suma 2 a cada elemento de la matriz. ¿Cómo cambiaría la matriz después de la suma?

  6. Dada una matriz M = matrix(c(1, 2, 3, 4, 5, 6), nrow = 2), suma el vector v = c(1, 2) a la matriz M.¿Qué operación ocurre en cada columna?

Potencia de una matriz.

R no dispone en su paquete base de una función para calcular la potencia n-ésima de una matriz. No obstante el paquete expm implementa el operador %^% construido con este objetivo. Su uso es muy simple: si queremos calcular la matriz An bastará con utilizar el comando A %^% n. Para que funcione esta instrucción debemos cargar previamente el paquete expm:

library(expm)
## Loading required package: Matrix
## 
## Attaching package: 'expm'
## The following object is masked from 'package:Matrix':
## 
##     expm
A_2 = matrix(1:4,nrow=2)
A_2
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4
A_2 %*% A_2
##      [,1] [,2]
## [1,]    7   15
## [2,]   10   22
A_2 %^% 2 
##      [,1] [,2]
## [1,]    7   15
## [2,]   10   22
A_2 %^% 5
##      [,1] [,2]
## [1,] 1069 2337
## [2,] 1558 3406

Matriz traspuesta

La traspuesta de una matriz se calcula simplemente con la función t()

t(A)
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9

Determinante

El determinante de una matriz cuadrada se calcula mediante la función det():

det(A)
## [1] -9.516197e-16

Preguntas

  1. Multiplica la matriz M definida como: M <- matrix(c(2, 4, 1, 3), nrow = 2) por sí misma utilizando la multiplicación matricial. ¿Cuál es el resultado?

  2. Dada la matriz: A <- matrix(c(3, 4, 5, 6), nrow = 2) utiliza el paquete expm para calcular la potencia 3 de la matriz. ¿Cuál es el resultado?

  3. Dada la matriz: A <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9), nrow = 3) calcula su determinante.

  4. Dada la matriz B <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2) calcula su matriz traspuesta. ¿Qué comando utilizarías y cómo sería el resultado?

  5. Dada la matriz C <- matrix(c(2, 0, 1, 3), nrow = 2) calcula su determinante. ¿Qué valor obtienes?

Sistemas de ecuaciones lineales.

R dispone de la función solve() para resolver sistemas de ecuaciones lineales. En particular, si el sistema a resolver es de la forma \(Ax=b\) (donde b podría ser también una matriz), su solución es \(x=A^{−1}b\), que en R se obtiene mediante solve(A,b).

# Definir la matriz A
A_3 <- matrix(
  c(1, 2, 3,
    1, 4, 6,
    2, -3, -5),
  nrow = 3,
  ncol = 3,
  byrow = FALSE
)

# Definir el vector b
b <- matrix(c(9, 0, 1), nrow = 3, byrow = FALSE)

# Resolver el sistema A_3 * x = b
x <- solve(A_3, b)

# Mostrar la solución
x
##      [,1]
## [1,]   29
## [2,]  -16
## [3,]   -2

Matriz inversa.

En particular, si b es la matriz identidad, la función anterior devuelve la matriz inversa de A:

A = matrix(ncol=3,c(2,4,3,5,7,1,2,2,3),byrow=TRUE)
A
##      [,1] [,2] [,3]
## [1,]    2    4    3
## [2,]    5    7    1
## [3,]    2    2    3
I=diag(1,nrow=3) # Matriz diagonal de dimension 3
B=solve(A,I) # Matriz inversa de A
B
##            [,1]          [,2]       [,3]
## [1,] -0.7307692  2.307692e-01  0.6538462
## [2,]  0.5000000  1.138690e-17 -0.5000000
## [3,]  0.1538462 -1.538462e-01  0.2307692
A%*%B
##              [,1]          [,2]          [,3]
## [1,] 1.000000e+00  0.000000e+00 -1.110223e-16
## [2,] 1.387779e-16  1.000000e+00 -1.110223e-16
## [3,] 1.942890e-16 -5.551115e-17  1.000000e+00

Podemos obviar la referencia a la matriz identidad I. Si utilizamos solamente solve(A) obtenemos también la inversa de A:

invA =solve(A)
A %*% invA
##              [,1]          [,2]          [,3]
## [1,] 1.000000e+00  0.000000e+00 -1.110223e-16
## [2,] 1.387779e-16  1.000000e+00 -1.110223e-16
## [3,] 1.942890e-16 -5.551115e-17  1.000000e+00
invA %*% A
##               [,1] [,2]         [,3]
## [1,]  1.000000e+00    0 3.330669e-16
## [2,]  0.000000e+00    1 0.000000e+00
## [3,] -5.551115e-17    0 1.000000e+00

Paquete pracma

Reducción de sistemas haciendo uso del paquete pracma.

#install.packages("pracma")
library(pracma)
## 
## Attaching package: 'pracma'
## The following objects are masked from 'package:expm':
## 
##     expm, logm, sqrtm
## The following objects are masked from 'package:Matrix':
## 
##     expm, lu, tril, triu
A <- matrix( c(1,2,3,1,4,6,2,-3,-5),byrow=F,nrow=3,ncol=3)
b <- matrix(c(9,0,1),nrow=3,ncol=1)

rref(cbind(A, b))
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0   29
## [2,]    0    1    0  -16
## [3,]    0    0    1   -2

Sistemas no cuadrados

A <- matrix(c(1,1,1,2,-1,-3), ncol=3, byrow=T)
b <- matrix(c(1,5),ncol=1,byrow=F)
rref(cbind(A,b))
##      [,1] [,2]       [,3] [,4]
## [1,]    1    0 -0.6666667    2
## [2,]    0    1  1.6666667   -1

Preguntas

  1. Resuelve el sistema de ecuaciones lineales dado por la matriz A = matrix(c(1, 2, 3, 1, 4, 6, 2, -3, -5), nrow = 3) y el vector b = matrix(c(9, 0, 1), nrow = 3) usando la función solve(). ¿Cuál es la solución?

  2. Calcula la inversa de la matriz A = matrix(c(2, 4, 3, 5, 7, 1, 2, 2, 3), nrow = 3). ¿Qué pasa cuando multiplicas A por su inversa?

  3. Utiliza la función rref() del paquete pracma para resolver el sistema no cuadrado dado por A = matrix(c(1, 1, 1, 2, -1, -3), nrow = 2, byrow = TRUE) y b = matrix(c(1, 5), nrow = 2). ¿Cuál es la solución?

  4. Resolver el siguiente sistema de ecuaciones lineales:

\[x+y+z=2\] \[2x+2y+2z=4\] \[3x+3y+3z=6\]