Broadcasting en numpy

Mariano Rivera

version: enero 2018

Broadcasting es una estrategia de numpy para operar con arreglos mutidimensionles de diferentes dimensiones.
Por ejemplo, suponga que desamos multiplicar en escalar $a$ por vector $x \in \mathbb{R}^n$:

$y = a x$

Esto equivale, matemáticamente, a:

$y_i = a x_i, \;\;\text{para}\,\; i=1,2,\ldots,n$
Es decir cada elemento del vector $x$ es multiplicado por el escalar $a$. Es como si, en términos computacionales, se creara un vector virtual con réplicas del escalar:

$\begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} = \begin{bmatrix} a \\ a \\ \vdots \\ a \end{bmatrix} \odot \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} = \begin{bmatrix} a y_1 \\ a y_2 \\ \vdots \\ a y_n \end{bmatrix}$

donde $\odot$ denota el producto de Hadamard (elemento a elemento), no confundir con el producto punto (dot) o producto Tensorial

Recordemos que producto de Hadamard es el default en numpy. Es decir

El producto de Hadamard entre dos vectores se calcula mediante:

import numpy as np
# arreglo de 3 números enteros aleatorios distribucion uniforme [low, high)
x = np.random.randint(low=1,high=10, size=3)  
y = np.random.randint(low=1,high=10, size=3)

print(x)
print(y)
print(x*y)  # producto de Hadamard

[2 7 7]
[2 7 8]
[ 4 49 56]

Producto punto:
$x \cdot y = \sum_i x_i y_i = \sum_i [x \odot y]_i$
Esto es:

print(x.dot(y))
print(sum(x*y))  

61
61

Regresemos al ejemplo que nos motivó, el producto de un escalar por un vector que es un ejemplo de broadcasting. Antes de explicar este proceso, explicaremos la convención usada por PYTHON y extendida en cómputo Tensorial que es ampliamente usado en implementaciones de Redes Neuronales de la modalidad de Aprendizaje Profundo (Deep Learning):

a) Una arreglo multidimensiosional tienen ejes (axes). Es decir, un vector tienen 1 eje, una matriz 2 ejes, y asi.

b) El tamaño de los ejes son la forma (shape). Asi una Matriz de $m$ renglones por $n$ columnas es, en la convención de Python, un Tensor de 2 ejes (tensor de orden 2) con forma de (m,n).

Ahora si, el broadcasting es un proceso automático implementado en las operaciones de numpy para eficientar los cálculos con arreglos multidimensionales donde, siempre que no haya ambiguedad, se procede como sigue:

• De ser necesario, toma el arreglo de menor forma (shape) y agrega 1’s por la izquierda a la forma para hacerlos compatibles.

• De ser necesario, replica el arreglo en los ejes (axes) con forma igual a 1, hasta que alcance la forma del arreglo de mayor forma.

Por ejemplo sumar X+Y

1. Asumamos X tiene shape ( 32, 256, 256) y Y tiene shape ( 30, 2 , 32, 256, 256)

2. Extiende los axes de X a (1,1, 32, 256, 256)

3. Replica X hasta alcanzar la shape de Y: ( 30, 2 , 32, 256, 256)

4. Suma elemento a elemento X+Y y el resultado tienen shape ( 30, 2 , 32, 256, 256)

El broadcasting se realiza en forma virtual, está implementado en memoria, por lo que es muy eficiente.

Ejemplo

import numpy as np
import pprint

# arreglo de 3 números enteros aleatorios distribucion uniforme [low, high)
X = np.random.randint(low=1,high=10, size=(3,3))  
Y = np.random.randint(low=1,high=10, size=(3,1))
Z = np.random.randint(low=1,high=10, size=(3))

print('X.shape =', X.shape) 
print('Y.shape =', Y.shape) 
print('Z.shape =', Z.shape) 

pp.pprint('matriz de 3x3, X=')
pp.pprint(X)

pp.pprint('vector columna de tamaño 3, Y= ')
pp.pprint(Y)

pp.pprint('vector columna de tamaño 3, Z= ')
pp.pprint(Z)

X.shape = (3, 3)
Y.shape = (3, 1)
Z.shape = (3,)
'matriz de 3x3, X='
array([[4, 1, 2],
       [8, 8, 5],
       [3, 7, 6]])
'vector columna de tamaño 3, Y= '
array([[7],
       [2],
       [1]])
'vector columna de tamaño 3, Z= '
array([7, 2, 3])

---

Calcular
$X+Y$

X es de (3,3), Y es (3,1) por lo que Y se replicará por columnas y se sumará:
$X+ \begin{bmatrix} Y, Y, Y \end{bmatrix}$

print('X+Y=')
pp.pprint(X+Y)  # suma elemento a elemento con Broadcasting

X+Y=
array([[11,  8,  9],
       [10, 10,  7],
       [ 4,  8,  7]])

Al vector Y no le es necesario añadir axes por la derecha (de hecho, numpy no lo hubiera hecho) pues su shape es (3,1).

---

Calcular
$X+Z$
X es de (3,3) y Z es (3). Luego Z se hará de (1,3), se replicará por renglones y se sumará:
$X+ \begin{bmatrix} Z \\ Z \\ Z \end{bmatrix}$

print('X+Z=')
pp.pprint(X+Z)  # suma elemento a elemento con Broadcasting

X+Z=
array([[11,  3,  5],
       [15, 10,  8],
       [10,  9,  9]])

Ejemplo: producto exterior

$Y Z^\top$

Y = np.array(range(0,4))
Y = Y.reshape((Y.shape[0],1))
Zt = np.array(range(0,4))
print('Y=\n', Y)
print('Zt =\n', Z)

Y=
 [[0]
 [1]
 [2]
 [3]]
Zt =
 [0 1 2 3]

Uno es renglón y el otro columna, asi que funciona el bradcasting (esto es correcto, matemáticamente)

Y*Zt  

array([[0, 0, 0, 0],
       [0, 1, 2, 3],
       [0, 2, 4, 6],
       [0, 3, 6, 9]])

Expliquemos que se hizo, paso a paso:

Y.shape = (4,1)

Z.shape = (4)

luego agrega ejes necesarios

Z.shape = (1,4)

y replica para poder sumar

$\begin{bmatrix} Y, Y, Y, Y \end{bmatrix} * \begin{bmatrix} Z \\ Z \\ Z \\ Z \end{bmatrix}$
que son arreglos de $4 \times 4$

Ahora intentemos calcular el producto punto punto invirtiendo los factores $Z^\top Y$:

Zt*Y  

array([[0, 0, 0, 0],
       [0, 1, 2, 3],
       [0, 2, 4, 6],
       [0, 3, 6, 9]])

¿Pero que es esto? es broadcasting, no es algebra lineal!

Debemos ser cuidadosos. Los pasos son muy similares al caso anterior, solo que ahora:

$\begin{bmatrix} Z \\ Z \\ Z \\ Z \end{bmatrix} * \begin{bmatrix} Y, Y, Y, Y \end{bmatrix}$
y el orden de los factores no altera el producto, al menos en este nivel: elemento a elemento.