Chapter 8 ÁLGEBRA MATRICIAL.

En este capítulo se desarrollará operaciones matemáticas con matrices. Y en la parte final se muestra como resolver un sistema de ecuaciones lineales.

8.1 Operaciones básicas.

En la siguiente tabla se muestran las operaciones y su sintaxis.

OPERACIÓN	SINTAXIS
Adición	+
Sustracción	-
Multiplicación por un escalar	*
Producto de Matrices	%*%
Potencia de una matriz	mtx.exp()
Producto exterior	%o%
Producto Kronecker	%x%

Primero, crearemos un par de matrices que nos permitirán desarrollar estas operaciones.

# Primera matriz.
A<-matrix(1:4, ncol = 2)
A

##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4

# Segunda matriz.
B<-matrix(c(4,2,0,1), ncol = 2)
B

##      [,1] [,2]
## [1,]    4    0
## [2,]    2    1

Si se desea sumar matrices entonces se usa el operador +.

# Sumando las 2 matrices. 
A+B

##      [,1] [,2]
## [1,]    5    3
## [2,]    4    5

si usted desea restar las matrices, tendrá que usar el operador -.

# Restando las 2 matrices. 
A-B

##      [,1] [,2]
## [1,]   -3    3
## [2,]    0    3

Si se desea multiplicar por un escalar a una matriz, entonces se usará el operador *. Por ejemplo, vamos a multiplicar 5 por la matriz A.

# Multiplicación de un escalar por una matriz. 
5*A

##      [,1] [,2]
## [1,]    5   15
## [2,]   10   20

Para el caso de multiplicación de matrices, se tendrá que usar el operador %*%.

# Multiplicación de matrices. 
A%*%B

##      [,1] [,2]
## [1,]   10    3
## [2,]   16    4

También podríamos multiplicar la matriz A tres veces.

# Multiplicación de A por A por A.
A%*%A%*%A

##      [,1] [,2]
## [1,]   37   81
## [2,]   54  118

La lógica indica que si deseamos multiplicar la matriz A 4 veces, tendremos que usar A%*%A%*%A%*%A. Y si deseamos multiplicar la matriz A 5 veces A%*%A%*%A%*%A%*%A. Pero esto puede resultar tedioso si desea usted desea elevar la matriz A a la 10 o a la 15 o a un exponente mayor. Entonces, se tiene que crear una función que permita desarrollar este cálculo.

Planteo la siguiente función para la solución de este tipo de problemas.

matriz_n<-function(x, n){
  y<-1
  t<-x
  while (y+1<=n) {
    t<-t%*%x
    y<-y+1
  }
  t
}

Para la construcción de esta función estoy usando funciones como while() que se verá en el siguiente capítulo, asimismo, en los capítulos finales se mostrará como crear funciones para solucionar cualquier tipo de problema.

Pongamos en práctica la función que se acaba de crear. Vamos a elevar la matriz A al cubo.

# La matriz A al cubo.
matriz_n(A, 3)

##      [,1] [,2]
## [1,]   37   81
## [2,]   54  118

El resultado es igual al que se haría manualmente.

A%*%A%*%A

##      [,1] [,2]
## [1,]   37   81
## [2,]   54  118

Esto es una muestra de cómo podríamos solucionar este problema creando nuestras propias funciones, lo que es el reflejo fiel de la palabra programar.

Pero el problema es sencillo y otros usuarios de R, ya han programado una función similar a la que hemos desarrollado y lo tienen documentado en un paquete.

Un paquete en R, es una compilación de funciones que han sido creadas por los usuarios y se pueden descargar desde el CRAM de R. Para instalar un paquete tendremos que usar la función install.package().

Como ejemplo, vamos a instalar el paquete Biodem que tiene una función que calcula la potencia de las matrices.

# Instalar el paquete Biodem.
install.packages("Biodem")

Los paquetes sólo se instalan una vez y nada más. El siguiente paso es cargar el paquete en nuestro entorno de R. Para esto se usa la función library().

library(Biodem)

## Warning: package 'Biodem' was built under R version 4.1.1

Se tendrá que cargar los paquetes cada vez que inicies sesión en R, es decir, cada vez que usted abra el software.

Una vez cargado el paquete llamaremos a la función mtx.exp() que nos permite calcular las potencias de las matrices. Vamos a calcular lo mismo que hicimos con la función que hemos creado, es decir, elevar la matriz A al cubo.

# Elevar la matriz al cubo. 
mtx.exp(A,3)

##      [,1] [,2]
## [1,]   37   81
## [2,]   54  118

Entonces, hay dos posibilidades. La primera es crear nuestra propia función o descargar un paquete en donde haya una función que permita realizar el cálculo que deseamos

Un punto interesante es el tiempo que realiza nuestro procesador para realizar el cálculo, pero esto lo veremos en los anexos.

Por otro lado, si nosotros deseamos calcular el producto exterior de las matrices.

