# Álgebra lineal para Deep Learning (Parte 2).
!!! Warn
Antes de nada, quiero aclarar que NO soy matemático y se de matemáticas lo justo para pasar el día.
Hay muchos conceptos que se me escapan y otros que trataré de manera muy superficial.
Matemáticos perdonadme, por favor. Y si detectais alguna metida de pata, avisarme.
> "To give the history of linear algebra is a task that is as important as it is difficult."
>
> -- [Nicolas Bourbaki]("https://www.curistoria.com/2010/10/los-tratados-de-matematicas-de-nicolas_14.html" target="_blank"), 1984 --
# LA OPERACIÓN
Históricamente ha existido una forma de simplificar los cálculos en el álgebra lineal, y sigue siendo indispensable en la física, ingeniería y en la informática. Esta operación es la multiplicación de matrices.
El descubrimiento de esta operación se atribuye a Augustin-Louis Cauchy y a Jacques Philippe Marie Binet alla por el 1812.
![
Augustin Louis Cauchy. [Wikimedia]("https://commons.wikimedia.org/w/index.php?curid=7059486" target="_blank")
]("images/Jacques_Binet.jpg")
Dando origen a teorema de Binet-Cauchy donde se define por primera vez la multiplicación de matrices.
Como curiosidad, podemos echar un ojo a los artículos de [Binet]("http://people.math.harvard.edu/~knill/graphgeometry/binet/binet1.pdf" target="_blank") y de [Cauchy]("http://people.math.harvard.edu/~knill/graphgeometry/binet/cauchy1.pdf" target="_blank") ;).
## ¿Por qué es tan importante?
Leyendo el artículo de Pete Warden [Why GEMM is at the heart of deep learning](https://petewarden.com/2015/04/20/why-gemm-is-at-the-heart-of-deep-learning/ target="_blank") podemos encontrar un dato muy esclarecedor.
Entre el 89% y el 95% del tiempo de calculo entrenando un modelo de Deep Learning se usa en la operación de Multiplicación de matrices.
Habréis visto que hemos usado el termino "GEMM", que significa GEneral Matrix Multiplication que son una serie de funciones preparadas para trabajar con matrices de la librería BLAS, que es el estándar para el cálculo numérico en cualquier aplicación informática.
La primera aparición de BLAS (Basic Linear Algebra Subprograms for *fortran* usage) fue allí [por el 1979](https://www.cs.utexas.edu/users/kincaid/blas.pdf target="_blank").
Bueno, pues después de esta pequeña introducción histórica, toca ir al lío.
Así que prepárate un buen café (o cerveza) y empecemos.
# Multiplicación de tensores de rango 2.
La multiplicación de tensores de rango 2 (también llamados matrices ;) ) está definida como:
$$ C_{i,j}=\sum_{k} A_{i,k}*B_{k,j} $$
o lo que es lo mismo, cada elemento de la matriz C es la multiplicación del vector linea de $A_{i}$ por el vector columna $B_{j}$.
Si todavía estás confuso, creo que esta imagen es bastante útil para entenderlo.
![Multiplicación de matrices Source [Wikipedia]("https://es.wikipedia.org/wiki/Multiplicaci%C3%B3n_de_matrices#/media/Archivo:Matrix_multiplication_diagram.svg" target="_blank")]("images/mm.png")
## Propiedades de la multiplicación de matrices.
Hay algunas características de la multiplicación de matrices que debemos grabarnos a fuego en la memoria:
- Las columnas de A deben de ser iguales a las lineas de B (En formato row-column)
- El tamaño de C es el numero de filas de A x el numero de columnas de B.
$$ A^{n,k} * B^{k,m} = C^{n,m} $$
- NO es conmutativa.
$$ A*B \neq B*A $$
## Implementando.
La versión más directa de la multiplicación de matrices es recorrer cada una de las filas de $A$ y cada una de las columnas de $B$ y por cada uno de los elementos de esos vectores multiplicarlos y acumular ese resultado.
````````` cpp
#pseudo-código
c <- inicializamos a 0 con tamaño {i,j}
for i=0 hasta A.lineas:
for j=0 hasta B.columnas
for k= 0 hasta A.columnas:
c[i,j]+=A[i,k]*B[k,j]
`````````
Ahora te toca implementarlo a ti:
````````` python
class MatMulNaive(Operation):
def __init__(self, A: Tensor,B:Tensor):
self.A = A
self.B = B
def forward(self):
cshape=(self.A.shape[0],self.B.shape[1])
C=Tensor([0 for x in range(np.prod(cshape))]).reshape(cshape)
for i in range(0, self.A.shape[0]):
# Tu código aquí
...
self.result = C
`````````
### Test unitarios
¿Ya lo tienes implementado?, ahora toca probarlo ;)
````````` python
from Tensor import Tensor
import Operation
A = Tensor([1,2,3,4,5,6]).reshape((2,3))
B = Tensor([1,2,3,4,5,6]).reshape((3,2))
C = Operation.MatMulNaive(A,B)
C.forward()
CR = Tensor ([22, 28,49, 64]).reshape((2,2))
assert(C.result == CR)
`````````
## Velocidades, coste y más.
Bueno, pues no ha sido tan difícil de implementar, ¿no?.
Entonces, ¿Por que se pone tan pesado todo el mundo con intentar hacer la multiplicación de matrices más rápida?
Puede que lo mejor sea que lo comprobéis en vuestras carnes...
````````` python
from Tensor import Tensor
import Operation
import numpy as np
from time import perf_counter
sizes = [2,5,10,50,100,300,500,750,1000,1500,2000,3000,5000]
for i in sizes:
A = Tensor([x for x in range(i*i) ]).reshape(i,i)
B = Tensor([x for x in range(i*i) ]).reshape(i,i)
C = Operation.MatMulNaive(A,B)
t0 = perf_counter()
C.forward()
t1 = perf_counter()
print("naive. size: {} cost: {} secs".format(i,t1-t0) )
for i in sizes:
A = np.array([x for x in range(i*i) ]).reshape(i,i)
B = np.array([x for x in range(i*i) ]).reshape(i,i)
t0 = perf_counter()
C = np.matmul(A,B)
t1 = perf_counter()
print("np. size: {} cost: {} secs".format(i,t1-t0) )
`````````
Para evitarnos tener que pasar un buen rato delante del ordenador esperando a que acabe esta prueba, aquí os dejo la gráfica de los costes temporales.
![Coste Multiplicación de matrices. La fuente xkcd no soporta acentos :( ]("images/costemm.png")
### ¿Por que hay tanta diferencia?
Bueno, pues nuestra implementación de la multiplicación de matrices tiene un coste algorítmico de $O(n^{3})$ y eso, como habéis podido comprobar es mucho.
Hay algoritmos mejores, como el de [Strassen](https://en.wikipedia.org/wiki/Strassen_algorithm target="_blank") que consigue un coste asintótico de $O(n^{2.807355})$ e incluso hay otros como el de [Coppersmith–Winograd](https://en.wikipedia.org/wiki/Coppersmith%E2%80%93Winograd_algorithm target="_blank") que consiguen un coste de $O(n^{2.375477})$. Aunque en la practica casi nunca se usa este último.
### ¿Y como podemos hacer que sea más rápido?
Aquí es donde entra BLAS.
Esta librería, o conjunto de operaciones, tiene unas versiones hiper-optimizadas para cada arquitectura.
Ya sea con BLAS genérica de [Netlib]("http://www.netlib.org/blas/" target="_blank") o [OpenBLAS]("https://www.openblas.net/" target="_blank"), de Intel con [MKL]("https://software.intel.com/content/www/us/en/develop/tools/math-kernel-library.html" target="_blank"), de AMD con [BLIS]("https://github.com/amd/blis" target="_blank") o en GPU con [cublas]("https://docs.nvidia.com/cuda/cublas/index.html" target="_blank") todas, o casi todas, las aplicaciones de cálculo se basan en esta librería.
Y ademas, casi todas tienen con capacidades para paralelizar los cálculos.
Por ejemplo, Numpy utiliza la versión que tengáis instalada en vuestro sistema operativo.
A no ser que alguien lo pida explicitamente, no voy a entrar en la implementación de la multiplicación de matrices con BLAS. Aunque, si tenéis curiosidad, en la versión C++ del repositorio la multiplicación de matrices está implementada con BLAS.
## ¿Y que pasa con los tensores de rango > 2?
Para explicar las distintas formas de multiplicar tensores de rango mayor de dos, deberíamos escribir un libro.
Pero estamos de enhorabuena, en el Deep Learning no se suelen usar la multiplicación tensorial como tal, si no, que se usa una forma de multiplicación llamada Multiplicación de modo N.
Que, simplificando es una multiplicación de matrices en grupo.
Si tenemos un tensor de tamaño $T^{2,3,2}$ en el modo 1, podríamos verlo como dos matrices de $A^{3,2}$ agrupadas.
Sí, se que es raro, pero nada más esclarecedor que un poco de código:
``` python
A = np.arange(12).reshape(2,2,3)
B = np.arange(12).reshape(2,3,2)
#aquí podéis ver que ha hecho 2 multiplicaciones de las matrices de 2x3 y 3x2
np.matmul(A,B)
[
[
[ 10, 13],
[ 28, 40]
],
[
[172, 193],
[244, 274]
]
]
```
O con TensorFlow
``` python
import tensorflow as tf
a = tf.constant(np.arange(1, 13, dtype=np.int32), shape=[2, 2, 3])
print(a)
# 3-D tensor
b = tf.constant(np.arange(13, 25, dtype=np.int32), shape=[2, 3, 2])
print(b)
# 3-D tensor
c = tf.matmul(a, b)
# a * b
print(c)
```
Lo bueno de esta multiplicación n-mode es que podemos multiplicar tensores de rango mayor de dos con matrices, simplemente haciendo broadcasting de la matriz.
``` python
import numpy as np
A=np.arange(12).reshape(2,2,3)
B=np.arange(6).reshape(3,2)
#Para multiplicar, B duplica para ser un tensor de 2x3x2
np.matmul(A,B)
[
[
[10, 13],
[28, 40]],
[
[46, 67],
[64, 94]
]
]
```
## Notación de Einstein (einsum)
Realmente las multiplicaciones de tensores se suelen programar con una "forma" de definir las funciones que se llama [Notación de Einstein]("https://es.wikipedia.org/wiki/Convenio_de_suma_de_Einstein" target="_blank").
No vamos a explicar eso aquí, pero está bien saber que existe y que funciona muy bien.
Como ejemplo vamos a definir la multiplicación de matrices:
``` python
import numpy as np
from time import perf_counter
sizes = [2,5,10,50,100,300,500,750,1000,1500,2000,3000,5000]
def einsmm(a,b): return np.einsum('ik,kj->ij', a, b)
for i in sizes:
A = np.array([x for x in range(i*i) ]).reshape(i,i)
B = np.array([x for x in range(i*i) ]).reshape(i,i)
t0 = perf_counter()
C = einsmm(A,B)
t1 = perf_counter()
print("np. size: {} cost: {} secs".format(i,t1-t0) )
```
Aquí podéis ver los resultados, y bueno... sacad vuestras conclusiones.
![Numpy matmul vs einsum La fuente xkcd no soporta acentos :( ]("images/einsum.png")
!!! Info:
recordad, si queréis que hable en profundidad sobre la einsum decirmelo por twitter.
# Ejercicio. Truco del sumatorio
Para terminar hoy vamos a calentarnos el coco un rato.
Hay un pequeño truco para realizar el sumatorio por filas o columnas usando multiplicación de matrices.
¿Sois capaces de sacarlo?
``` python
import numpy as np
#sumatorio por columnas
A = np.arange(12).reshape(3,4)
B = ?
C = Operation.matmul(A,B)
assert(C == np.sum(A,axis=1))
#sumatorio por filas
A = np.arange(12).reshape(3,4)
B = ?
C = Operation.matmul(B,A)
C == np.sum(A,axis=1)
```
# Unas palabras finales.
Hemos estado bastante tiempo con el álgebra lineal, y todavía nos quedan bastantes funciones por implementar. Pero para ir avanzando más rápido,dejaremos la implementación de las operaciones que necesitemos para cuando haga falta.
En el próximo post veremos que son, y que tipos hay de Grafos Computacionales y sobretodo... nos implementaremos nuestro propio Grafo Computacional.
# Dudas y preguntas.
Si tienes cualquier duda, pregunta o he metido la pata en algún sitio, puedes contactar conmigo en twitter.
-------------------------
