Notas para contenedor de docker:
Comando de docker para ejecución de la nota de forma local:
nota: cambiar <ruta a mi directorio> por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker.
docker run --rm -v <ruta a mi directorio>:/datos --cap-add SYS_ADMIN --privileged --name jupyterlab-numerical -p 8888:8888 -d palmoreck/jupyterlab_numerical:1.1.0password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_numericalDocumentación de la imagen de docker palmoreck/jupyterlab_numerical:1.1.0 en liga.
Para ejecutar el comando de linux perf, que se revisará en esta nota, obsérvese que es necesario correr el contenedor de docker descrito antes con las flags --cap-add SYS_ADMIN y --privileged. Además, sólo es posible obtener las métricas descritas en esta nota vía perf con un host que tenga un sistema operativo ubuntu
El cómputo matricial está construído sobre una jerarquía de operaciones del álgebra lineal:
Productos punto involucran operaciones escalares de suma y multiplicación (nivel BLAS 1).
La multiplicación matriz-vector está hecha de productos punto (nivel BLAS 2).
La multiplicación matriz-matriz utiliza colecciones de productos matriz-vector (nivel BLAS 3).
Las operaciones anteriores se describen en el álgebra lineal con la teoría de espacios vectoriales pero también es posible describirlas en una forma algorítmica. Ambas descripciones se complementan una a la otra.
Manejaremos nombres que en el Linear Algebra Package: LAPACK son utilizados para denotar operaciones con escalares, vectores o matrices. Ver Reference-LAPACK / lapack para su github.
Consideramos $x,y \in \mathbb{R}^n$. El producto punto entre $x$ y $y$ es $c = x^Ty = \displaystyle \sum_{i=1}^n x_iy_i$.
Ejemplo y algoritmo del producto punto:
In [1]:
c=0
n=5
x=[-1]*n
y=[1.5]*n
for i in range(n):
c += x[i]*y[i]
In [2]:
c
Out[2]:
Obs:
*Los flops que realiza un algoritmo es una forma de cuantificar el volumen de trabajo asociado con un cálculo. Un flop es una operación de punto flotante: suma, multiplicación o división. Por ejemplo, en la línea:
C[i][j] = C[i][j] + A[i][k]*B[k][j]se realizan $2$ flops. Los flops sólo representan una componente para categorizar a los algoritmos de acuerdo al trabajo que realizan, otras componentes son la transferencia o movimientos de datos, data movement/motion, ejecución secuencial o en paralelo y el data locality y data reuse que realizan.
Consideramos $\alpha \in \mathbb{R}, x,y \in \mathbb{R}^n$. El nombre lo recibe por scalar alpha x plus y. En LAPACK saxpy se escribe en forma update:
$$y=\alpha x + y \therefore y_i = \alpha x_i + y_i \forall i=1,...,n$$Obs:
El símbolo $=$ no se utiliza como igualdad de expresiones sino para denotar asignación (como en computación al escribir un algoritmo).
También encontramos en LAPACK caxpy o daxpy para el caso complejo y para números en doble precisión respectivamente.
Ésta operación realiza un trabajo de $\mathcal{O}(n)$ sobre una cantidad de datos $\mathcal{O}(n)$.
Ejemplo y algoritmo de saxpy:
In [3]:
alpha=2
n=5
x=[-2]*n
y=[0]*n
for i in range(n):
y[i] += alpha*x[i]
In [4]:
y
Out[4]:
o en una forma update:
In [5]:
alpha=2
n=5
x=[-2]*n
y=[3,4,-1,0,1]
for i in range(n):
y[i] += alpha*x[i]
In [6]:
y
Out[6]:
Comentario: La operación de producto punto y saxpy son algoritmos catalogados como de nivel BLAS 1, ver BLAS: Basic Linear Algebra Subprograms. Éstos algoritmos se caracterizan por involucrar una cantidad de trabajo lineal sobre una cantidad lineal de datos. Ver level 1 para más ejemplos de este tipo de algoritmos.
Consideramos $A \in \mathbb{R}^{m \times n}, x \in \mathbb{R}^n, y \in \mathbb{R}^m$. La operación $y = y + Ax$ es una operación generalizada saxpy, por ello se denomina gaxpy pero en LAPACK podemos encontrarla con nombres como sgemv, dgemv, cgemv o zgemv para los casos de precisión simple, doble o números complejos respectivamente. Hay diferentes formas de visualizar y escribir el algoritmo de multiplicación matriz-vector. Por ejemplo para una matriz $A$ con entradas:
In [7]:
m=2
n=5
A=[[1.2]*n if i%2==0 else [1]*n for i in range(m)]
In [8]:
A
Out[8]:
In [9]:
A[0][0]
Out[9]:
In [10]:
A[1][n-1]
Out[10]:
se tiene:
Ejemplo
In [11]:
x=[2]*n
y=[0]*m
for i in range(m):
for j in range(n):
y[i]+=A[i][j]*x[j]
In [12]:
y
Out[12]:
Si $y$ tiene valores distintos de $0$, se realiza un update:
In [13]:
x=[2]*n
y=[-1]*m
for i in range(m):
for j in range(n):
y[i]+=A[i][j]*x[j]
In [14]:
y
Out[14]:
Comentarios:
for i in range(m)
y[i]+=A[i,:]*x #producto puntodonde: $A[i,:]$ es el $i$-ésimo renglón de $A$. Así podemos reescribir de forma más compacta este algoritmo.
También puede escribirse al algoritmo gaxpy en una forma orientada por columnas:
Este algoritmo ayuda a visualizar al producto matriz-vector como una combinación lineal de las columnas de $A$:
$$Ax = \displaystyle \sum_{j=1}^n A_jx_j$$con $A_j$ la $j$-ésima columna de $A$.
Ejemplo
In [15]:
x=[2]*n
y=[0]*m
for j in range(n):
for i in range(m):
y[i]+=A[i][j]*x[j]
In [16]:
y
Out[16]:
Comentarios:
El algoritmo de multiplicación matriz-vector (versión row o column oriented) involucra $\mathcal{O}(mn)$ operaciones o una cantidad cuadrática de trabajo, que podemos entender como "si duplicamos cada dimensión de $A$ entonces la cantidad de trabajo se incrementa por un factor de $4$". Tal número de operaciones trabajan sobre una matriz o sobre una cantidad cuadrática de datos. A los algoritmos que realizan una cantidad cuadrática de trabajo sobre una cantidad cuadrática de datos se les cataloga de nivel BLAS 2. Ver level 2 para más ejemplos de algoritmos en el álgebra lineal en esta categoría.
En el algoritmo gaxpy column oriented el acceso a la matriz $A$ es por columna.
La versión column oriented se puede analizar desde el punto de vista puramente algorítmico como un intercambio entre las líneas con los índices $i$ y $j$ de cada loop y un acceso a los datos de la matriz por columna. O bien, se puede analizar desde el álgebra lineal indicando que el vector $y$ está en el espacio generado por las columnas de $A$ y cuyas coordenadas son dadas por las entradas del vector $x$:
Una ejemplo de visualización del espacio generado por las columnas de $A \in \mathbb{R}^{3 \times 2}$, llamado rango o imagen de $A$, $Im(A)$, es el siguiente:
En este dibujo los vectores $b, r(x) \in \mathbb{R}^3$ no están en $Im(A) \subset \mathbb{R}^3$ pero $Ax$ en azul sí. $Im(A)$ en el dibujo es un plano en el espacio $\mathbb{R}^3$.
for j in range(n)
y+=A[:,j] * x[j]sin embargo como hemos visto en Python con su implementación más común CPython, no es posible realizar tal indexado:
In [17]:
x=[2]*n
y=[0]*m
for j in range(n):
y+=A[:,j]*x[j]
a menos que incorporemos alguna paquetería que permita la vectorización y el uso de índices para extracción de columnas (o renglones) de $A$. Entre las razones que encontramos en Python el no soporte para la vectorización está que el bytecode de Python no está optimizado para vectorización. Un ejemplo de un paquete que permite realizar operaciones de forma vectorizada es numpy:
In [18]:
import numpy as np
In [19]:
x = 2*np.ones(n)
y = np.zeros(m)
In [20]:
x
Out[20]:
In [21]:
y
Out[21]:
In [23]:
A=np.array([[1.2,1.2,1.2,1.2,1.2],[1,1,1,1,1]])
In [24]:
A
Out[24]:
In [25]:
for j in range(n):
y+=A[:,j]*x[j]
In [26]:
y
Out[26]:
Asimismo, el algoritmo gaxpy row oriented puede escribirse de forma más compacta haciendo uso de la definición de producto punto estándar: $x^Ty$ para dos vectores columna $x$ y $y$. En el caso de una matriz $A$ se tiene:
for i=1:m
y[i]+=A[i,:]^T*xdonde: $A[i,:]$ es el $i$-ésimo renglón de $A$. En Python:
In [76]:
x = 2*np.ones(n)
y = np.zeros(m)
A=np.array([[1.2,1.2,1.2,1.2,1.2],[1,1,1,1,1]])
In [81]:
for i in range(m):
y[i]+=A[i,:].dot(x)
In [82]:
y
Out[82]:
en donde se utilizó la función numpy.dot de numpy.
Comentarios:
numpy soporta tales operaciones.*El incremento en IPC vía la vectorización se enfoca a que las instrucciones tiendan a completar trabajo útil por ciclo pues podría darse la situación en la que se tiene: a high rate of instructions, but a low rate of actual work completed. Recuérdese que un ciclo involucra los pasos de leer una instrucción, determinar acciones a realizar por tal instrucción y ejecutar las acciones, ver liga.
En Linux existe la herramienta perf que nos ayuda a calcular métricas de desempeño de la CPU o de los cores. Para la ejecución de los siguientes ejemplos es indispensable correr el contenedor de docker descrito al inicio de la nota en un sistema ubuntu.
Ejemplo de cálculo de la norma $2$ al cuadrado de un vector
In [27]:
%%file norm_square.py
n=10**5
vector=list(range(n))
norm=0
for v in vector:
norm+=v*v
In [28]:
%%bash
sudo perf stat -S -e cycles,instructions,cache-references,cache-misses -r 20 python3 norm_square.py
Comentarios:
Con perf se repiten las mediciones utilizando la flag -r. Con la flag -e se enlistan las métricas a calcular. -S llama a sync antes de iniciar la ejecución del programa.
En el ejemplo anterior además de los ciclos e instrucciones se calculan los cache references y los cache misses. Ver 2.1.Un_poco_de_historia_y_generalidades para definiciones de cache references y misses. Esencialmente el sistema de memoria buscó el número que aparece en cache-references por datos o instrucciones en el caché y de éstas búsquedas, en cache-misses, se repora el número (y porcentaje) de búsquedas fallidas que no estuvieron en memoria.
Un factor que incrementa el número de cache-misses es la fragmentación de datos, ver Fragmentation). La fragmentación incrementa el número de transferencias de memoria hacia la CPU y además imposibilita la vectorización porque el caché no está lleno. El caché puede llenarse sólo al tener bloques contiguos de memoria alojados para los datos (recuérdese que la interconexión o bus sólo puede mover chunks de memoria contigua, ver 2.1.Un_poco_de_historia_y_generalidades). Python en su implementación común CPython, tiene data fragmentation por ejemplo al usar listas. Las listas de Python alojan locaciones donde se pueden encontrar los datos y no los datos en sí. Al utilizar tales estructuras para operaciones matriciales el runtime de Python debe realizar lookups por índices y al no encontrarse de forma contigua los datos, el overhead en la transferencia de datos es mayor lo que causa que se tengan mayor número de cache-misses y que los cores tengan que esperar hasta que los datos estén disponibles en el caché.
Obsérvese que se reporta el IPC al lado de instructions con nombre insn per cycle. Un alto IPC indica una alta transferencia de instrucciones de la unidad de memoria hacia la unidad de cómputo y un menor IPC indica más stall cycles, ver pipeline stall. Sin embargo como se mencionó antes, se debe evaluar si a high rate of instructions indica un high rate of actual work completed (p.ej. un loop tiene una alta tasa de IPC pero no siempre se realiza trabajo útil). Ver sección CPU statistics en la liga. perf también tiene una métrica para medir el número de ciclos que se utilizaron para esperar a ejecutar instrucciones (los cores están stalled). Las métricas son: stalled-cycles-frontend y stalled-cycles-backend.
También podemos obtener las estadísticas por core. Obsérvese que la salida anterior de perf calculó el IPC a partir del número de instrucciones y ciclos ahí reportados. El mismo cálculo se realiza a continuación:
In [29]:
%%bash
sudo perf stat -S -a --per-core -e cycles,instructions,cache-references,cache-misses -r 20 python3 norm_square.py
Otras métricas pueden ser obtenidas si ejecutamos perf sólo con la flag -r:
In [30]:
%%bash
sudo perf stat -S -r 20 python3 norm_square.py
Comentarios:
La métrica de task-clock indica cuántos clock cycles tomó nuestra tarea y se reporta en milisegundos. Obsérvese que también se indica cuántos CPU's fueron utilizados. En el output anterior no es exactamente 1 pues el programa no sólo involucra trabajo de CPU sino también de alojamiento de memoria.
context-switches y cpu-migrations indican cuánto se detuvo el programa para realizar:
la idea es que nuestro programa tenga un número pequeño de éstas métricas.
page-faults se relaciona con el alojamiento de memoria a partir del kernel del sistema operativo. Ocasiona que se detenga la ejecución del programa y por tanto deseamos que ésta métrica no sea muy grande. Ver Page fault.
Ver 2.1.Un_poco_de_historia_y_generalidades para el branching. Esencialmente interesa que el número de branch misses sea pequeño.
Podemos generar la salida anterior por core:
In [31]:
%%bash
sudo perf stat -S -a --per-core -r 20 python3 norm_square.py
Para contrastar las salidas anteriores con un mejor uso del caché y de los cores se utiliza a continuación el paquete de numpy:
In [33]:
%%file norm_square_numpy.py
import numpy as np
n=10**5
vector=np.arange(n)
vector.dot(vector)
In [32]:
%%bash
sudo perf stat -S -e cycles,instructions,cache-references,cache-misses -r 20 python3 norm_square_numpy.py
Comentarios:
Obsérvese que es más tardado este programa para el número $n$ que se está utilizando en comparación con el programa anterior sin vectorizar.
El número de instrucciones es menor al reportado anteriormente y parece tener un mayor número $%$ de cache-misses. Obsérvese que para este ejemplo es más ilustrativo observar las estadísticas por core y realizar conclusiones:
In [33]:
%%bash
sudo perf stat -S -a --per-core -e cycles,instructions,cache-references,cache-misses -r 20 python3 norm_square_numpy.py
Comentarios:
Se observa que la métrica IPC está más uniforme a través de los cores lo que indica un mejor load balancing.
Además, el $\%$ de cache-misses es menor que en el caso no vectorizado a diferencia del output anterior que indicaba un mayor $\%$ de cache-misses.
También podemos obtener las salidas anteriores de perf al usar numpy:
In [34]:
%%bash
sudo perf stat -S -r 20 python3 norm_square_numpy.py
In [35]:
%%bash
sudo perf stat -S -a --per-core -r 20 python3 norm_square_numpy.py
Utilizaremos perf en lo que sigue para evaluar las métricas de IPC, cache-references, cache-misses y medir el tiempo de ejecución. Además, esto permitirá comparar los tiempos de ejecución del algoritmo gaxpy (nivel BLAS 2), con vectorización y sin vectorización.
Consideramos una matriz $A \in \mathbb{R}^{10^4 \times 10^4}$ con entradas pseudoaleatorias. No se sugiere ejecutar los siguientes ejemplos para máquinas que tengan menos de 8gb de memoria:
In [22]:
np.random.seed(2020)
m=10**4
n=10**4
A=np.random.rand(m,n)
file='A.txt'
np.savetxt(file,A)
Arhivo sin vectorizar y utilizando listas para construir a la matriz y a los vectores:
In [8]:
%%file mult_matrix_vector.py
m=10**4
n=10**4
x=[2.5]*n
y=[0]*m
A = []
file='A.txt'
with open(file,'r') as f:
for l in f:
A.append([float(k) for k in l.replace('\n','').replace(' ',',').split(',')])
for j in range(n):
for i in range(m):
y[i]+=A[i][j]*x[j]
Archivo vectorizando con numpy:
In [9]:
%%file mult_matrix_vector_numpy.py
import numpy as np
m=10**4
n=10**4
x = 2.5*np.ones(n)
y = np.zeros(m)
file='A.txt'
A = np.loadtxt(file)
for j in np.arange(n):
y+=A[:,j]*x[j]
In [10]:
%%bash
sudo perf stat -S -e cycles,instructions,cache-references,cache-misses -r 2 python3 mult_matrix_vector.py
Métricas por core:
In [11]:
%%bash
sudo perf stat -S -a --per-core -e cycles,instructions,cache-references,cache-misses -r 2 python3 mult_matrix_vector.py
In [12]:
%%bash
sudo perf stat -S -e cycles,instructions,cache-references,cache-misses -r 2 python3 mult_matrix_vector_numpy.py
In [13]:
%%bash
sudo perf stat -S -a --per-core -e cycles,instructions,cache-references,cache-misses -r 2 python3 mult_matrix_vector_numpy.py
Comentarios:
No se observa una gran diferencia entre el caso vectorizado y no vectorizado salvo que en el caso vectorizando se tiene un $\%$ mayor de cache-misses. Sería conveniente realizar más repeticiones de esta misma operación o con diferentes operaciones de nivel 2 BLAS para realizar conclusiones más sólidas. Ver level 2 para más ejemplos de algoritmos en el álgebra lineal en esta categoría.
Otra opción que se tiene es utilizar la siguiente implementación ya que numpy soporta operaciones del tipo matriz-vector sin necesidad de utilizar un for loop:
In [14]:
%%file mult_matrix_vector_numpy_2.py
import numpy as np
m=10**4
n=10**4
x = 2.5*np.ones(n)
y = np.zeros(m)
file='A.txt'
A = np.loadtxt(file)
y=A@x
In [15]:
%%bash
sudo perf stat -S -e cycles,instructions,cache-references,cache-misses -r 2 python3 mult_matrix_vector_numpy_2.py
In [16]:
%%bash
sudo perf stat -S -a --per-core -e cycles,instructions,cache-references,cache-misses -r 2 python3 mult_matrix_vector_numpy_2.py
O con una versión gaxpy row oriented
In [5]:
%%file mult_matrix_vector_numpy_row_oriented.py
import numpy as np
m=10**4
n=10**4
x = 2.5*np.ones(n)
y = np.zeros(m)
file='A.txt'
A = np.loadtxt(file)
for i in np.arange(m):
y[i]+=A[i,:].dot(x)
In [6]:
%%bash
sudo perf stat -S -e cycles,instructions,cache-references,cache-misses -r 2 python3 mult_matrix_vector_numpy_row_oriented.py
In [7]:
%%bash
sudo perf stat -S -a --per-core -e cycles,instructions,cache-references,cache-misses -r 2 python3 mult_matrix_vector_numpy_row_oriented.py
Otras comparaciones que podemos realizar entre los tiempos de ejecución, cache-references y cache-misses utilizando la versión vectorizada y no vectorizada, es con el algoritmo de multiplicación de matrices $C = C + AB$ con $A \in \mathbb{R}^{m \times r}, B \in \mathbb{R}^{r \times n}, C \in \mathbb{R}^{m \times n}$. Este algoritmo se cataloga como de nivel 3 de BLAS pues realiza una cantidad de trabajo cúbica sobre una cantidad cuadrática de datos. Ver level 3 para más ejemplos de algoritmos en el álgebra lineal en esta categoría.
El algoritmo de multiplicación de matrices puede escribirse en diferentes versiones. Por ejemplo la versión i,j,k es:
for i in range(m):
for j in range(n):
for k in range(r):
C[i][j]+=A[i][k]*B[k][j]y la versión j,k,i es:
for j in range(n):
for k in range(r):
for i in range(m):
C[i][j]+=A[i][k]*B[k][j]Comentarios:
Cualquiera de las versiones (hay $3!$ de ellas) involucra una cantidad de trabajo del orden $\mathcal{O}(mnr)$, la cual es cúbica. Es posible interpretar ésta cantidad de trabajo con una frase del tipo "si duplicamos cada dimensión de $A$ entonces la cantidad de trabajo se incrementa por un factor de $8$".
Distintas versiones tienen distinto patrón de acceso a los datos de $A, B$ y $C$. Por ejemplo, para la variante i,j,k el inner loop realiza un producto punto que requiere acceso a renglones de $A$ y columnas de $B$ en una forma izquierda-derecha y abajo-arriba respectivamente. La variante j,k,i involucra una operación saxpy en el inner loop y acceso por columnas a la matriz $C$ y a $A$. Un resúmen de lo anterior se presenta en el siguiente cuadro:
Consideramos matrices $A, B \in \mathbb{R}^{10^4 \times 10^4}$ con entradas pseudoaleatorias. No se sugiere ejecutar los siguientes ejemplos para máquinas que tengan menos de 8gb de memoria. Si se ejecutó el ejemplo de gaxpy de nivel 2 de BLAS, no es necesario volver a ejecutar la siguiente celda para crear a la matriz A
In [21]:
np.random.seed(2020)
m=10**4
r=10**4
A=np.random.rand(m,r)
fileA='A.txt'
np.savetxt(fileA,A)
In [44]:
np.random.seed(2021)
r=10**4
n=10**4
B=np.random.rand(r,n)
fileB='B.txt'
np.savetxt(fileB,B)
In [45]:
%%file mult_matrix_matrix.py
m=10**4
r=10**4
n=10**4
A = []
B = []
fileA='A.txt'
fileB='B.txt'
with open(fileA,'r') as f:
for l in f:
A.append([float(k) for k in l.replace('\n','').replace(' ',',').split(',')])
with open(fileB,'r') as f:
for l in f:
B.append([float(k) for k in l.replace('\n','').replace(' ',',').split(',')])
C=[[0]*n for i in range(m)]
for i in range(m):
for j in range(n):
for k in range(r):
C[i][j]+=A[i][k]*B[k][j]
In [46]:
%%file mult_matrix_matrix_numpy.py
import numpy as np
m=10**4
r=10**4
n=10**4
fileA='A.txt'
fileB='B.txt'
A = np.loadtxt(fileA)
B = np.loadtxt(fileB)
C = A@B
In [47]:
%%bash
sudo perf stat -S -a --per-core -e cycles,instructions,cache-references,cache-misses -r 1 python3 mult_matrix_matrix.py
... 36 min y todavía no terminaba... cambiar a r 1
In [1]:
%%bash
sudo perf stat -S -a -e cycles,instructions,cache-references,cache-misses -r 1 python3 mult_matrix_matrix.py
In [3]:
%%bash
sudo perf stat -S -a -e cycles,instructions,cache-references,cache-misses -r 2 python3 mult_matrix_matrix_numpy.py
In [2]:
%%bash
sudo perf stat -S -a --per-core -e cycles,instructions,cache-references,cache-misses -r 2 python3 mult_matrix_matrix_numpy.py
Si utilizamos la versión i,j,k podemos reescribir la multiplicación de matrices en una forma:
for i in range(m):
for j in range(n)
C[i][j]+= A[i,:].dot(B[:,j])y hemos vectorizado el inner loop para usar productos punto que es una operación nivel 1 de BLAS.
Además numpy nos provee de funcionalidad para reescribir ésta forma de productos punto en una como sigue:
for i in range(m):
C[i,:]+= A[i,:]@BObs: obsérvese que en esta reescritura en el loop se realizan $m$ operaciones gaxpy de nivel 2. Esta versión como se verá a continuación es más rápida que la que utiliza operaciones de nivel 1 de BLAS:
Ejemplo
Utilizamos matrices $A, B \in \mathbb{R}^{10^3 \times 10^3}$ pseudoaleatorias:
In [19]:
np.random.seed(2020)
m=10**3
r=10**3
A=np.random.rand(m,r)
fileA_10_3='A_10_3.txt'
np.savetxt(fileA_10_3,A)
In [20]:
np.random.seed(2021)
m=10**3
r=10**3
B=np.random.rand(m,r)
fileB_10_3='B_10_3.txt'
np.savetxt(fileB_10_3,B)
In [23]:
%%file mult_matrix_matrix_numpy_dot_product.py
import numpy as np
m=10**3
n=10**3
fileA_10_3='A_10_3.txt'
fileB_10_3='B_10_3.txt'
A = np.loadtxt(fileA_10_3)
B = np.loadtxt(fileB_10_3)
C = np.zeros((m,n))
for i in np.arange(m):
for j in np.arange(n):
C[i][j]+= A[i,:].dot(B[:,j])
In [26]:
%%bash
sudo perf stat -S -a -e cycles,instructions,cache-references,cache-misses -r 5 python3 mult_matrix_matrix_numpy_dot_product.py
In [27]:
%%bash
sudo perf stat -S -a --per-core -e cycles,instructions,cache-references,cache-misses -r 5 python3 mult_matrix_matrix_numpy_dot_product.py
In [30]:
%%file mult_matrix_matrix_numpy_dot_product_gaxpy.py
import numpy as np
m=10**3
n=10**3
fileA_10_3='A_10_3.txt'
fileB_10_3='B_10_3.txt'
A = np.loadtxt(fileA_10_3)
B = np.loadtxt(fileB_10_3)
C = np.zeros((m,n))
for i in np.arange(m):
C[i,:] = A[i,:]@B
In [31]:
%%bash
sudo perf stat -S -a -e cycles,instructions,cache-references,cache-misses -r 5 python3 mult_matrix_matrix_numpy_dot_product_gaxpy.py
In [32]:
%%bash
sudo perf stat -S -a --per-core -e cycles,instructions,cache-references,cache-misses -r 5 python3 mult_matrix_matrix_numpy_dot_product_gaxpy.py
Comentarios:
Obsérvese que en las implementaciones anteriores ambas versiones utilizan la vectorización mediante numpy.
A continuación se muestra una tabla que presenta el número de flops realizados por distintas operaciones del álgebra lineal:
perf indican que la versión de gaxpy es más rápida que la versión con productos punto (nivel 1 de BLAS) para la multiplicación de matrices ya que aprovecha mejor el data reuse y data locality. Esto se logra pues gaxpy disminuye el número de cache misses y por tanto el tráfico hacia y desde el caché, el llamado data movement/motion. El data reuse y el data locality como se vio en el ejemplo pasado se incrementa al utilizar operaciones de nivel BLAS mayores. También disminuye el data motion hacia y desde el caché lo que disminuye el tráfico de datos o data movement/motion pues se atiende uno de los cuellos de botella que se revisó en 2.1.Un_poco_de_historia_y_generalidades conocido con el nombre de bottleneck de Von Neumann.
Por lo anterior hay algoritmos que se diseñan con el objetivo de aprovechar lo más posible el data reuse y el data locality y se nombran blocking algorithms. Entre los blocking algorithms encontramos a los que trabajan con matrices por bloques pues son más eficientes al utilizar un nivel $3$ de BLAS. Una matriz $A \in \mathbb{R}^{m \times n}$ por bloques se puede escribir como:
donde: $m_1 + \dots + m_q = m$, $n_1 + \dots + n_r = n$. Con esta definición se llama a la matriz $A$, una matriz por bloques de tamaño $q \times r$.
Comentarios:
Hay que tener en cuenta que trabajar con matrices por bloques utiliza mayor cantidad de memoria que un algoritmo que opera a un nivel escalar esto puede ser benéfico o no dependiendo del problema, máquina o arquitectura en la que se esté trabajando.
Trabajar con blocking algorithms permite simplificar notación matemática. Su escritura ayuda a pensar en una implementación con cómputo en paralelo.
Entre las operaciones que son posibles realizar a un nivel por bloques se encuentran:
con $\mu \in \mathbb{R}$.
Obs: el número de columnas de los bloques $A_{11}, A_{21}, A_{31}$ deben coincidir con el número de renglones de $B_{11}, B_{12}$. Así como con los bloques $A_{12}, A_{22}, A_{32}$ y $B_{21}, B_{22}$.
Considérese que se particionan a las matrices $A, B, C$ como sigue:
Esto es, $A,B,C$ son $N \times N$ matrices por bloques con $\ell \times \ell$ bloques. Entonces, para índices $\alpha = 1,2,\dots, N$ y $\beta = 1,2,\dots,N$ se tiene que el bloque $\alpha \beta$ de $C$ se obtiene: $C_{\alpha \beta} = \displaystyle \sum_{\gamma=1}^N A_{\alpha \gamma} B_{\gamma \beta}$.
El algoritmo es:
for alpha in np.arange(N)
i = np.arange((alpha - 1)*l + 1,alpha*l + 1) #hay que revisar si están bien definidos los índices
for beta in np.arange(N)
j = np.arange((beta - 1)*l + 1, beta*l + 1) #hay que revisar si están bien definidos los índices
for gamma in np.arange(N)
k = np.arange((gamma - 1)*l + 1, gamma*l + 1) #hay que revisar si están bien los índices
C[i][:,j]+= A[i][:,k]*B[k][:,j]Ejercicio:
Implementar multiplicación por bloques en Python y R. Perfilar tal implementación.
In [34]:
import os
In [37]:
os.remove(fileA)
os.remove(fileB)
os.remove(fileA_10_3)
os.remove(fileB_10_3)
Referencias:
E. Anderson, Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. Mckenney and D. Sorensen, LAPACK Users Guide, Society for Industrial and Applied Mathematics, Philadelphia, PA, third ed., 1999.
G. H. Golub, C. F. Van Loan, Matrix Computations, John Hopkins University Press, 2013.
Otras referencias útiles para perf son:
perf-tools para tools que se apoyan de perf.