# Producto exterior de la matriz A.
A%o%A

## , , 1, 1
## 
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4
## 
## , , 2, 1
## 
##      [,1] [,2]
## [1,]    2    6
## [2,]    4    8
## 
## , , 1, 2
## 
##      [,1] [,2]
## [1,]    3    9
## [2,]    6   12
## 
## , , 2, 2
## 
##      [,1] [,2]
## [1,]    4   12
## [2,]    8   16

El producto exterior de matrices nos servirá para el cálculo del producto Kronecker. En ese sentido, para hallar el producto Kronecker se tendrá que usar el operador %x%.

# El producto Kronecker de la matriz A con ella misma. 
A%x%A

##      [,1] [,2] [,3] [,4]
## [1,]    1    3    3    9
## [2,]    2    4    6   12
## [3,]    2    6    4   12
## [4,]    4    8    8   16

8.2 Principales operaciones con matrices.

En la siguiente tabla se muestran las principales operaciones de las matrices y su sintaxis en el software R.

OPERACIÓN	SINTAXIS
Transpuesta	t()
Diagonal	diag()
Traza	sum(diag())
Determinante	det()
Inversa	solve()
Descomposición	qr()
Rango	qr()$rank
Descomposición de cholesky	chol()
Varianza	var()

Si nosotros queremos calcular la transpuesta de la matriz A, entonces, usaremos la función t().

# La matriz A. 
A

##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4

# Su transpuesta.
t(A)

##      [,1] [,2]
## [1,]    1    2
## [2,]    3    4

Asimismo, si se quere calcular la diagonal de la matriz A, tendremos que usar la función diag().

# La diagonal de la matriz A.
diag(A)

## [1] 1 4

Como sabemos la traza de una matriz es la suma de la diagonal. Entonces, si queremos calcular la traza de la matriz, sólo tendríamos que sumar los valores de los elementos de la diagonal.

# Calculando la traza de la matriz. 
sum(diag(A))

## [1] 5

Pero si usted desea crear la función para la traza, lo podría hacer de la siguiente manera.

# Creando la función que calcule la traza. 
traza<-function(x){
  sum(diag(x))
}

# Aplicando la función que calcula la traza.
traza(A)

## [1] 5

Lo que nos resulta en el mismo resultado.

Uno de los calculos muy importantes en matrices cuadradas es el cálculo de la determinante. En R para calcular la determinante con la función det().

# Calculando la determinante de la matriz A. 
det(A)

## [1] -2

Asimismo, si deseamos calcular la inversa de la matriz, entonces, usaremos la función solve().

# Calculando la inversa de la matriz A.
solve(A)

##      [,1] [,2]
## [1,]   -2  1.5
## [2,]    1 -0.5

En álgebra matricial también es muy importante descomponer las matrices (factorizar), en R se puede descomponer matrices con el método de descomposicón QR el cual es el producto de la matriz Q (matriz ortogonal) y la matriz R (matriz triángulo superior). Para lo cual se usa la función qr().

# Calculando la descomposición de una matriz. 
qr(A)

## $qr
##            [,1]       [,2]
## [1,] -2.2360680 -4.9193496
## [2,]  0.8944272 -0.8944272
## 
## $rank
## [1] 2
## 
## $qraux
## [1] 1.4472136 0.8944272
## 
## $pivot
## [1] 1 2
## 
## attr(,"class")
## [1] "qr"

Primero si nos damos cuenta el resultado es una lista que contiene 4 objetos.

class(qr(A))

## [1] "qr"

Al ser una lista podemos llamar a cada uno de los objetos con el doble corchete [[]] o con el dólar $ usando el nombre de los objetos.

El objeto qr es una matriz que contiene a la matriz R, es el triángulo superior de esta matriz y a la matriz Q, que es el triángulo superior pero de forma compacta.

Si nosotros queremos ver cuales son las matrices Q y R específicamente, tendremos que usar las funciones qr.Q() y qr.R(), en donde el argumento de las dos funciones es el objeto qr.

# Calculando la matriz Q.
qr.Q(qr(A))

##            [,1]       [,2]
## [1,] -0.4472136 -0.8944272
## [2,] -0.8944272  0.4472136

# Calculando la matriz R.
qr.R(qr(A))

##           [,1]       [,2]
## [1,] -2.236068 -4.9193496
## [2,]  0.000000 -0.8944272

Perfecto, hemos obtenido la matriz Q que es una matriz ortogonal y la matriz R que es una matriz triángulo superior.

Asimismo, podemos hacer la comprobación. Sabemos que ha descompuesto nuestra matriz A en dos matrices en: Q y R. Entonces si multplicamos la matriz Q por la matriz R deberíamos de obtener la matriz A.

# Comprobando que la descomposición es la correcta. 
qr.Q(qr(A))%*%qr.R(qr(A))

##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4

