In [134]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from IPython.display import Image # Esto es para desplegar imágenes en la libreta
La regresión softmax (o tambien conocida como regresión logística multinomial) es el último de los algoritmos basado en modelos lineales generalizados que cubriremos en el curso de reconocimiento de patrones. Para ejemplificar su uso, vamos a utilizar una base de datos bastante comun, MNIST.
MNIST es una base de datos de digitos escritos a mano, en formato de $20 \times 20$ pixeles. La base completa puede obtenerse en la página de Yan LeCun (http://yann.lecun.com/exdb/mnist/).
Nosotros en realidad vamos a utilizar una base de datos reducida de la original y con imágenes de calidad más reducida ($16 \times 16$ pixeles por imagen). Numpy provée un método para guardad objetos tipo numpy en un solo archivo, utilizando el método de compresión gunzip. Los datos ya se encuentran preprocesados y empaquetados en un archivo llamado digitos.npz.
In [135]:
data = np.load("digitos.npz")
print data.keys()
En este caso, data es un objeto contenedor de numpy cuyas llaves son X_valida, X_entrena, T_valida, T_entrena. Cada una de estas son a su vez objetos tipo ndarray de numpy, los cuales contienen valores de entrada y salida, tantopara entrenamiento como para validación. No se preocupen, esto de entrenamiento y validación lo vamos a ver más adelante en la clase.
Cada renglon de x es una imagen desenrrollada, esto es los 256 datos de una imágen de $16 \times 16$ pixeles. Por otra parte, cada renglon de y es un vector de 10 posiciones, donde todos los valores son ceros, salvo uno, que es el que define la clase de la imagen.
Para darse una mejor idea, ejecuta el siguiente script varias veces.
In [136]:
x = data['X_entrena']
y = data['T_entrena']
a = np.random.randint(0, y.shape[0])
print "-- x es de dimensiones ", x.shape
print "-- y es de dimensiones ", y.shape
print "\ny si escogemos la imagen ", a, "veremos"
plt.imshow(x[a,:].reshape(16,16), cmap=plt.gray())
plt.axis('off')
plt.show()
print "la cual es un", list(y[a,:]).index(1)
print"\n\nY si miramos lo que contiene, veremos que"
print"x[a,:] = "
print x[a,:]
print "y[a,:] = "
print y[a,:]
O bien, ejecuta este script para ver un grupo grande de imágenes (puedes hacer más grande la imagen para verla mejor).
In [137]:
indices = np.arange(y.shape[0])
np.random.shuffle(indices)
ind = indices[0:100].reshape(10,10)
imagen = np.ones((10 * 16 + 4*11, 10 * 16 + 4*11))
for i in range(10):
for j in range(10):
imagen[4 + i * 20: 20 + i * 20, 4 + j * 20: 20 + j * 20] = x[ind[i, j], :].reshape(16,16)
plt.imshow(imagen, cmap=plt.gray())
plt.axis('off')
plt.title(u"Ejemplos aleatorios de imágenes a clasificar")
Out[137]:
Por último, vamos a hacer una función que agregue la hilera de unos para la x extendida, y simplificar más adelante su uso.
In [138]:
def extendida(x):
"""
Agrega una columna de unos a x
"""
return np.c_[np.ones((x.shape[0], 1)), x]
En la regresión softmax en lugar de estimar solamente un vector de parámetros, lo que tenemos que hacer es estimar una matriz de parámetros $\theta$ tal que:
$$ \theta = (\theta^{(1)}, \ldots, \theta^{(K)}) $$donde $\theta^{(k)} = (\theta_0^{(k)}, \ldots, \theta_n^{(k)})^T$ es el vector columna que parametriza la clase $k$. De esta manera, $\theta$ es ahora una matriz de dimensión $n+1 \times K$. El aporte lineal a cada clase de un objeto $x^(i)$ está dado por
$$ z^{i} = \theta^T x^{(i)}, $$el cual es de dimensión $K \times 1$ (un valor por cada clase). La probabilidad de pertenecer a la clase $k$ está dada por:
$$ \hat{y}_k^{(i)} = softmax_k(z^{(i)}) = \frac{\exp(z_k^{(i)})}{\sum_{r=1}^K \exp(z_r^{(i)})}. $$Como se puede ver, $\hat{y}^{(i)}$ en realizad es el vector de las exponenciales de $z^{(i)}$ normalizadas.
In [139]:
def softmax(z):
"""
Calculo de la regresión softmax
@param z: ndarray de dimensión (T, K) donde z[i, :] es el vector de aportes lineales de el objeto i
@return: un ndarray de dimensión (T, K) donde cada columna es el calculo softmax de su respectivo vector de entrada.
"""
#--------------------------------------------------------------------------------
e = np.exp(z)
return e/e.sum(axis = 1).reshape(-1,1)
#--------------------------------------------------------------------------------
def test_softmax():
z = np.array([[1, 0.1, -10], [20, 100, 2.5], [2, 2.1, 2.2]])
y_hat = softmax(z)
y_selec = y_hat.argmax(axis=1)
y_hat_acc = y_hat.sum(axis=1)
assert abs(y_hat_acc[0] - 1) < 1e-10
assert abs(y_hat_acc[1] - 1) < 1e-10
assert abs(y_hat_acc[2] - 1) < 1e-10
assert y_selec[0] == 0
assert y_selec[2] == 2
return "Paso la prueba"
print test_softmax()
Y ahora es necesario implementar la función de costo, la cual es:
$$ J(\theta) = -\frac{1}{T}\sum_{i=1}^T \sum_{k=1}^K y_k^{(i)} \log(\hat{y}_k^{(i)}), $$donde $y_k^{(i)}$ es un valor de 0 o 1 dependiendo si el objeto $i$ pertenece a la clase $k$ o no, mientras que $\hat{y}_k^{(i)}$ es la probabilidad que el objeto $i$ pertenezca a la clase $k$ conociendo $x^{(i)}$ y parametrizado por $\theta$. Esto se vio en clase y no se realizará el repaso aqui en la libreta.
In [140]:
def costo(theta, x, y):
"""
Calcula el costo para la regresión softmax parametrizada por theta,
con el conjunto de datos dado por (x, y)
@param theta: ndarray de dimensión (n+1, K) con los parámetros
@param x: ndarray de dimensión (T, n+1) con los datos
@param y: ndarray de dimensión (T, K) con la clase por cada dato
@return: Un valor flotante
"""
T, K = y.shape
n = x.shape[1] - 1
#--------------------------------------------------------------------------------
# codigo del profe
# Y_hat = softmax(np.dot(x,theta))
# selecY = Y_hat[y > 0.5]
# return (- np.log(selecY + 0.0001).sum() / T)
#codigo mio
y_estimada = softmax(np.dot(x,theta))
mtz = y * np.log(y_estimada)
return -1.0/T * mtz.sum()
#--------------------------------------------------------------------------------
def test_costo():
x = np.array([[1, 0, 0], [1, 0, 1], [1, 1, 0], [1, 1, 1]])
y = np.eye(4)
theta = np.array([[3, -4, -4],[-1, -1, 3], [.01, 3, -10], [-5, 5, 5]]).T
assert costo(theta, x, y) < 0.1
return "Paso la prueba"
print test_costo()
In [141]:
def predice(theta, x):
"""
Prediccion de y_hat a partir de la matriz theta para los valores de x
@param theta: ndarray de dimensión (n+1, K) con los parámetros
@param x: ndarray de dimensión (T, n+1) con los datos
@return: ndarray de dimensión (T, K) con la clase por cada dato (unos y ceros)
"""
#--------------------------------------------------------------------------------
h = softmax(np.dot(x,theta))
y_estimada = np.zeros_like(h)
indices_argmax = h.argmax(axis=1)
j = 0
for i in indices_argmax:
y_estimada[j,i]=1
j+=1
# print y_estimada
return y_estimada
#--------------------------------------------------------------------------------
def prueba_prediccion():
x = np.array([[1, 0, 0], [1, 0, 1], [1, 1, 0], [1, 1, 1]])
y = np.eye(4)
theta = np.array([[3, -4, -4],[-1, -1, 3], [.01, 3, -10], [-5, 5, 5]]).T
assert abs((y - predice(theta, x)).sum()) < 1e-12
return "Paso la prueba"
print prueba_prediccion()
Y por último tenemos que implementar el gradiente para poder utilizar los métodos de optimización (ya sea por descenso de gradiente o por algún método de optimización.
El gradiente se obtiene a partir de las derivadas parciales:
$$ \frac{\partial J(\theta)}{\partial \theta_s^{(r)}} = - \frac{1}{T} \sum_{i = 1}^T \left(y_k^{(i)} - \sum_{j=1}^K \hat{y}_j^{i}\right) x_s^{(i)}, $$donde $k$ es la clase a la que pertenece el ejemplo $i$. Esto implica que cuando $k = r$, entonces solamente tenemos $1 - \hat{y}_k^{(i)}$ para este dato. Si $k \neq r$ entonces tendremos $-\hat{y}_r^{(i)}$.
Esto se puede resolver en forma matricial como
$$ \nabla J(\theta) = - \frac{1}{T} X^T (Y - \hat{Y}) $$
In [142]:
def gradiente(theta, x, y):
"""
Calculo del gradiente para el problema de regresión softmax
@param theta: ndarray de dimensión (n+1, K) con los parámetros
@param x: ndarray de dimensión (T, n+1) con los datos
@param y: ndarray de dimensión (T, K) con la clase por cada dato
@return: Un ndarray de mismas dimensiones que theta
"""
#--------------------------------------------------------------------------------
T = y.shape[0]
y_estimada = softmax(np.dot(x,theta))
return - np.dot(x.T, (y-y_estimada)) / T
#--------------------------------------------------------------------------------
def prueba_gradiente():
x = np.array([[1, 0, 0], [1, 0, 1], [1, 1, 0], [1, 1, 1]])
y = np.eye(4)
theta = np.array([[3, -4, -4],[-1, -1, 3], [.01, 3, -10], [-5, 5, 5]]).T
g = gradiente(theta, x, y)
assert np.abs(g).max() < 0.05
return "Paso la prueba"
print prueba_gradiente()
Ahora si, ya nos encontramos en capacidad para realizar el aprendizaje para la regresión logística.
In [143]:
def dg_softmax_lotes(theta, x, y, alpha=None, max_epoch=10000, epsilon=1e-3, errores=False):
"""
Descenso de gradiente por lotes para la clasificación softmax
"""
historial = np.zeros((max_epoch)) if errores else None
for epoch in xrange(max_epoch):
#--------------------------------------------------------------------------------
theta = theta - alpha * (gradiente(theta, x, y))
if errores:
historial[epoch] = costo(theta, x, y)
if np.max(np.abs(gradiente(theta, x, y))) <= epsilon:
return theta, historial
#--------------------------------------------------------------------------------
return theta, historial
Pero para utilizar el descenso de gradiente hay que ajustar un valor de $\alpha$.
In [144]:
# Ajusta un valor de alpha razonable
alfita = 1.5
T, K = y.shape
n = x.shape[1]
theta = 0.1 * (np.random.random((n + 1, K)) - 0.5)
theta, e_hist = dg_softmax_lotes(theta, extendida(x), y, alpha=alfita, max_epoch=100, errores=True)
plt.plot(e_hist)
Out[144]:
y para probarlo vamos a aprender a clasificar a los digitos de nuestra base de datos
In [145]:
theta = 0.1 * (np.random.random((n + 1, K)) - 0.5)
theta, e_hist = dg_softmax_lotes(theta, extendida(x), y, alpha=alfita, max_epoch=1000)
print "El costo de la solución final es de ", costo(theta, extendida(x), y)
y_estimada = predice(theta, extendida(x))
errores = np.where(y.argmax(axis=1) == y_estimada.argmax(axis=1), 0, 1)
print "\nLos datos utilizados para el aprendizaje mal clasificados son el ", 100 * errores.mean(),"%"
# Esto solo es para hacerla más emocionante
x_test = data['X_valida']
y_test = data['T_valida']
y_estimada_T = predice(theta, extendida(x_test))
errores = np.where(y_test.argmax(axis=1) == y_estimada_T.argmax(axis=1), 0, 1)
print "\nY con los datos de prueba el error es del ", 100 * errores.mean(),"%"
¿Será esta la mejor solución? ¿Será una buena solución? Por esto no hay que preocuparse mucho todavía, lo vamos a revisar más adelante en el curso. Se espera con la regresión softmax poder clasificar correctamente más del 95% de los datos de entrenamiento. Bueno lo que nos falta ahora es agregar al aprendizaje tipo softmax la regularización.
In [146]:
def costo_reg(theta, x, y, lammbda):
"""
Calcula el costo para la regresión softmax parametrizada por theta,
con el conjunto de datos dado por (x, y)
@param theta: ndarray de dimensión (n+1, K) con los parámetros
@param x: ndarray de dimensión (T, n+1) con los datos
@param y: ndarray de dimensión (T, K) con la clase por cada dato
@param lammbda: Valor de la regularización
@return: Un valor flotante
"""
#--------------------------------------------------------------------------------
# AGREGA AQUI TU CÓDIGO
#--------------------------------------------------------------------------------
T = x.shape[0]
return costo(theta, x, y) + (lammbda/(2.0*T)) * np.square(theta[1:,:]).sum()
#--------------------------------------------------------------------------------
def gradiente_reg(theta, x, y, lammbda):
"""
Calculo del gradiente para el problema de regresión softmax
@param theta: ndarray de dimensión (n+1, K) con los parámetros
@param x: ndarray de dimensión (T, n+1) con los datos
@param y: ndarray de dimensión (T, K) con la clase por cada dato
@param lammbda: Valor de la regularización
@return: Un ndarray de mismas dimensiones que theta
"""
#--------------------------------------------------------------------------------
# AGREGA AQUI TU CÓDIGO
#--------------------------------------------------------------------------------
T = x.shape[0]
hipotesis = softmax(np.dot(x,theta))
identidad = np.eye(theta.shape[0])
identidad[0,0] = 0
return (-1.0/T) * (np.dot(x.T, y-hipotesis)) + (lammbda/T * identidad.dot(theta))
#--------------------------------------------------------------------------------
def dg_softmax_lotes_reg(theta, x, y, lammbda, alpha=None, max_epoch=10000, epsilon=1e-3, errores=False):
"""
Descenso de gradiente por lotes para la clasificación softmax con regularización
"""
#--------------------------------------------------------------------------------
# AGREGA AQUI TU CÓDIGO
#--------------------------------------------------------------------------------
historial = np.zeros((max_epoch)) if errores else None
for epoch in xrange(max_epoch):
#--------------------------------------------------------------------------------
theta = theta - alpha * (gradiente_reg(theta, x, y, lammbda))
if errores:
historial[epoch] = costo_reg(theta, x, y, lammbda)
if np.max(np.abs(gradiente_reg(theta, x, y, lammbda))) <= epsilon:
return theta, historial
#--------------------------------------------------------------------------------
return theta, historial
#--------------------------------------------------------------------------------
Y para probar:
In [147]:
theta5 = 0.1 * (np.random.random((n + 1, K)) - 0.5)
theta5, e_hist = dg_softmax_lotes_reg(theta5, extendida(x), y, lammbda=5, alpha=1, max_epoch=1000)
print "Para lammbda = 5"
print "-----------------"
print "\nEl costo de la solución final es de ", costo(theta, extendida(x), y)
y_estimada = predice(theta5, extendida(x))
errores = np.where(y.argmax(axis=1) == y_estimada.argmax(axis=1), 0, 1)
print "\nLos datos utilizados para el aprendizaje mal clasificados son el ", 100 * errores.mean(),"%"
# Esto solo es para hacerla más emocionante
x_test = data['X_valida']
y_test = data['T_valida']
y_estimada_T = predice(theta5, extendida(x_test))
errores = np.where(y_test.argmax(axis=1) == y_estimada_T.argmax(axis=1), 0, 1)
print "\nY con los datos de pureba el error es del ", 100 * errores.mean(),"%\n\n"
#----------------------------------------------------------------------------------------------------------
theta20 = 0.1 * (np.random.random((n + 1, K)) - 0.5)
theta20, e_hist = dg_softmax_lotes_reg(theta20, extendida(x), y, lammbda=20, alpha=1, max_epoch=1000)
print "Para lammbda = 20"
print "-----------------"
print "\nEl costo de la solución final es de ", costo(theta, extendida(x), y)
y_estimada = predice(theta20, extendida(x))
errores = np.where(y.argmax(axis=1) == y_estimada.argmax(axis=1), 0, 1)
print "\nLos datos utilizados para el aprendizaje mal clasificados son el ", 100 * errores.mean(),"%"
# Esto solo es para hacerla más emocionante
x_test = data['X_valida']
y_test = data['T_valida']
y_estimada_T = predice(theta20, extendida(x_test))
errores = np.where(y_test.argmax(axis=1) == y_estimada_T.argmax(axis=1), 0, 1)
print "\nY con los datos de pureba el error es del ", 100 * errores.mean(),"%"
Con la regularización y $\lambda = 5$ se espera que el error en los datos de prueba sea ligeramente mejor que sin regularización, mientras que con $\lambda = 20$ será ligeramente peor. Obviamente, con los datos de aprendizaje no se puede hacer mejor que sin regularización.
¿Y los valores de $\theta$ como cambian con la regularización?
In [ ]: