Disclaimer:
Antes de nada, quiero aclarar que NO soy matemático y sé 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 detectáis alguna metida de pata, ponedlo en los comentarios.

“To give the history of linear algebra is a task that is as important as it is difficult.”

Nicolas Bourbaki, 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 allá por el 1812.

Dando origen al 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 y de Cauchy.

¿Por qué es tan importante?

Leyendo el artículo de Pete Warden, “Why GEMM is at the heart of deep learning“, podemos encontrar un dato muy esclarecedor.

Entre el 89% y el 95% del tiempo de cálculo 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. 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.

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.

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 número de filas de A por el número de columnas de B. \(A^{n,k} * B^{k,m} = C^{n,m}\)
• NO es conmutativa. \(AB \neq BA\)

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 multiplicar cada uno de los elementos de esos vectores y acumular ese resultado.

#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:

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 😉

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 qué 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…

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.

¿Por qué 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 que consigue un coste asintótico de \(O(n^{2.807355})\) e incluso hay otros como el de Coppersmith–Winograd que consiguen un coste de \(O(n^{2.375477})\). Aunque en la práctica casi nunca se usa este último.

Y, ¿cómo 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 u OpenBLAS, de Intel con MKL, de AMD con BLIS o en GPU con cublas. Todas, o casi todas, las aplicaciones de cálculo se basan en esta librería.

Y además, casi todas vienen 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 en los comentarios, 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 suele 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í, sé que es raro, pero nada más esclarecedor que un poco de código:

import numpy as np

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)

Cuyo resultado será:

array([[[ 10,  13],
        [ 28,  40]],

       [[172, 193],
        [244, 274]]])

O con TensorFlow

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)

El resultado de lo anterior será:

tf.Tensor(
[[[ 1  2  3]
  [ 4  5  6]]

 [[ 7  8  9]
  [10 11 12]]], shape=(2, 2, 3), dtype=int32)
tf.Tensor(
[[[13 14]
  [15 16]
  [17 18]]

 [[19 20]
  [21 22]
  [23 24]]], shape=(2, 3, 2), dtype=int32)
tf.Tensor(
[[[ 94 100]
  [229 244]]

 [[508 532]
  [697 730]]], shape=(2, 2, 2), dtype=int32)

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.

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)

Que dará el siguiente resultado:

array([[[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.

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:

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

El resultado de lo anterior saldrá:

np. size: 2 cost: 3.647299945441773e-05 secs
np. size: 5 cost: 1.3099000170768704e-05 secs
np. size: 10 cost: 3.490800008876249e-05 secs
np. size: 50 cost: 0.00016806599978735903 secs
np. size: 100 cost: 0.0010476679999555927 secs
np. size: 300 cost: 0.025637499000367825 secs
np. size: 500 cost: 0.1418691549997675 secs
np. size: 750 cost: 0.43895936499939125 secs
np. size: 1000 cost: 1.0441638820002481 secs
np. size: 1500 cost: 3.556194990000222 secs
np. size: 2000 cost: 8.025744180000402 secs
np. size: 3000 cost: 26.964918182000474 secs
np. size: 5000 cost: 133.58672909400048 secs

Aquí podéis ver los resultados, y, junto con lo anterior,… sacad vuestras conclusiones.

Info: recordad, si queréis que hable en profundidad sobre la einsum ponedlo en los comentarios.

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?

# ATENCIÓN, DEBES MODIFICAR ESTE CÓDIGO
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, sobre todo, 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 o mejor, puedes usar los comentarios.