Ahora vamos a usar un algoritmo llamado regresion logística para poder clasificar objetos o fenómenos en dos clases diferentes. Esto lo lograremos tomando un conjunto de entrenamiento que está marcado como clase uno o como clase dos y lo probaremos sobre datos que no conocemos con antelación. De esta manera podemos ver si nuestro modelo funciona sobre estos datos. Los resultados los comparamos con las predicciones hasta el final. Para ello usaremos regresión logística. Hay que recordar que el modelo de regresión logística únicamente puede separa datos linealmente separables, al igual que el perceptrón

En este artículo vamos a trabajo con numpy y python. Así que si gustan, de entrada debemos importar:

import numpy as np 

Notación y problema

En este ejemplo vamos a hacer la clasificación sobre el problema: ¿un dato nuevo petenece a la clase 1 o a la case dos?. Si x es el valor que se desea obtener en la clasificación, entonces:

\[ y= \begin{cases} 1 & \text{Pertenece a clase 1} \\ 0 & \text{Pertenece a clase 2} \end{cases} \]

Ahora que hemos definido las clases, debemos tomar en cuenta representar de forma matemática las características de la foto. Vamos a representar en forma de vector los pixeles en una división RGB de la imagen y la vamos a concatenar. Esto nos arrojará un gigantesco vector que contendrá todos los valores de los pixeles de la imagen en una dimensión. Con estos valores definamos nuestra notación.

Cada ejemplo es un par ordenado

\[ (x,y), x\in \mathbb{R}^{n_x}, y\in\{0,1\} \]

 donde x es un vector n dimensional de características y y es uno o cero. Para nuestro conjunto de entrenamiento se definirán m ejemplos entrenamientos tal que onbtendremos los siguientes ejemplos de entrenamiento:

\[M=\{(x^1,y^1), (x^2,y^2), \dots, (x^m,y^m)\} \]

Para facilitar el trabajo y recordando que cada x es un vector, podemos unir estos vectores en una matriz, que llamaremos X, y que tiene la forma:

\[X=[x^1,x^2,\dots,x^m], \text{ donde } X\in\mathbb{R}^{n_x\times m}\]

Ahora, las clases también se guardarán en un vector. Si recordamos, que cada x tiene asociado una y escalar, entonces el vector de todos los y es:

\[Y=[y^1,y^2,\dots,y^m], \text{ donde } Y\in\mathbb{R}^{1\times m}\]

Vamos a crear un set artificial de datos para pode mostrar lo expuesto arriba. Este set será usado para exponer todo el trabajo posterior que vamos a realizar.

X = np.matrix(
        [[1.,3.,2.,2.,0.],
        [2.,0.,3.,1.,1.],
        [3.,2.,2.,1.,2.],
        [0.,1.,2.,8.,2.],
        [2.,4.,3.,2.,4.],
        [3.,1.,4.,0.,3.],
        [4.,3.,1.,2.,2.],
        [0.,2.,1.,1.,3.],

        [8.,7., -2.,12.,5.],
        [7.,5.,  0.,21.,4.],
        [8.,0.,-12.,18.,4.],
        [9.,5., -2.,18.,6.],
        [6.,4.,  0.,22.,5.],
        [7.,1., -1.,11.,3.],
        [8.,5.,  1.,23.,6.],
        [3.,5., -3.,12.,4.]]
    )

Y = np.matrix([[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]])

 

Ahora que hemos definido la notación que vamos a usar, es hora de plantear el problema. Nosotros queremos encontrar una función que data una entrada x nos pueda generar una estimación de y, tal que:

 

\[\hat{y}=P(y=1|x) \text{ y los parámetros } w\in\mathbb{R}^{n_x}, b\in\mathbb{R}\]
 

La Regresión Logística

 

 ¿Pero qué función logrará definir mejor el comportamiento de nuestra salida para definir la probabilidad de y? Es posible generar una función lineal con la multiplicación de wTx+b , pero esto tendría algunos problemas. El primero de ellos es que al ser lineal puede adquirir valores negativos y positivos superiores a 1. Pero al ser una probabilidad la que necesitamos regresar, esto no tendría sentido. La función que necesitamos generar debe encontrarse entre 0<=y<=0, por lo que usaremos una función sigmoide. Esta se comporta de la siguiente manera:

Como se puede ver nunca llega a tocar el eje de las y a 0 o 1, y se mueve entre estos dos valores. Lo cual la hace perfecta para expresar una probabilidad para nuestra decisión binaria. Veamos ahora cómo es la forma de una función sigmoide:

\[S(z) = \frac{1}{1+e^{-z}} \text{ donde } z=w^Tx+b\]
 

Como ya tenemos la función que nos permitirá clasificar correctamente nuestras fotografías, nuestro modelo debe encontrar los valores correctos de w y b que nos permitan encontrar el valor correcto de la función S. Usando la biblioteca numpy la implementación de nuestra función logística quedaría cómo:

def sigmoid(z):
    return 1/(1+np.exp(-z))

El vector de pesos w se inicializa de manera arbitraria. En este caso lo iniciaremos con ceros, pero también se puede hacer de manera aleatoria.

W = np.zeros(shape=(1,X.shape[1]))

La función de costo!

Ahora bien, en nuestros estimadores W y b se escogen en un inicio de manera aleatoria, pero esto no parece ser muy útil. ¿Cómo saber qué tan bien o mal logran predecir nuestros estimadoras la clase de nuestra foto? Para ellos requerimos una función de error o loss function. Definir es una función así no es trivial. Por ejemplo, hay que tomar en cuenta que nuestro problema no es convexo.  Para este problema se utilizará la función de error:

\[\mathscr{L}(\hat{y},y) = - \big(y \log \hat{y} + (l-y) \log (1-\hat{y})\big)\]

 Esta función únicamente se usará sobre una entrada del conjunto de datos. Ahora, para obtener una imagen general sobre el error de los parámetros con respecto al total de los datos, usaremos la función de costo, que es definida cómo:

\[\begin{align}J(W,b)&=\frac{1}{m}\sum^m_{i=1}\mathscr{L}(\hat{y}^i,y^i)\\ &=-\frac{1}{m}\sum^m_{i=1}\big[ y^i \log \hat{y}^i + (1-y^i) \log (1-\hat{y}^i)\big] \end{align}\]

El objetivo de nuestra red neuronal va a ser minimizar el costo total del modelo, para así acercarnos lo más posible a nuestras clases dadas y. Si logramos ajustar los parámetros de tal forma que predigan y, entonces estaremos entrenando un modelo bastante bueno para el uso con otros datos no vistos. ¿Cómo hacerlo? Lo vemos a continuación. 

El algoritmo de la gradiente descendiente. 

Ya que hemos explicado bien nuestra función de predicción y nuestra función de costo queda claro el problema para resolver el modelo: buscamos los parámetros w y b que minimicen J(w,b) . Lo que queremos lograr es descender en nuestro mundo hasta un mínimo que esperemos sea el mínimo global. El espacio definido por la función de costo y el cual queremos descender tendrá una figura como a continuación se muestra:

 

 Tendremos dos ejes donde theta 1 representa los valores de W y theta 0 a los valores de b, y la altura es el valor que adquiere la función de costo J con esta combinación de valores. La función de costo que se escogió es buena para este problema por que es convexa y esto nos ayudará a encontrar la solución global al problema. 

¿Cuál es la idea detrás de este algoritmo? La idea central es en cada paso del algoritmo ir descendiendo por el relieve de la función. Sabremos cuantos pasos y en qué dirección descender por la derivada de J respecto a w y b.  Esto se expresa en las siguientes ecuaciones:

\[\begin{align}w&=w-\alpha \frac{\partial J(w,b)}{\partial w}\\ b&=b-\frac{\partial J(w,b)}{\partial b}\end{align}\]

 Cada una de estas actualizaciones se realizará en cada ciclo, hasta llegar a la convergencia. La derivada de J(w,b) en cada función es la pendiente de la función de costo, y nos indicará en qué dirección bajar y cuanto. El alfa que multiplica la derivada es un valor dado que permite acelerar o decrecer esos pasos. Es útil poder modificar este valor, ya que para cada problema vamos a necesitar diferentes velocidades de descenso. Si bien este tema es complejo lo hemos abordado de una manera sencilla. Si quieren más información quiero recomendar leer el artículo sobre la gradiente descendiente aquí

Gradiente descendiente de la función logística. 

Ahora bien, vamos a calcular la gradiente descendiente para nuestra función logistica. Recordemos las formulas para la regresión logística, donde y estimada es una sigmoide de z, y L es la función de pérdida. 

\[ \begin{align} z&=w^Tx+b\\ \hat{y}&=a=S(z)\\ {\cal L}(a,y)&=-(y \log (a)+(1-y) \log (1-a)) \end{align} \]

 Si imaginamos dos características sobre las que vamos (w1 y w2) a trabajar, entonces tendremos que las operaciones de propagación serían de la siguiente forma:

\[ \begin{align} z=w_1x_1 + w_2x_2+b \to a=S(z) \to {\cal L}(a,y) \end{align} \]

Ahora se intentará hacer la propagación para atrás, reduciendo la función de pérdida L(a,y) y encontrar los valores de w1 y w2 que lo logran. Para lograrlo se calculará la regla la derivada de z respecto a L(a,y). Pero dado a que L(a,y) es derivable respecto a a, entonces por regla de la cadena obtendremos: 

\[ \begin{align} \frac{\partial {\cal L}}{z} &= \frac{\partial a}{z} \frac{\partial {\cal L}}{a} \\ &=a(1-a) \Big( -\frac{y}{a}+\frac{1-y}{1-a} \Big) \\ &=a-y \end{align} \]

Con la derivada calculada, únicamente necesitamos cuánto vamos a cambiar w1, w2 y b. Primer calculamos las derivadas parciales dz respecto a todas las variables. Y una vez obtenido esto, lo multiplicamos por alfa, que es el factor de aprendizaje (una constante definida por nosotros de manera arbitraria) y lo restamos a las respectivas variables. 

 

\[ \begin{align} \frac{\partial {\cal L}}{w_1}&=x_1\frac{\partial {\cal L}}{z} \\ \frac{\partial {\cal L}}{w_2}&=x_2\frac{\partial {\cal L}}{z} \\ \frac{\partial {\cal L}}{b}&=\frac{\partial {\cal L}}{z} \\ w_1 &= w_1 - \alpha \frac{\partial {\cal L}}{w_1}\\ w_1 &= w_1 - \alpha \frac{\partial {\cal L}}{w_2}\\ b &= b - \alpha \frac{\partial {\cal L}}{b}\\ \end{align} \]

Vectorización de la gradiente descendiente de la regresión logística.

 Para pasar a código lo explicado en el artículo pasado, vamos a usar la biblioteca numpy que servirá para mejorar la velocidad de cálculo usando operaciones matriciales  evitando resolver una por una las variables. Para obtener el valor de todas las z en el vector Z podemos usar:

Z = np.dot(w.T,X)+b

Donde w es una matriz que contiene en sus m columnas (una columna por ejemplo de entrenamiento) y una matriz w, que contiene los pesos en cada columna (una columna por ejemplo). b es un escalar. Para sumar un escalar con un vector (Z es un vector resultado de la multiplicación de matrices) numpy genera un vector del escalar; esta operación se llama Broadcasting. Las matrices que usamos tienen la siguiente forma, tal cuál se han descrito:

\[ X= \begin{bmatrix} \vdots & \vdots & & \vdots \\ x_1 & x_2 & \dots  & x_m \\ \vdots & \vdots & & \vdots \\ \end{bmatrix} \]
 

\[ W^T= \begin{bmatrix} \vdots & \vdots & & \vdots \\ w_1 & w_2 & \dots  & w_m \\ \vdots & \vdots & & \vdots \\ \end{bmatrix} \]
 

La operación descrita en el código es el producto de las dos matrices por el vector b que se ha expandido del escalar b, tal que:

\[ \begin{bmatrix} z_1, z_2, \dots, z_m \end{bmatrix} = w^T X^T \begin{bmatrix} b,b, \dots, b \end{bmatrix} \]
 

De igual manera, podemos calcular de manera vectorial todas las operaciones de la gradiente descendiente. Definimos los valores a y z cómo vectores:

\[ A= \begin{bmatrix} a_1, a_2, \dots, a_m \end{bmatrix} \\ Y= \begin{bmatrix} y_1,y_2, \dots, y_m \end{bmatrix} \\  \frac{\partial {\cal L}}{Z} = A-Y \]
 

 Ahora, para calcular la derivada de b y w, vamos a usar también operaciones vectoriales:

\[ \frac{\partial {\cal L}}{b} = \frac{1}{m} \sum^m_{i=1}\frac{\partial {\cal L}}{z_i}\]
 

Que expresado en código python es:

db = (1/m)*np.sum(dz)

Y dw cómo:

\[ \partial w = \frac{1}{m} X \frac{\partial {\cal L}}{Z}^T\]
  

Para resumir, cada iteración de gradiente descendiente vectorizada tendrá la siguiente forma:

\[ \begin{align} Z &= w^T X+b \\ A &= S(Z) \\ \frac{\partial {\cal L}}{Z} &= A-Y \\ \frac{\partial {\cal L}}{w} &= \frac{1}{m} X \frac{\partial {\cal L}}{Z}^T\\ \frac{\partial {\cal L}}{b} &= \frac{1}{m} \sum^m_{i=1}\frac{\partial {\cal L}}{z_i} \\ w&=w-\alpha \frac{\partial {\cal L}}{w} \\ b&=b-\alpha \frac{\partial {\cal L}}{b} \end{align} \]

La gran ventaja de este método es que aumenta significativamente la velocidad de cálculo y permite distribuir el trabajo sobre varias CPU's y GPU's.

Lo anterior es posible de la sigueinte forma en python:

def train(X,Y,alpha=.1,iter=100):
    W = np.zeros(shape=(1,X.shape[1]))
    b = 1
    for a in range(iter):
        
        # Feed Forward
        z = np.matmul(X,W.T) + b
        a = Yhat = sigmoid(z.T)
        m = a.shape[1]

        # Loss calculation
        loss = -(np.multiply(Y,np.log(Yhat))+np.multiply((1-Y),np.log(1-Yhat)))
        loss = (1/m)*np.sum(loss)
        print(loss)
        
        # Gradient Decent
        # Calculate derivatives
        dz = a-Y
        dw = (1/m)*(X.T@dz.T).T
        db = (1/m)*np.sum(dz)

        
        # Apply changes
        W = W - alpha*dw
        b = b - alpha*db
    return {'W':W,'b':b}

El la iteración se realizará n veces, para reducir el costo del modelo. El resultado es una reducción del costo del modelo. Veamos el comportamiento sobre 10 iteraiciones:

0.8132616875182228
0.7172191897298279
0.4655356015304304
0.311442095470515
0.22082850671820856
0.16998968532821118
0.14528730059496492
0.13459595441291947
0.1282303305476174
0.12281237720910451

Entre menor es el costo, más se ajustan los valores de W y b para poder realizar una predicción acertada. 

Implementación

 A continuación vamos a presentar el código completo, para que lo prueben y modifiquen en casa. 

import numpy as np

# Author: Manuel Mager
# Creamos un data set artificial X = np.matrix( [[1.,3.,2.,2.,0.], [2.,0.,3.,1.,1.], [3.,2.,2.,1.,2.], [0.,1.,2.,8.,2.], [2.,4.,3.,2.,4.], [3.,1.,4.,0.,3.], [4.,3.,1.,2.,2.], [0.,2.,1.,1.,3.], [8.,7., -2.,12.,5.], [7.,5., 0.,21.,4.], [8.,0.,-12.,18.,4.], [9.,5., -2.,18.,6.], [6.,4., 0.,22.,5.], [7.,1., -1.,11.,3.], [8.,5., 1.,23.,6.], [3.,5., -3.,12.,4.]] ) Y = np.matrix([[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]]) def sigmoid(z): return 1/(1+np.exp(-z))

def train(X,Y,alpha=.1,iter=100): W = np.zeros(shape=(1,X.shape[1])) b = 1 for a in range(iter): # Feed Forward z = np.matmul(X,W.T) + b a = Yhat = sigmoid(z.T) m = a.shape[1] # Loss calculation loss = -(np.multiply(Y,np.log(Yhat))+np.multiply((1-Y),np.log(1-Yhat))) loss = (1/m)*np.sum(loss) print(loss) # Gradient Decent # Calculate derivatives dz = a-Y dw = (1/m)*(X.T@dz.T).T db = (1/m)*np.sum(dz) # Apply changes W = W - alpha*dw b = b - alpha*db return {'W':W,'b':b}

def predict(X, model): W = model['W'] b = model['b'] z = np.matmul(X,W.T) + b a = Yhat = sigmoid(z.T) m = a.shape[1] print(a) for ai in np.matrix.tolist(a)[0]: print(ai) if ai>.5: print('We predict', 1) else: print('We predict', 0)

if __name__ == "__main__": model = train(X,Y) Xpred =[[12.,8.,-4.,22.,5.], [7. ,9.,-2.,8.,2.], [-1.,2., 2.,0.,1.], [2. ,3., 3.,2.,5.]] predict(Xpred, model)

Agradecemos cualquier comentario. Las traspuetas fueron adaptadas para poder funcionar mejor con los datos ingresados. Para realizar la preducción simplemente se realizar el paso de forward propagatin y se estima sobre los resultados de las sigmoides. 

Esperamos que este tutorial les ha gustado. 

 

Share This