En efecto, tenemos la matriz A, por lo cual el método de descomposición queda comprobado.

Ahora si deseamos saber el rango de la matriz, usaremos la misma función qr() y llamaremos al objeto rank.

# Calculando el rango de la matriz.
qr(A)$rank

## [1] 2

Que nos indica que es de rango 2.

También podemos descomponer las matrices por el método de cholesky, pero lo primordial para poder desarrollar este método es que la matriz se simétrica y definida positiva. Por lo que para mostrar un ejemplo, vamos a crear una matriz simétrica y definida positiva.

C<-matrix(c(4,1,1,4), ncol = 2)
C

##      [,1] [,2]
## [1,]    4    1
## [2,]    1    4

El método de Cholesky dice que descompondrá a la matriz en dos matrices, la primera en una matriz triángulo inferior y la segunda es la transpuesta de la matriz triángulo inferior.

\[X=R*R'\]

En el software R se usa la función chol(), pero el resultado que arroja es la matriz triángulo superior, ya que usa la siguiente fórmula para descomponer.

\[X=R'*R\]

Donde: $R$ es la matriz triángulo inferior.

Es así que el software R arroja $R'$, por eso se convierte en una matriz triángulo superior.

# Determinando la matriz triángulo superior por la descompisición
# de Cholesky.
chol(C)

##      [,1]     [,2]
## [1,]    2 0.500000
## [2,]    0 1.936492

Comprobando que el resultado de la factorización es el correcto.

# Comprobando.
t(chol(C))%*%chol(C)

##      [,1] [,2]
## [1,]    4    1
## [2,]    1    4

Y en efecto, obtenemos la matriz original, por lo que la descomposición de Cholesky es la correcta.

Por último, para el cálculo de la matriz de varianzas y covarianzas se usará la función var(). Vamos a calcular la matriz de varianas y covarianzas de la matriz A.

var(A)

##      [,1] [,2]
## [1,]  0.5  0.5
## [2,]  0.5  0.5

8.3 Solución a sistemas de ecuación.

El software R también nos permite resolver sistemas de ecuaciones. Por ejemplo, si tenemos el siguiente sistema de ecuaciones.

\[2x+3y=1\]

\[3x-7y=2\]

Para poder desarrollarlo se puede hacer mediante el cálculo matricial, lo que se tendrá que hacer es convertir las dos ecuaciones en matrices.

\[\begin{pmatrix} 2 & 3\\ 3 & 7 \end{pmatrix}\begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} 1\\ 2 \end{pmatrix}\]

Las matrices que nos interesan son la primera y la tercera. Las creamos.

# La primera matriz.
D<-matrix(c(2,3,3,-7), nrow = 2, byrow = TRUE)
D

##      [,1] [,2]
## [1,]    2    3
## [2,]    3   -7

# La segunda Matriz.
v<-matrix(1:2, ncol = 1)
v

##      [,1]
## [1,]    1
## [2,]    2

Como se puede ver las matrices son iguales a las ecuaciones que se muestran.

Para poder obtener respuesta se hace usa la función solve().

# Para solucionar el sistema de ecuaciones.
solucion<-solve(D,v)    
solucion

##             [,1]
## [1,]  0.56521739
## [2,] -0.04347826

Guardamos el resultado en el objeto solucion y podemos ver que nos arroja los valores de x y y que resuelven el sistema de ecuaciones.

Podemos cambiar los nombres de las filas con el fin de obtener un resultado más accesible a los ojos.

dimnames(solucion)<-list(c("x", "y"), NULL)
solucion

##          [,1]
## x  0.56521739
## y -0.04347826

Perfecto, ahora se muestran los resultados de x y y que resuelven las ecuaciones.

8.4 Valores y vectores propios.

Si usted desea calcular, los valores y vectores propios de una matriz tendrá que usar la función eigen().

Para el ejemplo se calculará los valores y vectores propios de la matriz C, la que se usó para la factorización con el método de cholesky.

# Para el cálculo de valores y vectores propios.
eigen(C)

## eigen() decomposition
## $values
## [1] 5 3
## 
## $vectors
##           [,1]       [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068  0.7071068

El objeto que hemos obtenido es una lista, por lo cual, si sólo deseamos los valores propios se tendría que usar la siguiente sintaxis eigen()$values. Asimismo, si sólo deseo los vectores propios, la sintaxis que se tendría que usar seríaeigen()$vectors.

# Sólo para los valores propios.
eigen(C)$values

## [1] 5 3

# Sólo para los vectores propios.
eigen(C)$vectors

##           [,1]       [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068  0.7071068

En los últimos capítulos se verá aplicaciones en economía, específicamente en la maximización de la utilidad del consumidor, minimización de los costos de producción, modelos de regresión lineal, entre otras, usando matrices. Por el momento no podemos hacerlo, ya que aún no hemos desarrollado el capítulo de calculo diferencial.