Sistema de reconocimiento de caracteres manuscritos usando redes neuronales convolucionales implementado en Python

i Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería Electrónica, Robótica y Mecatrónica Sistema de reconocimiento de caracteres manuscritos usando Redes Neuronales Convolucionales implementado en Python Dep. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Autor: Alejandro Trujillo Quevedo Tutor: Alejandro del Real Torres Sevilla, 2017 ii iii Trabajo Fin de Grado Grado en Ingeniería Electrónica, Robótica y Mecatrónica Sistema de reconocimiento de caracteres manuscritos usando Redes Neuronales Convolucionales implementado en Python Autor: Alejandro Trujillo Quevedo Tutor: Alejandro del Real Torres Profesor asociado Dep. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 iv v Trabajo Fin de Grado: Sistema de reconocimiento de caracteres manuscritos usando Redes Neuronales Convolucionales implementado en Python Autor: Alejandro Trujillo Quevedo Tutor: Alejandro del Real Torres El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: vi Sevilla, 2017 El Secretario del Tribunal vii A mi familia A mis maestros viii ix Agradecimientos Me gustaría, antes de todo, agradecer a mi tutor, Alejandro del Real Torres, por la oportunidad brindada durante estos meses, así como la ayuda y el apoyo prestados durante la realización de este proyecto. A mi familia por su inquebrantable confianza y su inestimable apoyo durante todos y cada uno de estos años, en los que me alentaban a no rendirme y seguir adelante. A mis amigos, los de toda la vida y las grandes amistades formadas durante estos años, por el apoyo incondicional y la ayuda prestada cuando más lo necesitaba. xvi 6 Implementación en Python 63 6.1. Librerias 63 6.2. Diferencias con respecto a MATLAB 64 6.3. Soluciones y Optimización 66 6.4. Resultados obtenidos 70 7 Conclusiones 71 8 Líneas de Trabajo Futuras 73 Referencias 75 Apéndice I: Códigos de Regresión Lineal 77 Apéndice II: Códigos de Regresión Logística 81 Apéndice III: Códigos de Regresión Softmax 83 Apéndice IV: Códigos de la Red Neuronal Multicapa 85 Apéndice V: Códigos de la Red Neuronal Convolucional 91 Apéndice VI: Códigos en Python tras Optimización 101 xvii ÍNDICE DE TABLAS Tabla 6–1: Comparación de matrices según el método de indexación 66 Tabla 6–2: Matriz de 4 dimensiones 68 xviii xix ÍNDICE DE FIGURAS Figura 3-1: Valores RMS obtenidos y tiempo de ejecución. 30 Figura 3-2: Gráfica con precios reales y predichos. 31 Figura 3-3: Valores RMS obtenidos y tiempo de ejecución. 32 Figura 3-4: Gráfica con precios reales y predichos. 32 Figura 3-5: Imágenes del set de datos MNIST 34 Figura 3-6: Resultados de la Regresión Logística 34 Figura 3-7: Resultados de la Regresión Softmax 38 Figura 4-1: Esquema de una neurona 39 Figura 4-2: Diagrama de una neurona artificial 40 Figura 4-3: Comparación de funciones de activación 41 Figura 4-4: Red neuronal con varias neuronas 41 Figura 4-5: Red neuronal multicapa 43 Figura 4-6: Resultados obtenidos con Red Neuronal Multicapa 49 Figura 5-1: Convolución 51 Figura 5-2: Convolución completa 52 Figura 5-3: Pooling 53 Figura 5-4: Capa convolucional y de pooling aplicadas a una imagen 56 Figura 5-5: Red Neuronal Convolucional 56 Figura 5-6: Resultados obtenidos con una CNN 61 Figura 6-1: Carga del set MNIST 64 Figura 6-2: Función sub2ind en Python 65 Figura 6-3: Carga de datos del MNIST con librería mnist 66 Figura 6-4: Carga de un archivo .mat 67 Figura 6-5: Función para redimensionar matriz 2D a 4D 67 Figura 6-6: Pooling con bucles anidados 69 Figura 6-7: Almacenamiento de parámetros si la precisión supera el 97% 70 Figura 6-8: Resultados obtenidos en Python 70 xx xxi Notación a Escalar a Vector A Matriz ∈ Perteneciente a ℜ𝑛 Conjuntos de números reales de dimensión n e Número e 𝑥𝑗 Elemento j-ésimo del vector 𝒙 𝒙(𝑖) Columna i-ésima de la matriz 𝑿 (por comodidad se denotará esta matriz como 𝒙) 𝑨𝑇 Matriz 𝑨 traspuesta ∑𝑦(𝑖) 𝑖 Sumatorio para todo 𝑖 de los elementos 𝑦(𝑖) del vector 𝒚 ∇𝜃 Gradiente con respecto a 𝜽 𝜕𝐽(𝜽) 𝜕𝜃𝑖 Derivada parcial de la función 𝐽(𝜽) con respecto al i-esimo elemento del vector 𝜽 𝜎(z) Función sigmoide 𝑃(𝑦=1|𝒙) Probabilidad de que y valga 1 conociendo 𝒙 1{𝑎 𝑡𝑟𝑢𝑒 𝑠𝑡𝑎𝑡𝑒𝑚𝑒𝑛𝑡} Función indicadora que vale 1 si el “statement” es cierto y 0 en caso contrario 𝑓′(𝑧) Derivada de la función 𝑓(𝑧) ⨀ Función producto de Hadamard ∗ Función convolución xxii 23 1 INTRODUCCIÓN Debido al deseo de mejorar sus condiciones de vida, el ser humano ha descubierto la manera de sustituir el esfuerzo físico o mental que acarreaba la realización de ciertas tareas mediante la creación de aparatos. De acuerdo con lo anterior, la creación de los ordenadores ha permitido que el avance tecnológico tenga un avance sin precedentes en las últimas décadas. Esto ha desembocado en la resolución de cada vez más problemas de manera algorítmica, mediante bases de datos, programación top-down, programación orientada a objetos y demás métodos. Sin embargo, para resolver problemas tales como reconocer imágenes, analizar patrones o descubrir similitudes en objetos, es necesario realizar un cambio de paradigma y emplear soluciones no convencionales. Con el objetivo de resolver situaciones como las planteadas anteriormente, los científicos empezaron a estudiar las capacidades humanas cerebrales, haciendo de éstas las bases de la inteligencia artificial. Ésto permitió el desarrollo de redes neuronales, algoritmos genéticos y lógica difusa entre otros. De esta manera, el objetivo de las redes neuronales no es la resolución de problemas de índole compleja haciendo uso de una secuencia de pasos, sino como una evolución de un sistema computacional inspirado en el cerebro humano. 1.1 Historia y Estado del Arte La historia de las Redes Neuronales Artificiales proviene de varios campos que buscaron desarrollar la idea de encontrar un modelo matemático que asimilara la forma y el comportamiento de las neuronas humanas. A lo largo de los siglos XIX y XX se empezó a desarrollar teorías de aprendizaje, visión y acondicionamiento por parte de un equipo de físicos, psicólogos y neurofisiólogos. Sin embargo, aún no se disponía de un modelo específico para la neurona. En la década de los 40 se probó que las redes neuronales podían realizar cualquier operación aritmética o lógica. Sin embargo, no fue hasta la década de los 50 cuando se inventó el perceptrón, así como las reglas de aprendizaje. Desde ese momento se empezaron a desarrollar más algoritmos y la implementación de éstos creció, sin embargo, solo podían resolver una limitada serie de problemas. Ésto, unido a la falta de computación digital por parte de los ordenadores, resultó en el abandono de este campo durante un tiempo. El resurgir de la investigación en esta área fue gracias al algoritmo de retropropagación (Backpropagation algorithm). Desde entonces se han escrito muchos artículos y se han encontrado múltiples aplicaciones. Actualmente, la gran capacidad computacional de los ordenadores permite explorar nuevos conceptos, arquitecturas y reglas de aprendizaje. Con el fin de entender el estado del arte en el campo de las redes neuronales, se van a explicar algunas de las aplicaciones de éstas actualmente: • Reconocimiento de Caracteres: Esta idea tuvo mucho auge y se hizo muy popular hace algunos años, pues permite reconocer texto manuscrito. • Compresión de imágenes: Las redes neuronales pueden recibir y procesar grandes cantidades de información a la vez, siendo esto útil en la compresión de imágenes. • Mercado de Valores: Estas redes pueden examinar una gran cantidad de información de forma rápida y organizar todo de manera que se pueda hacer un adecuado estudio de proyección y predecir el valor de las acciones. 24 • Medicina: una de las áreas que más ha ganado la atención es en la detección de afecciones cardiopulmonares, es decir compara muchos modelos distintos para identificar similitud en patrones y síntomas de la enfermedad. Estos sistemas ayudan a los médicos con el diagnóstico por el análisis de los síntomas reportados y las resonancias magnéticas y rayos x. También se han usado para dispositivos analizadores de habla para ayudar a personas con sordera profunda, monitorización de cirugías, predicción de reacciones adversas a un medicamento, entendimiento de causas de ataques epilépticos. • Militar: Las redes neuronales juegan un papel importante en el campo de batalla, especialmente en aviones de combate y tanques que son equipados con cámaras digitales de alta resolución que funciona conectado a un computador que continuamente explora el exterior de posibles amenazas. De igual manera se pueden emplear para clasificar señales de radar o creación de armas inteligentes. • Diagnóstico de Maquinas: A nivel industrial cuando una de estas máquinas presenta fallos automáticamente las apaga cuando esto ocurre. • Análisis de Firmas: Las redes neuronales pueden ser empleadas para la comparación de firmas generadas. Esta ha sido una de las primeras aplicaciones implementada a gran escala en USA. • A nivel Administrativo: Identificaciones de candidatos para posiciones específicas, optimización de plazas y horarios en líneas de vuelo, minería de datos. Como se puede observar, la multitud de usos que estas redes tienen en casi cualquier campo pone de manifiesto la gran importancia de las mismas, además de hacer ver que las líneas de investigación son cada vez mayores y están llevadas a cabo por las empresas de mayor importancia internacional. Por último, destacar el logro más importante de este tipo de inteligencia artificial en los últimos años: la derrota del campeón mundial de Go a manos de AlphaGo, la inteligencia artificial de Google. Go es un juego de tablero que nació en China, y que tiene más de 2.500 años de antigüedad, cuyo objetivo final no es otro más que ser capaz de rodear con unas piedras un área mayor que la del oponente. En el caso de AlphaGo, los ingenieros de DeepMind enviaron al sistema más de 30 millones de movimientos realizados por los mayores expertos del juego Go. Luego, con una segunda red neuronal conocida como “aprendizaje de refuerzo”, hicieron que el sistema jugara contra sí mismo miles de veces. A cada partida, el “software” iba aprendiendo y se hacía cada vez más inteligente en tiempo real. Los creadores del programa hicieron que el sistema además analizara los resultados de las partidas contra sí mismo para predecir futuros resultados y planificar nuevas estrategias. Es decir, inventaron un sistema capaz de analizar el presente e intuir el futuro, un programa que toma decisiones de forma muy parecida a como lo hacemos los humanos. En la segunda partida, el surcoreano Lee Sedol, campeón mundial del juego Go, se levantó de la silla desesperado ante un movimiento incomprensible, pero que fue un golpe magistral. Desencajado, necesitó 15 minutos para improvisar una respuesta, la cual no sirvió de nada, pues cuatro horas más tarde perdió la partida. La máquina volvería a ganarle en otras dos ocasiones. 25 1.2 Motivación Esta investigación tiene como motivación adentrarse en el mundo de la visión por computador, la inteligencia artificial y, más concretamente, las redes neuronales. Dado que este campo es altamente complejo e inabarcable en su totalidad, se ha decidido empezar por los conceptos básicos y, una vez adquiridos, se procederá a crear de una red capaz de identificar y clasificar caracteres manuscritos, en concreto el MNIST (Modified National Institute of Standards and Technology). Esta base de datos es un análogo a “hello world” para la visión por computador. Desde su lanzamiento en 1999, esta base de datos ha servido como referencia para algoritmos de clasificación. Mi motivación personal es la de introducirme en el campo de la inteligencia artificial, el cual ha captado mi atención debido a la extensa cantidad de logros conseguidos en un periodo relativamente corto de tiempo. Además de todo esto, me motiva la idea de saber que al final del proyecto comprenderé gran parte de los fundamentos que usan la mayoría de servicios del Smartphone, redes sociales o incluso alguna que otra aplicación. Acciones que realizamos todos cada día, tras las cuales están presentes las redes neuronales sin ni siquiera percatarnos de ello. 1.3 Objetivos El principal objetivo del proyecto es realizar un software capaz de reconocer caracteres manuscritos con una precisión deseable por encima del 95%. Para alcanzar esta meta, se va a proceder a realizar un tutorial proporcionado por la universidad de Stanford, el cual implementa en MATLAB, desde los conocimientos más básicos hasta un apartado de reconocimiento de las imágenes del MNIST. Posteriormente, con el fin de comprobar que los conceptos están totalmente comprendidos y dado que la mayoría de sistemas actuales no usan este lenguaje, se implementará en Python. Además, esto permitirá comprobar las diferencias entre los dos lenguajes y ver cuál es más óptimo. Por tanto, el orden lógico de los objetivos a alcanzar es: • Implementar un algoritmo de aprendizaje y optimización simple para entender los conceptos básicos. • Realizar un algoritmo de aprendizaje capaz de clasificar dos clases. • Ampliar el concepto anterior para 𝑛 clases. • Crear la primera red multicapa. • Implementar una Red Neuronal Convolucional capaz de clasificar caracteres manuscritos. • Exportar la idea anterior a Python. • Optimizar, en la medida de lo posible el rendimiento y el tiempo de entrenamiento de la red. 1.4 Estructura del proyecto La estructura a seguir en el presente Proyecto es intrínsecamente relacionada con los objetivos descritos en la sección anterior. Este primer apartado es una introducción en el que se ve el porqué de realizar un estudio sobre redes neuronales, el estado del arte y los objetivos a alcanzar. Posteriormente, se explicarán los algoritmos de aprendizaje y optimización más sencillos, que servirán de base para el objetivo final, y se implementarán, viendo los resultados obtenidos y la posible optimización del código. Tras esto, se diseñará la primera red multicapa, la cual es de vital importancia para el reconocimiento de imágenes y se explicará uno de los algoritmos de optimización. A continuación, se implementará una red neuronal convolucional capaz de detectar caracteres manuscritos en MATLAB usando los conceptos anteriores. Tras obtener buenos resultados, se migrará el código a Python para comprobar cuál de los dos lenguajes es más óptimo para este tipo de tareas. Por último, se presentarán unas conclusiones generales y vías para el desarrollo posterior de este proyecto. 32 Escrito en código de MATLAB, queda: y_hat=(theta'*X-y); f=1/2*(sum(y_hat.^2)); g=X*y_hat'; Ejecutando de nuevo el código ex1a_linreg, pero pasándole esta vez a la función de optimización la función linear_regression_vec (el código completo está en el Apéndice I – Código 3), la cual calcula el coste y el gradiente según el código anterior, se obtiene como resultado: Figura 3-3: Valores RMS y tiempo de ejecución. Figura 3-4: Gráfica con precios reales y predichos Como se puede observar, los resultados obtenidos son prácticamente idénticos, sin embargo, se puede comprobar que el tiempo de ejecución es significativamente menor que en el caso anterior. 33 3.2 Logistic Regression 3.2.1 Teoría y conceptos Gracias a la regresión lineal se pueden predecir cantidades continuas como el precio de las casas como una función lineal que depende de los datos de entrada. Sin embargo, en algunas ocasiones resulta de interés predecir variables discretas, es decir, para predecir por ejemplo si un conjunto de píxeles representa un “0” o un “1”. Los problemas de este tipo son problemas de clasificación y la regresión logística es un algoritmo de clasificación bastante simple para aprender a tomar dichas decisiones. En la regresión lineal, se probó a predecir los valores de 𝑦(𝑖) usando una función lineal 𝒚=ℎ𝜃(𝒙)=𝜽𝑇𝒙. Este tipo de funciones no son las adecuadas para predecir valores binarios (𝑦(𝑖)∈{0,1}), por lo que se usa una hipótesis diferente para predecir la probabilidad de pertenecer a una clase frente a la contraria. Concretamente, se usará una función de la forma: 𝑃(𝒚=1|𝒙)=ℎ𝜃(𝒙)=1 1+𝑒−𝜽𝑇𝒙≡𝜎(𝜽𝑇𝒙) 𝑃(𝒚=0|𝒙)=1−𝑃(𝒚=1|𝒙)=1−ℎ𝜃(𝒙) (3-6) La función 𝜎(𝑧)≡1 1+𝑒−𝑧, también llamada función “sigmoide” o “logística”, contiene todos los valores de 𝜽𝑇𝒙 dentro del rango [0,1] para poder ser interpretados como probabilidad. El objetivo es encontrar un valor de 𝜽 tal que la probabilidad 𝑃(𝒚=1|𝒙) sea alta cuando 𝒙 pertenece a la clase “1” y pequeña cuando pertenezca a la clase “0”. La función coste que se usará con un set de entrenamiento cuyos targets son binarios es: 𝐽(𝜽)=−∑(𝑦(𝑖)log(ℎ𝜃(𝒙(𝑖)))+(1−𝑦(𝑖))log(1−ℎ𝜃(𝒙(𝑖)))) 𝑖 (3-7) En la función anterior, solo uno de los dos términos del sumatorio es distinto de cero para cada ejemplo del set de entrenamiento (dependiendo si la clase de 𝑦(𝑖) es “1” o “0”). Una vez que se tiene la función coste, es necesario aprender a clasificar los datos minimizando 𝐽(𝜽) para encontrar una 𝜽 óptima. Una vez hecho esto, se pueden clasificar nuevos datos como “1” o “0” comprobando las probabilidades de pertenecer a cada una de las clases. Si 𝑃(𝒚=1|𝒙)>𝑃(𝒚=0|𝒙), entonces se clasifica como “1”, “0” en caso contrario. Esto último es lo mismo que comprobar ℎ𝜃(𝒙)>0.5. Para minimizar la función coste se pueden utilizar las mismas herramientas que con la regresión lineal. Hay que proveer a la función optimizadora el cálculo del coste y del gradiente para cualquier valor de 𝜽. En este caso, las componentes del gradiente se pueden calcular como: 𝜕𝐽(𝜽) 𝜕𝜃𝑗=∑𝑥𝑗(𝑖) 𝑖(ℎ𝜃(𝒙(𝑖))−𝑦(𝑖)) (3-8) Escrito en forma vectorial, el gradiente se puede expresar como: ∇θ𝐽(𝜽)=∑𝒙(𝑖) 𝑖(ℎ𝜃(𝒙(𝑖))−𝑦(𝑖)) (3-9) Esta última ecuación es esencialmente la misma que la del gradiente para la regresión lineal, solo que ahora ℎ𝜃(𝒙)=𝜎(𝜽𝑇𝒙), por lo que, haciendo 𝒚=𝜎(𝜽𝑇𝒙)−𝒚 la equivalencia de la ecuación (3-5) sigue siendo válida. 34 3.2.2 Resultados experimentales A continuación, se va a implementar una regresión logística en MATLAB usando como ejemplo el código ex1b_logreg.m dentro de la carpeta ex1/. El código completo se encuentra disponible en el Apéndice II – Código 1, sin embargo, se procederá a explicar dicho código grosso modo y se insertarán aquellas partes del mismo que sean más interesantes, así como de los resultados obtenidos. En este ejercicio, se realizará un algoritmo capaz de clasificar imágenes del set de datos MNIST en las clases “0” y “1”. Obviamente, para que esto sea posible solo se pueden tomar imágenes del MNIST que correspondan con estos dos valores. Algunas de las imágenes del set de datos son: Figura 3-5: Imágenes del set de datos MNIST Lo primero que se realiza en el código es cargar dichas imágenes mediante la función ex1_load_mnist, la cual almacena las imágenes en los sets de entrenamiento y de prueba. Las imágenes tienen una intensidad de pixel dentro del rango [0,1] y se han pasado a formato vector. Por tanto, el i-ésimo elemento de la j-ésima columna de la matriz 𝒙 representa el i-ésimo pixel de la j-ésima imagen, es decir, cada columna de la matriz de entrada es una imagen del set de datos convertida a formato vectorial Tras obtener los sets de entrenamiento y prueba, se inicializa el vector 𝜽, se define el máximo de iteraciones y se llama a la función encargada de entrenar los parámetros, a la que se le pasa una función que calcula el coste y el gradiente, el vector 𝜽, las opciones y el set de entrenamiento. Una vez más, no se va a explicar en esta sección el algoritmo de optimización, solo la función encargada de calcular el coste y el gradiente. Dicha función es logistic_regression_vec y se puede encontrar completa en el Apéndice II – Código 2. En este caso 𝑦=𝜎(𝜽𝑇𝒙) y, por tanto, la función coste y el gradiente se pueden calcular con las siguientes líneas de código: y_hat=logsig(theta'*X); f=-sum(y.*log(y_hat)+(1-y).*log(1-y_hat)); g=X*(y_hat-y)'; Destacar que la función logsig es la función sigmoide en MATLAB. Los puntos delante de las operaciones son para que esas operaciones se realicen elemento a elemento evitando así, por tanto, la necesidad de utilizar bucles que ralentizarían dicho código. Los resultados obtenidos en este ejemplo son los siguientes: Figura 3-6: Resultados de Regresión Logística El resultado obtenido es del 100% debido a que es relativamente sencillo clasificar en “0” y “1”, pues ambos tienen a ser muy diferentes. En futuros ejercicios será mucho más complicado conseguir resultados tan perfectos como éste. 35 3.3 Softmax Regression 3.3.1 Teoría y conceptos La regresión softmax, también llamada regresión logística multinomial, es una generalización de la regresión logística para el caso en el que se manejen más de dos clases. En la regresión logística, se asumió que los objetivos eran binarios 𝑦(𝑖)∈{0,1} y se utilizó el algoritmo anterior para distinguir entre dos tipos de dígitos manuscritos. La regresión softmax permite que 𝑦(𝑖)∈{1,…,𝐾}, donde 𝐾 es el número de clases. En el caso que nos abarca, el set MNIST tiene 10 clases diferentes. Dado una entrada 𝒙, el objetivo sería estimar la probabilidad 𝑃(𝒚=𝑘|𝒙) para cada valor de 𝑘=1,…,𝐾, por lo que la salida será un vector de 𝐾 elementos cuya suma es igual a 1. Por tanto: ℎ𝜃(𝒙)=[𝑃(𝒚=1|𝒙;𝜽) 𝑃(𝒚=2|𝒙;𝜽) ⋮ 𝑃(𝒚=𝐾|𝒙;𝜽)]=1 ∑𝑒𝜽(𝑗)𝑇𝒙 𝐾 𝑗=1 [ 𝑒𝜽(1)𝑇𝒙 𝑒𝜽(2)𝑇𝒙 ⋮ 𝑒𝜽(𝐾)𝑇𝒙 ] (3-10) Aquí 𝜽(1),𝜽(2),…,𝜽(𝐾)∈ℜ𝑛 son los parámetros de nuestro modelo. Asimismo, el termino 1 ∑𝑒𝜽(𝑗)𝑇𝒙 𝐾 𝑗=1 normaliza la distribución de manera que la suma sea uno. Por conveniencia, se usará 𝜽 para denotar a todos los parámetros del modelo. Al implementar una regresión softmax, es útil representar 𝜽 como una matriz de 𝑛×𝐾 obtenida de la concatenación en columnas de 𝜽(1),𝜽(2),…,𝜽(𝐾), es decir: 𝜽=[| | | | 𝜽(1)𝜽(2)… 𝜽(𝐾) | | | | ] (3-11) A continuación, se describirá la función coste que se usará para la regresión softmax. En la ecuación (3-11), 1{∙} es una “función indicadora”, de manera que 1{𝑎 𝑡𝑟𝑢𝑒 𝑠𝑡𝑎𝑡𝑒𝑚𝑒𝑛𝑡}=1 y 1{𝑎 𝑓𝑎𝑙𝑠𝑒 𝑠𝑡𝑎𝑡𝑒𝑚𝑒𝑛𝑡}=0. Por ejemplo 1{2+2=4}=1. Teniendo esto en cuenta, la función coste es: 𝐽(𝜽)=−[∑∑1{𝑦(𝑖)=𝑘}log 𝑒𝜽(𝑘)𝑇𝒙(𝑖) ∑𝑒𝜽(𝑗)𝑇𝒙(𝑖) 𝐾 𝑗=1 𝐾 𝑘=1 𝑚 𝑖=1 ] (3-12) En la regresión softmax se tiene que: 𝑃(𝑦(𝑖)=𝑘|𝒙(𝑖);𝜽)=𝑒𝜽(𝑘)𝑇𝒙(𝑖) ∑𝑒𝜽(𝑗)𝑇𝒙(𝑖) 𝐾 𝑗=1 (3-13) 36 Debido a que no se puede obtener el mínimo de la función coste de forma analítica, es necesario utilizar una vez más un algoritmo de optimización iterativo. Este algoritmo necesita, a su vez el gradiente, que se puede expresar como: ∇𝜽(𝑘)𝐽(𝜽)=−∑[𝒙(𝑖)(1{𝑦(𝑖)=𝑘}−𝑃(𝑦(𝑖)=𝑘|𝒙(𝑖);𝜽))] 𝑚 𝑖=1 (3-14) En este caso, ∇𝜽(𝑘)𝐽(𝜽) es un vector tal que el j-ésimo elemento es 𝜕𝐽(𝜽) 𝜕𝜃𝑗(𝑘), es decir, la derivada parcial de 𝐽(𝜽) con respecto al j-ésimo elemento de 𝜽(𝑘). Con esto, ya se puede calcular el gradiente e introducírselo al algoritmo de optimización para minimizar 𝐽(𝜽). 3.3.2 Resultados experimentales A continuación, se va a implementar una regresión softmax en MATLAB usando como ejemplo el código ex1c_softmax.m dentro de la carpeta ex1/. El código completo se encuentra disponible en el Apéndice III – Código 1, sin embargo, se procederá a explicar dicho código grosso modo y se insertarán aquellas partes del mismo que sean más interesantes, así como de los resultados obtenidos. En este ejercicio se procederá a entrenar un clasificador capaz de manejar las 10 clases del set MNIST. El código es muy parecido al utilizado para la regresión logística, solo que esta vez se cargarán los sets de entrenamiento y prueba del MNIST completos. A los objetivos 𝑦(𝑖) se les ha añadido un 1, de forma que 𝑦(𝑖)∈{1,…,10} (este cambio permite usar 𝑦(𝑖) como índices dentro de una matriz). El código ejecuta las mismas operaciones que la regresión logística: carga los sets de entrenamiento y prueba y llama a la función minFunc, a la que se le pasa la función softmax_regression_vec (se puede ver en el Apéndice III – Código 2), que calcula el coste y el gradiente. Si expandimos la ecuación (3-12), se obtiene los siguientes términos: 𝐽(𝜽)=−[∑∑1{𝑦(𝑖)=𝑘}log 𝑒𝜽(𝑘)𝑇𝒙(𝑖) ∑𝑒𝜽(𝑗)𝑇𝒙(𝑖) 𝐾 𝑗=1 𝐾 𝑘=1 𝑚 𝑖=1 ]= =−[1{𝑦(1)=1}log 𝑒𝜽(1)𝑇𝒙(1) 𝑒𝜽(1)𝑇𝒙(1)+⋯+𝑒𝜽(𝐾)𝑇𝒙(1) +1{𝑦(1)=2}log 𝑒𝜽(2)𝑇𝒙(1) 𝑒𝜽(1)𝑇𝒙(1)+⋯+𝑒𝜽(𝐾)𝑇𝒙(1)+⋯ +1{𝑦(𝑚)=𝐾}log 𝑒𝜽(𝐾)𝑇𝒙(𝑚) 𝑒𝜽(1)𝑇𝒙(𝑚)+⋯+𝑒𝜽(𝐾)𝑇𝒙(𝑚)] (3-15) 37 Si uno se fija, solo se sumará para cada 𝑦(𝑖) aquel término k que cumpla la igualdad y el resto se descartarán. Para implementar dichos términos sin utilizar bucles, se procederá a explicar la analogía entre la expresión anterior, en concreto el cociente dentro del logaritmo, y las siguientes operaciones matemáticas: 𝜽𝑇𝒙=[𝜽(1)𝑇𝒙(1)𝜽(1)𝑇𝒙(2)… 𝜽(1)𝑇𝒙(𝑚) 𝜽(2)𝑇𝒙(1)𝜽(2)𝑇𝒙(2)… 𝜽(2)𝑇𝒙(𝑚) ⋮ ⋮ ⋱ ⋮ 𝜽(𝐾)𝑇𝒙(1)𝜽(𝐾)𝑇𝒙(2)… 𝜽(𝐾)𝑇𝒙(𝑚)] (3-16) Por tanto, la multiplicación matricial anterior permite, tras elevar 𝑒 a dicha matriz, tener en una matriz todos los términos que corresponden al dividendo del cociente en la ecuación (3-15). Además, es fácil darse cuenta que el divisor de cada término perteneciente a la i-ésima columna no es más que la suma de todos los elementos de la i-ésima columna de la ecuación (3-16) (tras elevar 𝑒 a dicho valor, obviamente). Por tanto, el cociente se puede realizar mediante el comando bsxfun, el cual permite realizar la división anterior. Una vez obtenido el logaritmo del cociente de la expresión, cuyo resultado es una matriz, solo resta realizar la suma de aquellos elementos de la matriz que cumplen 1{𝑦(𝑖)=𝑘}. El hecho de haber expresado 𝑦(𝑖)∈ {1,…,10} de dicha forma permite utilizarlo como índices mediante el comando sub2ind. Este comando obtiene los índices de una matriz de acuerdo a los argumentos, que actúan indicando la fila y columna. Para su mejor compresión, se supondrá que se tiene la matriz 𝑨=[1 3 4 10 9 8 2 7 5] y que se usa el comando ind=sub2ind(size(A), [1 3 2], 1:3);. Con dicho comando las filas son [1 3 2] y las columnas [1 2 3], por lo que se obtendrán los índices de los elementos 𝐴11,𝐴32 y 𝐴23. Si se ejecuta el comando A(ind), se obtiene un vector con valores [1 7 8]. Por tanto, con este comando se pueden obtener aquellos índices que cumplen 1{𝑦(𝑖)=𝑘}, haciendo de 𝒚 las filas y de 1:size(A) las columnas. A continuación, se va a mostrar el código en MATLAB: y_hat=exp(theta'*X); % It’s the numerator y_hat=[y_hat; ones(1,size(y_hat,2))]; % Last row set to 1 to have 10 classes y_hat_sum=sum(y_hat,1); % It’s the denominator y_hat_sum=y_hat_sum-1; % As last row doesn’t count, substract it from the total sum coc=bsxfun(@rdivide,y_hat,y_hat_sum); A=log(coc); ind=sub2ind(size(A),y,1:size(y_hat,2)); A(end,:)=0; % Last row set to 0 because the prob that y belongs to 10 is null f=-sum(A(ind)); B=zeros(size(A)); B(ind)=1; g=-X*(B-coc)'; g=g(:,1:end-1); % Last row doesn’t count g=g(:); % make gradient a vector for minFunc El hecho de establecer 𝑦(𝑖)∈{1,…,10} tiene la ventaja de poder usarse como índices, pero la desventaja es que ahora 𝑃(𝑦(𝑖)=10|𝒙(𝑖);𝜽)=0. Por tanto, se trabajará como si fueran 10 clases (de ahí el incluir una fila más en y_hat) y posteriormente se eliminará dicha clase del coste y del gradiente, tal y como se observa en el código anterior. 38 A continuación, se va a proceder a mostrar los resultados obtenidos al ejecutar el script: Figura 3-7: Resultados Regresión Softmax Como se puede observar en la imagen anterior, el resultado es superior al 87%, lo cual no está nada mal teniendo en cuenta que este es el algoritmo más simple que se puede crear con dicho propósito. Además, es notable el hecho de que el tiempo de simulación ya no es tan pequeño como en los casos anteriores. Si se hubieran realizado las operaciones usando bucles, el tiempo se habría incrementado de manera exponencial. Una vez entendidos los algoritmos más básicos, es necesario aumentar el nivel y adentrarse en las redes neuronales. 39 4 RED NEURONAL MULTICAPA 4.1 La biología como Fuente Para comprender las redes neuronales artificiales, primero se debe conocer cómo funcionan las naturales y saber las partes de las que se componen. Una neurona está compuesta por las dendritas, ramificaciones que recogen los impulsos eléctricos y los conducen al cuerpo de la célula, también llamado soma; el cuerpo o soma, que procesa la información recibida y el axón, que conduce la señal proveniente del cuerpo a otras neuronas. La conexión entre neuronas se realiza mediante la sinapsis. Ésta puede ser más o menos fuerte dependiendo del proceso químico encargado. A continuación, se muestra una imagen de una neurona para facilitar su comprensión: Figura 4-1: Esquema de una neurona Las Redes Neuronales Artificiales no se aproximan a la complejidad del cerebro, pero tienen los elementos clave que las asimilan a las biológicas para resolver problemas concretos Aunque los circuitos electrónicos son mucho más rápidos que las neuronas biológicas, la capacidad de estas últimas para trabajar en paralelo hace que sean capaces de resolver cualquier tarea de una manera más rápida que un ordenador. Por ello, las Redes Neuronales Artificiales son ideales para su implementación en VLSI, dispositivos ópticos y procesadores en paralelo. 4.2 La neurona artificial Si se considera un problema de aprendizaje supervisado en el que se tiene acceso a determinados ejemplos para el entrenamiento (𝒙(𝑖),𝑦(𝑖)), las redes neuronales permiten definir la forma de la función ℎ𝑊,𝑏(𝒙), la cual puede ser compleja o incluso no lineal. 40 Antes de describir redes neuronales, se empezará describiendo la red neuronal más sencilla posible, es decir, aquella que se corresponde con una sola “neurona”. El diagrama de una neurona es el siguiente: Figura 4-2: Diagrama de una neurona artificial Esta neurona es una unidad computacional que toma como entradas 𝑥1,𝑥2,𝑥3,𝑥4,…,𝑥𝑛 y produce una salida de acuerdo con la ecuación 𝑦=ℎ𝑊,𝑏(𝒙)=𝑓(𝑾𝑇𝒙+𝑏)=𝑓(∑ 𝑊𝑖𝑥𝑖+𝑏 𝑛𝑖=1 ), donde 𝑓 (𝜑 en la figura 4-2) es la función de activación y su dominio e imagen son 𝑓:ℜ→ℜ. Se tomará 𝑓(∙) como la función sigmoide: 𝑓(𝑧)=1 1+𝑒−𝑧 (4-1) Por tanto, nuestra neurona simple se corresponde exactamente con el algoritmo de regresión logística visto en el apartado 3.2. Aunque se va a utilizar la función sigmoide, cabe destacar que otra elección bastante común de 𝑓 es la tangente hiberbólica o tanh: 𝑓(𝑧)=tanh (𝑧)=𝑒𝑧−𝑒−𝑧 𝑒𝑧+𝑒−𝑧 (4-2) Además, investigaciones recientes han encontrado una función de activación diferente: la función lineal rectificada (rectified linear function o ReLU en inglés). Esta función trabaja bien en redes neuronales con Deep Learning y es diferente de la sigmoide y la tanh pues no está limitada y no es continuamente diferenciable. Esta función viene dada por: 𝑓(𝑧)=max (0,𝑧) (4-3) 41 A continuación, se muestran en una sola gráfica las tres funciones mencionadas anteriormente para su mejor visualización y comparación: Figura 4-3: Comparación de funciones de activación La función tanh (𝑧) es una versión reescalada de la sigmoide, pues el rango de salida es de [−1,1] en vez de ser de [0,1]. La función lineal rectificada es una función lineal definida a trozos que satura el valor a 0 cuando la entrada 𝑧 es menor que dicho valor. Por último, comentar una igualdad que será de utilidad en el futuro: si 𝑓(𝑧)=1/(1+𝑒𝑧) es la función sigmoide, entonces su derivada viene dada por 𝑓′(𝑧)=𝑓(𝑧)(1−𝑓(𝑧)). Si, por el contrario, 𝑓(𝑧) es la función tanh, entonces su derivada resulta ser 𝑓′(𝑧)=1−(𝑓(𝑧))2. La función lineal rectificada tiene un gradiente de 0 cuando 𝑧≤0 y 1 en cualquier otro caso. 4.3 Modelo de Red Neuronal Una red neuronal se crea interconectando las neuronas simples que se vió en el apartado anterior de manera que la salida de una neurona puede ser la entrada de otra. Por ejemplo, he aquí una pequeña red neuronal: Figura 4-4: Red neuronal con varias neuronas 48 Tras esto, realiza el código de propagación hacia adelante o forward propagation, almacenando en cada una de las activaciones en una variable tipo celda. Para esto, se ha seguido la ecuación (4-6) y escrito en código queda como: for l = 1 : numHidden if l > 1 hAct{l} = stack{l}.W * hAct{l - 1} + repmat(stack{l}.b, 1,numSamples); else hAct{l} = stack{l}.W * data + repmat(stack{l}.b, 1, numSamples); end hAct{l} = sigmoid(hAct{l}); end A continuación, se calcula la activación de la capa de salida y se calcula la ecuación (4-18) de la misma manera que se hizo en el apartado 3.3.2 teniendo en cuenta que ahora no existe una clase 0 (porque se ha aumentado la etiqueta de la clase en una unidad), por lo que no es necesario realizar las operaciones para obtener una probabilidad nula en la clase 0. Antes de calcular el coste, calcula las probabilidades de pertenecer a cada clase y obtiene el máximo para conocer la clase final predicha. Si el flag po está activo, se vuelve a la función principal (esto sirve en el caso de que la red esté entrenada y se quiera conocer la salida producida). l = numHidden+1; y_hat = stack{l}.W * hAct{l - 1} + repmat(stack{l}.b, 1, numSamples); y_hat = exp(y_hat); hAct{l} = bsxfun(@rdivide, y_hat, sum(y_hat, 1)); [pred_prob pred_labels] = max(hAct{l}); pred_prob = hAct{l}; %% return here if only predictions desired. if po cost = -1; ceCost = -1; wCost = -1; numCorrect = -1; grad = []; return; end; %% compute cost y_hat = log(hAct{numHidden+1}); index = sub2ind(size(y_hat), labels', 1:numSamples); ceCost = -sum(y_hat(index))/(size(data,2)); Una vez calculado el coste, solo resta calcular el gradiente mediante el método de la retropropagación descrito en el apartado anterior. Las ecuaciones (4-14), (4-15) y (4-16) implementadas en MATLAB quedan como for l = numHidden+1 : -1 : 1 if l > numHidden gradFunc = ones(size(gradInput)); else gradFunc = hAct{l} .* (1 - hAct{l}); end gradOutput = gradInput .* gradFunc; if l > 1 gradStack{l}.W = gradOutput * hAct{l-1}'./(size(data,2)); else gradStack{l}.W = gradOutput * data'./(size(data,2)); end gradStack{l}.b = sum(gradOutput, 2)/(size(data,2)); gradInput = stack{l}.W' * gradOutput; end 49 El último paso es añadir el término que depende de 𝜆 a la ecuación del coste y al gradiente de los pesos en cada una de las capas. Tras esto, se pasa el gradiente a forma vectorial y la función minFunc es la que se encarga de actualizar los pesos gracias al coste y el gradiente proporcionados. %% compute weight penalty cost and gradient for non-bias terms wCost = 0; for l = 1:numHidden+1 wCost = wCost + .5 * ei.lambda * sum(stack{l}.W(:) .^ 2); end cost = ceCost + wCost; % Computing the gradient of the weight decay. for l = numHidden : -1 : 1 gradStack{l}.W = gradStack{l}.W + ei.lambda * stack{l}.W; end Una vez explicado el código, es el momento de la ejecución para poder comprobar los resultados obtenidos mediante la implementación de esta red neuronal: Figura 4-6: Resultados obtenidos con Red Neuronal Multicapa Como se puede comprobar el resultado actual está muy por encima del obtenido mediante el uso de la regresión softmax de forma aislada. La precisión obtenida al probar el set de entrenamiento es perfecta y la del de prueba es del 97.08%, la cual es bastante alta del mismo modo. El único punto negativo es que ahora el código tarda bastante más en ejecutarse en comparación, sin embargo, los resultados obtenidos hacen que valga la pena el uso de esta red y el consecuente tiempo de espera. 50 51 5 REDES NEURONALES CONVOLUCIONALES 5.1 Extracción de características usando la convolución Hasta ahora, a lo largo del presente proyecto se ha trabajado con imágenes contenidas en el set de datos MNIST, las cuales tienen una resolución bastante baja. Esto no es lo común a la hora de usar redes neuronales para clasificar imágenes en una clase, ya que la resolución actual de las imágenes asciende como mínimo al HD (1280×720 píxeles). Por tanto, es necesario realizar un cambio de paradigma para desarrollar métodos que permitan el uso de set de datos más realistas, cuyas imágenes sean de resoluciones superiores. Además de lo anterior, existe un problema añadido: en el diseño anterior se decidió conectar todas y cada una de las unidades de la capa oculta a todas las unidades de entrada. Al trabajar con imágenes cuya resolución es relativamente baja, resulta computacionalmente factible aprender características usando la imagen completa, sin embargo, aprender características que abarquen la imagen en su totalidad usando resoluciones mayores (96×96, por ejemplo) es muy costoso hablando en términos de cálculo computacional. Con la resolución de ejemplo anterior, se tendrían alrededor de 104 unidades de entrada y, asumiendo que se quisieran aprender 100 características, el número de parámetros oscila en 106. Aparte de esto, la propagación, tanto hacia adelante como hacia atrás, resultaría ser 102 veces más lenta comparada con imágenes de 28×28. Una solución bastante simple al problema descrito es restringir las conexiones entre las unidades de la capa oculta y las de entrada, de manera que cada unidad de entrada queda conectada solamente a un pequeño subconjunto de las unidades de entrada. Concretamente, cada unidad de la capa oculta se conectará a una región continua de pixeles de la entrada. La idea de tener redes locales también se inspira en cómo funciona el sistema visual biológico, ya que las neuronas del córtex visual poseen campos receptivos que responden solo ante estímulos localizados en una región específica. Las imágenes tienen la propiedad natural de ser “estacionarias”, es decir, que los datos estadísticos de una parte de la imagen son los mismos que en cualquier otra parte. Esto sugiere que las características que se aprendan en una parte de la imagen pueden ser aplicadas a otras partes de dicha imagen o incluso en todas las partes de ella. De forma más específica, tras haber aprendido características sobre pequeñas áreas (de 8×8, por ejemplo) llamadas filtros o kernels, se puede aplicar este detector de características en cualquier lugar de la imagen. Por tanto, se puede tomar este kernel y realizar la convolución con la imagen para obtener diferentes valores de activación en cada región de la imagen. A continuación, se ofrece una imagen visual para comprender mejor la convolución: Figura 5-1: Convolución 52 Matemáticamente, dada una imagen grande 𝒙𝒍 de dimensiones 𝑟 × 𝑐 y 𝑘 filtros 𝒙𝒔 de dimensiones 𝑎 × 𝑏, entonces se generarán 𝑘 mapas de rasgos que en total forman un conjunto de (𝑟−𝑎+1) × (𝑐−𝑏+1) × 𝑘. Si se expresa cada uno de lo filtros en forma de 𝑾 y 𝒃, para cada subconjunto 𝒙𝒔𝒂𝒎𝒑𝒍𝒆 de 𝑎 × 𝑏 de la imagen 𝒙𝒍, se calcula 𝑓𝑠𝑎𝑚𝑝𝑙𝑒=𝜎(𝑾(1)𝒙𝒔𝒂𝒎𝒑𝒍𝒆+𝒃(1)). Para implementar la convolución en MATLAB es necesario calcular 𝜎(𝑾𝒙(𝒓,𝒄)+𝒃) para todos los (𝒓,𝒄) válidos. Válido (valid) significa que el kernel está completamente contenido en la imagen; lo cual es lo contrario a completo (full), donde se puede extender el kernel al exterior de la imagen asumiendo estos valores como nulos. Figura 5-2: Convolución completa En la ecuación anterior 𝑾 y 𝒃 son los pesos y las bias aprendidos que conectan la capa de entrada con la capa oculta y 𝒙(𝒓,𝒄) es el subconjunto de 𝑎 × 𝑏 que empieza en la esquina superior izquierda de la imagen. Para realizar dicha ecuación, un primer método sería realizar un bucle que recorra todos los subconjuntos y calcule 𝜎(𝑾𝒙(𝒓,𝒄)+𝒃), sin embargo, resulta que en la práctica este método es válido pero muy lento. Un segundo método es calcular 𝑾𝒙(𝒓,𝒄) para todo (𝒓,𝒄), sumar 𝒃 a los valores calculados y aplicar la función sigmoide al resultado obtenido. Este método no parece significativamente mejor que el anterior pues es necesario un bucle para calcular el primer término, sin embargo, dicho bucle se puede reemplazar por una función propia de MATLAB optimizada para la convolución: conv2, acelerando bastante el proceso. A la hora de usar la función anterior es necesario tener tos aspectos muy en cuenta. El primero de ellos es que la convolución a realizar es en 2𝐷, pero en el código se manejarán 4 dimensiones: número de imágenes, número de filtros, filas de la imagen y columnas de la imagen. Por ello, es necesario realizar la convolución de forma independiente para cada filtro aplicado a cada una de las imágenes, usando las filas y columnas de la imagen pertinente como las dos dimensiones sobre las que se va a convolucionar. Esto significa que hay que usar dos bucles: uno para el número de imágenes y otro para el número de filtros. Dentro de los dos bucles se realizará la convolución entre la i-ésima imagen y el j-ésimo filtro. El segundo de los aspectos es que, por la definición matemática de la convolución, el filtro es rotado 180º antes de ejecutar ejecutar la convolución. Esto lo hace MATLAB intrínsecamente en la función, por lo que para usar el filtro deseado es necesario rotarlo antes de pasárselo a la función conv2. (123 456 789)→(9 8 7 6 5 4 3 2 1) (5-1) 53 A continuación, se muestra el código necesario para hacer la convolución de los 𝑘 filtros para cada una de las imágenes: for imageNum = 1:numImages for filterNum = 1:numFilters % convolution of image with feature matrix convolvedImage = zeros(convDim, convDim); % Obtain the feature (filterDim x filterDim) needed during the convolution filter=W(:,:,filterNum); % Flip the feature matrix because of the definition of convolution, as explained later filter = rot90(filter,2); % Obtain the image im = images(:, :, imageNum); % Convolve "filter" with "im", adding the result to convolvedImage convolvedImage=conv2(im,filter,'valid'); % Add the bias unit % Then, apply the sigmoid function to get the hidden activation convolvedImage=convolvedImage+b(filterNum); convolvedImage=logsig(convolvedImage); convolvedFeatures(:, :, filterNum, imageNum) = convolvedImage; end end 5.2 Pooling Tras extraer las características usando la convolución, el siguiente paso sería utilizarlas para la clasificación. En teoría, se podrían usar todas las características extraídas con un clasificador tipo softmax, pero esto resultaría todo un reto computacionalmente hablando. Si se considera, por ejemplo, imágenes de 96×96 (no es una alta resolución) y 400 filtros de 8×8, a la salida de la capa convolucional se tendría un vector (tras convertir de matriz a vector) de dimensión (96−8+1)2∙400=3168400. Aprender a clasificar con entradas de mas de 3 millones de elementos puede ser inviable e incluso propenso al sobreajuste u over-fitting. Para manejar esta situación, primero se decidió obtener características mediante convolución pues las imágenes tienen propiedades “estacionarias” que permiten usar características útiles en una parte de la imagen en cualquier otra parte. Por tanto, para describir una imagen grande, resulta natural analizar en que zonas de la imagen éstas características están más presentes. Por ejemplo, se puede calcular la media (o el máximo) valor de una característica a lo largo de una región de la imagen. La salida que se obtiene es mucho menor en cuanto a dimensión se refiere (comparada con los de la salida de la capa de convolución) y pueden mejorar los resultados (menos sobreajuste). La operación que hace esto posible se llama “pooling” e incluso a veces “mean pooling” o “max pooling” (dependiendo de la operación usada). La siguiente imagen muestra cómo se realiza la operación de pooling a lo largo de una imagen cualquiera: Figura 5-3: Pooling 54 De una manera más formal, tras obtener las características resultantes de la convolución tal y como se describió en el apartado 5.1, se decide realizar la operación de pooling sobre regiones de dimensiones 𝑚 × 𝑛. Por tanto, se divide la salida de la capa de convolución en regiones disjuntas de 𝑚 × 𝑛 y se toma la media o el máximo, permitiendo así obtener unas nuevas características que ya sí se pueden usar para la clasificación. Para realizar la operación pooling en MATLAB se puede utilizar de nuevo la función conv2, cuyos argumentos ahora serían la salida de la capa convolucional y un filtro que sea una matriz de unos. Una vez realizado este paso, se toman de la matriz los elementos correspondientes y se les hace la media. Cabe destacar que la convolución anterior debe ser valid. En código quedaría: W=ones(poolDim,poolDim); for imageNum=1:numImages for filterNum=1:numFilters filter=W; pooledImage=zeros(convolvedDim / poolDim, convolvedDim / poolDim); im=convolvedFeatures(:,:,filterNum,imageNum); convolvedImage=conv2(im,filter,'valid'); pooledImage=convolvedImage(1:poolDim:end, 1:poolDim:end)./(poolDim^2); pooledFeatures(:,:,filterNum,imageNum)=pooledImage; end end 5.3 Optimización: Stochastic Gradient Descent Algunos métodos de optimización como el LBFGS (el que se usó en los apartados anteriores) usan el set de entrenamiento completo para actualizar los parámetros en cada iteración y tiende a converger de manera eficiente hacia un mínimo absoluto. Sin embargo, en la práctica el cálculo del coste y del gradiente para todo el set de entrenamiento puede resultar muy lento e incluso imposible para una sola máquina si el set de datos es demasiado grande para poder guardarse en memoria. Además, otro problema con estos métodos de optimización es que no hay una forma fácil de incorporar nuevos datos “en línea”. El método del Stochastic Gradient Descent (Descenso por Gradiente Estocástico) o SGD resuelve estos dos problemas siguiendo la dirección negativa del gradiente de la función objetivo tras haber visto solo uno o unos pocos ejemplos de entrenamiento. El alto coste de retropropagar el set completo de entrenamiento hace necesario el uso del SGD. El algoritmo estándar del gradient descent actualiza los parámetros 𝜽 de forma: 𝜽=𝜽−𝛼∇𝜃𝐸[𝐽(𝜽)] (5-2) En la ecuación anterior la esperanza se aproxima por el gradiente de la función coste usando el set ompleto de entrenamiento. El método del SGD simplemente elimina dicha esperanza y calcula el gradiente usando solo uno o unos pocos ejemplos de entrenamiento. La nueva actualización viene ahora dada por. 𝜽=𝜽−𝛼∇𝜃𝐽(𝜽;𝒙(𝑖),𝑦(𝑖)) (5-3) El par (𝒙(𝑖),𝑦(𝑖)) es un par del set de entrenamiento. 55 Generalmente, cada actualización en el SGD se calcula con vista a unos pocos ejemplos del set de entrenamiento (llamados minibatch) en lugar de solo a un ejemplo. La razón tras esto es doble: en primer lugar, se reduce la varianza en la actualización de los parámetros y esto permite una convergencia más estable, y, en segundo lugar, permite el uso de álgebra lineal para el eficiente cálculo del coste y del gradiente. El tamaño típico de un minibatch es de 256, aunque el tamaño óptimo puede variar según la aplicación y la arquitectura. En el SGD el ratio de aprendizaje 𝛼 normalmente es mucho más pequeño que el usado en el método del gradient descent pues la varianza en la actualización es menor. Escoger el valor adecuado para el ratio de aprendizaje y sus actualizaciones a lo largo del entrenamiento puede ser realmente complicado. Un método estándar que funciona bien en la práctica es usar un ratio pequeño y constante que proporcione una convergencia estable en la epoch (uso completo del set de entrenamiento mediante minibatches) inicial y reducir a la mitad su valor a medida que la convergencia se ralentiza. Otro método bastante utilizado es actualizar el ratio de aprendizaje en cada iteración 𝑡 como 𝑎/(𝑏+𝑡), donde 𝑎 es el ratio inicial y 𝑏 la iteración en la que comienta dicha actualización. Los métodos más sofisticados incluyen el uso de backtracking (vuelta atrás) a lo largo de una línea para encontrar la actualización óptima. Por último, pero no menos importante, en el SGD importa el orden en el que los datos son presentados al algoritmo. Si los datos son presentados en un orden cualquiera, ésto puede introducir un offset en el gradiente y conducir a una convergencia más pobre. Un buen método para evitar esto es, por lo general, “barajar” de forma aleatoria los datos antes de cada epoch de entrenamiento. Si la función coste tiene la forma análoga a un largo desfiladero poco profundo que conduce al mínimo y con paredes empinadas a los lados, el SGD tenderá a oscilar a lo ancho del estrecho desfiladero, pues la dirección negativa al gradiente apuntará a alguna de las paredes en lugar de al mínimo. Esta forma, que puede parecer muy poco probable, es bastante común en arquitecturas más complejas y, por tanto, la convergencia puede llegar a ser muy lenta. El método Momentum es más óptimo en este caso pues apunta al mínimo más rápido a lo largo del desfiladero. La actualización de los parámetros con este método es: 𝒗=𝛾𝒗−𝛼∇𝜃𝐽(𝜽;𝒙(𝑖),𝑦(𝑖)) 𝜽=𝜽−𝒗 (5-4) En la ecuación anterior 𝒗 es el vector velocidad actual, cuyas dimensiones coinciden con las de 𝜽. El ratio de aprendizaje 𝛼 se toma como se describió previamente, aunque, al usar este método, 𝛼 tal vez deba ser más pequeño dado que la magnitud del gradiente será mayor. Por último, 𝛾∈(0,1] determina por cuantas iteraciones se incorporarán los gradientes anteriores. Generalmente, 𝛾 se establece a 0.5 hasta que el aprendizaje se estabiliza, momento en el que se aumenta a 0.9. 5.4 Red Neuronal Convolucional Una Red Neuronal Convolucional o Convolutional Neural Network (CNN) es una red compuesta por una o más capas convolucionales (seguidas normalmente de una capa de pooling) y con una o varias capas fully connected (totalmente conectadas, como las del apartado 4) tras éstas. La arquitectura de una CNN está diseñada para aprovechar la entrada como una sola estructura 2D. Esto se consigue con conexiones locales y pesos asociados a ellas, seguidos por alguna operación de pooling cuyas salidas son características invariantes a las traslaciones. Otra ventaja de estas redes es que son más fáciles de entrenar y pueden tener incluso menos parámetros que las redes del apartado 4 con el mismo número de redes ocultas. La entrada de una capa convolucional es una imagen de dimensión 𝑚 × 𝑚 × 𝑟, donde 𝑚 es el alto y el ancho de la imagen y 𝑟 es el número de canales, el cual valdrá 3 en una imagen RGB. La capa convolucional contará con 𝑘 filtros (o kernels) de dimensión 𝑛 × 𝑛 × 𝑞, donde 𝑛 es más pequeña que la dimensión de la imagen y 𝑞 puede ser del mismo valor que el número de canales 𝑟 o menor y puede variar según el kernel. A la salida de esta capa se tendrán 𝑘 mapas de características (feature maps) de dimensión 𝑚−𝑛+1 tal y como se vió en el apartado 5.1. Posteriormente, se submuestrea dicha salida mediante alguna operación de pooling a lo largo de regiones continuas de dimensión 𝑝 × 𝑝, donde 𝑝 vale normalmente 2 para imágenes pequeñas como las del 56 MNIST y no más de 5 por lo general para imágenes grandes. Además, antes o después de la capa de pooling se añade un término de bias y también una función de activación como la sigmoide a cada uno de los mapas de características. Tras estas capas, puede haber cualquier número de capas fully connected, las cuales son idénticas a las vistas en el apartado 4. A continuación, se muestra la salida tras aplicar a una imagen una capa convolucional y otra de pooling. Las unidades del mismo color tienen los pesos asociados: Figura 5-4: Capa convolucional y de pooling aplicadas a una imagen Se va a ilustrar, también, una de las configuraciones de CNN más comunes: Figura 5-5: Red Neuronal Convolucional 57 Una vez se ha entendido como funciona una CNN y el método que se va a emplear para entrenarla, solo resta calcular los gradientes de la ecuación (5-4), para lo que se volverá a usar el algoritmo de retropropagación, el cual permite calcular dichos gradientes de una forma rápida y bastante simple. Si se supone que 𝜹(𝑙+1) es el término de error de la capa (𝑙+1) de la red con una función coste 𝐽(𝑾,𝒃;𝒙,𝑦), donde (𝑾,𝒃) son los parámetros y (𝒙,𝑦) es el par del set de entrenamiento, y si, además, la capa 𝑙 esta completamente conectada a la (𝑙+1), entonces el error de la capa 𝑙 se calcula como: 𝜹(𝑙)=((𝑾(𝑙))𝑇𝜹(𝑙+1))⨀𝒇′(𝒛(𝑙)) (5-5) Y los gradientes quedan: ∇𝑾(𝑙)𝐽(𝑾,𝒃;𝒙,𝑦)=𝜹(𝑙+1)(𝒂(𝑙))𝑇 ∇𝒃(𝑙)𝐽(𝑾,𝒃;𝒙,𝑦)=𝜹(𝑙+1) (5-6) Sin embargo, si la capa 𝑙 es convolucional o de pooling, entonces el error se propaga según la siguiente ecuación: 𝜹𝑘(𝑙)=𝑢𝑝𝑠𝑎𝑚𝑝𝑙𝑒((𝑾𝑘(𝑙))𝑇𝜹𝑘(𝑙+1))⨀𝒇′(𝒛𝑘(𝑙)) (5-7) En la ecuación anterior 𝑘 es un índice que hace referencia al filtro correspondiente y 𝒇′(𝒛𝑘(𝑙)) es la derivada de la función de activación. La operación upsample debe propagar el error a través de la capa de pooling calculando el error para cada unidad de dicha capa, es decir, si se usa la operación mean pooling, entonces la operación upsample simplemente distribuye el error de forma uniforme entre aquellas unidades que sirvieron como entrada. Si, por el contrario, se usó max pooling, entonces aquella unidad de entrada que tenía el valor máximo es la que recibe todo el error. A continuación, se muestra un ejemplo práctico para comprender este concepto mejor: { 𝑚𝑒𝑎𝑛 𝑝𝑜𝑜𝑙𝑖𝑛𝑔:𝜹=(1 2 3 4) →(1 1 2 2 1 1 2 2 3 3 4 4 3 3 4 4) 𝑚𝑎𝑥 𝑝𝑜𝑜𝑙𝑖𝑛𝑔:𝜹=(1 2 3 4) →(1 0 0 2 0 0 0 0 0 0 0 0 3 0 0 4) (5-8) 64 6.2 Diferencias con respecto a MATLAB A pesar de que ambos códigos van a realizar las mismas funciones, existen grandes diferencias a la hora de implementarlo en Python con respecto a MATLAB. Estas diferencias son imperceptibles a menos que uno realice una depuración completa del código comprobando que los resultados en ambos deben ser idénticos. La primera diferencia es como Python carga los datos. Más allá de que el código sea bastante más complicado, la diferencia radica en los datos una vez cargados, ya que, por ejemplo, MATLAB carga los archivos binarios como imágenes cuya intensidad de píxeles están comprendida entre [0,1], mientras que en Python se cargan en el rango [0,255]. La solución a esto es bastante simple, pues solo hay que dividir la imagen entre 255. Figura 6-1: Carga del set MNIST Aparte de esto, existe una gran diferencia entre ambos lenguajes de programación, no solo al manejar este set de datos, sino al manejar matrices en general. Si uno se fija en la línea 52 de la figura anterior, se pueden reseñar tres cosas. La primera de ellas es que se pueden concatenar determinados comandos usando un punto, es decir, concatena las órdenes fromfile (toma los datos del archivo) y reshape (cambia las dimensiones de la matriz). El comando np. no cuenta como concatenación pues simplemente indica que se está usando la librería numpy. Esto evita el uso de más líneas y permite compactar el código. De no ser por esta concatenación, los comandos quedarían: img=np.fromfile(fimg, dtype=uint8) img=np.reshape(-1,rows,cols) Los dos restantes puntos reseñables residen en el comando reshape. En primer lugar, parece obvio, pero es necesario que las dimensiones de destino sean compatibles con las de origen. Para ilustrar esto, si se tiene una matriz de dimensiones 4×4×8=128, una posibilidad para redimensionarla sería 16×8=128, pero no se podría redimensionar a 15×4=60, ya que el número de elementos no es el mismo en el origen y destino. Esto también ocurre en MATLAB, sin embargo, es necesario explicar el significado del -1 en el comando reshape. Tomando como ejemplo las matrices anteriores, el comando reshape(-1,8) proporcionará una matriz de 8 columnas y el número de filas necesarias para que el número de elementos entre origen y destino permanezca constante, es decir, 16. Por último, pero no menos importante, existe una gran diferencia a la hora de manejar matrices de más de dos dimensiones. Mientras que en MATLAB las dimensiones de una imagen se representan mediante (rows,columns,channels), en Python las dimensiones superiores a dos se apilan hacia la izquierda, de modo que queda (channels,rows,columns). Esto hay que tenerlo muy en cuenta, ya que a lo largo del código se usan operaciones con matrices de hasta 4 dimensiones (esta última en Python quedaría a la izquierda de channels). 65 Otra diferencia importante entre MATLAB y Python es que en MATLAB el índice inicial para los vectores y cada dimensión de las matrices era 1, mientras que en Python es 0. Si esto no se tiene en cuenta, se pueden ocasionar errores bastante graves, pues se podría elegir una imagen o filtro diferente a la que debería. Además, aparte de tener esto en cuenta a la hora de indexar matrices, también es importante recalcar que, dado que en MATLAB no se podía indexar con 0, aquellas labels que eran 0 debían establecerse a 10 (para poder usarse como índices en sub2ind, explicado en capítulo 3.3.2). En Python no es necesario realizar dicha operación pues el índice inicial es 0. Con respecto a la función sub2ind, ésta no existe en Python y no hay ninguna que se le asemeje para poder extraer las probabilidades correctas de la matriz A, tal como se hizo en los apartados 3.3.2, 4.5 y 5.5. La única solución encontrada es la realización de un bucle, guardando en un vector las probabilidades para después sumarlas, tal como indica la ecuación (4-18): Figura 6-2: Función sub2ind en Python Una diferencia mínima, pero que es importante recalcar es la operación “elevar un número a otro”. Frente a la mayoría de lenguajes de programación en los que el comando sería y=x^2, en Python el comando resulta ser y=x**2. Es una diferencia casi sin importancia, pero que cualquier persona nueva en este lenguaje suele ignorar y que puede causar errores en la programación si no se tiene cuidado. La última diferencia reseñable con respecto al código en MATLAB reside en la función reshape. Existen dos formar de redimensionar una matriz: ordenando los índices como C o como Fortran (F). El primero de ellos toma las filas sucesivamente de la matriz de origen y las va apilando en las filas de la matriz destino, mientras que el método F toma las columnas de la matriz origen y las apila como columnas en el destino. Para ilustrar esto: { 𝐶:𝑨=(1 2 3 4 5 6 7 8 910 11 12)3×4→𝑩=(123 456 789 10 11 12)4×3 𝐹:𝑨=(1 2 3 4 5 6 7 8 910 11 12)3×4→𝑩=(1611 510 4 9 3 8 2712)4×3 (6-1) En MATLAB el método por defecto es Fortran, mientras que en Python es C. En Python, para matrices 2D, es sencillo cambiar el método de indexado a Fortran. Por ejemplo, realizar la operación anterior en C resultaría en el comando B=np.reshape(A, (4,3)) y en F B=np.reshape(A, (4,3), order=’F’). El parámetro order permite cambiar entre ambos métodos (no es necesario usar order=’C’ en el primer comando porque viene por defecto). Sin embargo, si la matriz de origen tiene más de dos dimensiones, ambos métodos se complican, ya que en Python, a diferencia de MATLAB, la estructura es (channels,rows,columns). Esto se traduce en que al utilizar la orden reshape, independientemente del método de indexación, primero tomará elementos de la dimensión channels, y esto no es nada deseable porque se estarían mezclando los valores de las distintas imágenes y filtros, desorganizando todos los valores importantes para la red. 66 Si uno supone que tiene en Python una matriz 𝑨 de dimensiones 2×2×2, cuyos channels son: 𝑨{𝐶ℎ𝑎𝑛𝑛𝑒𝑙 1:(1 2 3 4) 𝐶ℎ𝑎𝑛𝑛𝑒𝑙 2:(5 6 7 8) (6-2) En la siguiente tabla se va representar lo que se obtiene en Python con los dos métodos de indexación, así como el resultado deseado al redimensionalizar la matriz como 4×2: MÉTODO DE INDEXACIÓN MATRIZ OBTENIDA MATRIZ DESEADA INDEXACIÓN POR FILAS (C) (1 2 3 4 5 6 7 8) (1 2 3 5 5 6 7 8) INDEXACIÓN POR COLUMNAS (F) (1 2 5 6 3 4 7 8) (1 5 3 7 2 6 4 8) Tabla 6-1: Comparación de matrices según el método de indexación Debido a lo anterior, para redimensionar una matriz multidimensional en Python indexando mediante Fortran, no se puede realizar el comando reshape, pues el resultado no es el deseado. En el siguiente apartado, se abordará una solución para este problema. Si se realiza el código utilizando la indexación C, el resultado de la red está bastante lejos de ser bueno, por lo que es necesario implementar el método Fortran. 6.3 Soluciones y Optimización En el apartado anterior se comentaron las diferencias más significativas que se encontraron a la hora de programar en Python frente a MATLAB. En el presente apartado se va a explicar cómo se llegó a conocer dichas diferencias, las soluciones adoptadas y las optimizaciones realizadas en el código para que tardara menos tiempo. En primer lugar, la carga de los datos no se realizaba según el código de la Figura 6-1, sino que se hacía mediante la librería mnist y se utilizaba el siguiente código: Figura 6-3: Carga de datos del MNIST con librería mnist Es cierto que si se realiza de esta manera el código queda más compacto y comprensible, sin embargo, el gran inconveniente que presenta es que en cada ejecución ha de descargar dicho set de datos. Esto ralentiza bastante el código, sobre todo cuando se dispone de los archivos binarios ya descargados, y hace imposible su ejecución si no existe una conexión a internet. Por tanto, debido a estos motivos se decidió cambiar esto y utilizar el código 67 de la Figura 6-1 para optimizar el proceso y que no dependa de una conexión a internet. Todas las diferencias explicadas en el apartado anterior se desconocían cuando se empezó a realizar la implementación en Python. Dado que al ejecutar el código los resultados eran completamente diferentes a los esperados, para comprobar los fallos existentes se procedió a realizar una exhaustiva depuración. Para ello, se generó un vector theta en MATLAB, se guardó con el comando save y se pasó este archivo a Python, donde con un sencillo código (el de la Figura 6-4) se cargaba dicho vector. Teniendo los mismos datos iniciales en MATLAB y en Python, se podía ir línea por línea comprobando los resultados en ambos programas para comprobar los puntos de fallo. Figura 6-4: Carga de un archivo .mat Gracias a estas sencillas líneas de código, se pudo depurar el programa y conocer aquellos puntos en los que el resultado obtenido no coincidía con el de MATLAB. De esta manera se pudo conocer, por ejemplo, que el método por defecto de indexación de Python es C. Hasta ahora se ha hablado de cómo se descubrió el problema que presentaba la función reshape para matrices multidimensionales, pero no se ha hablado de cómo solucionar dicho problema. La primera idea que se tuvo para resolver esto fue el uso de una función creada a mano en la que, mediante el uso de bucles, permitiera pasar de una dimensión de origen específica a una de destíno específica. Por ejemplo, si se quería redimensionar una matriz de 2 dimensiones a una de 4, se usaba la siguiente función: Figura 6-5: Función para redimensionar matriz 2D a 4D Se creó una función específica para cada reshape necesario de una forma análoga a la mostrada en la Figura anterior. Sin embargo, esto realmente es muy poco óptimo ya que el uso de tantos bucles ralentizaba el código significativamente, aunque el resultado final era el correcto. Lo deseable sería poder expresar la matriz como en MATLAB, es decir, (rows,columns,channels) para poder utilizar order=’F’ en reshape y que realizara la operación correctamente. Tras una extensa investigación sobre las diferentes funciones relacionadas con reshape, se encontró una que permitía realizar el comportamiento buscado: swapaxes. La función swapaxes permite intercambiar las dimensiones de las matrices, de manera que si se tiene una matriz de 4×4×2 y se utiliza la orden swapaxes(0,2) se obtiene una matriz de 2×4×4. Para ilustrar mejor el funcionamiento de esta función aplicada al problema que nos atañe, si se supone una matriz 𝑨 de dimensiones 2×2×2×3, es decir, (rows,columns,numFilters,numImages) y suponinedo que 𝑨 vale: 68 FILTER 1 FILTER 2 IMAGE 1 (1 2 3 4) (5 6 7 8) IMAGE 2 (910 11 12) (13 14 15 16) IMAGE 3 (17 18 19 20) (21 22 23 24) Tabla 6-2: Matriz de 4 dimensiones La matriz que se desea obtener mediante el método de indexación Fortran es: 𝑻= ( 1 9 17 311 19 210 18 412 20 513 21 715 23 614 22 816 24 ) (6-3) La matriz obtenida en Python mediante el comando B=np.reshape(A,(-1,3),order=’F’) resulta ser: 𝑩= ( 119 14 9 7 22 17 15 4 523 12 13 220 21 10 8 318 16 11 624 ) (6-4) Dicha matriz, como se explicó anteriormente mezcla los elementos de las distintas imágenes y filtros, por lo que no es correcto el resultado. Sin embargo, si se pone en práctica el comando swapaxes de manera que se cambien las dimensiones en Python a (rows,columns,numFilters,numImages) y luego se redimensionaliza, es decir, se usa el comando C=A.swapaxes(0,2).swapaxes(1,3).swapaxes(2,3).reshape((-1,3),order=’F’), se obtiene: 69 𝑪= ( 1 9 17 311 19 210 18 412 20 513 21 715 23 614 22 816 24 ) (6-5) Como se puede observar, la matriz ahora obtenida se corresponde con la matriz deseada de la ecuación (6-3). Gracias a este comando, se puede realizar sin problemas una indexación tipo Fortran para matrices multidimensionales de una forma bastante compacta concatenando las órdenes. Este método se usará para aquellos casos en los que haya que redimensionalizar matrices, sustituyendo así a las funciones con bucles antes explicadas. Una vez ha quedado explicado todo lo relativo a la función reshape solo resta explicar cómo se ha realizado la operación pooling, pues no es tan sencillo como en MATLAB. El comando equivalente a pooledImage = convolved(1:poolDim:end, 1:poolDim:end) en Python sería pooledImage = convolved[0:None:poolDim] [0:None:poolDim], sin embargo, ésto no funciona y da un resultado no deseable. Para solventar dicho problema, en primera instancia se pensó en dos bucles anidados que recogieran de la matriz convolved los valores adecuados y los almacenara en pooledImage tal y como se observa en la siguiente imagen: Figura 6-6: Pooling con bucles anidados El resultado era el que debería ser, pero, debido a que el uso de dichos bucles ralentizaba el código de manera significativa, se buscaron otras formas de realizar la operación de pooling. Asímismo. se probó también con una función perteneciente a una librería, pero, ya que el tiempo de ejecución era el mismo que con los dos bucles anidados, se descartó dicha función. La solución final adoptada vino de la mano del comando swapaxes explicado anteriormente. Se almacenó en un vector los índices necesarios a recoger de la matriz convolved y se almacenó en pooledImage los valores correspondientes a dichos índices mediante los comandos: ind=range(convDim) ind=ind[0:None:poolDim] pooledImage=convolved[ind].swapaxes(0,1)[ind].swapaxes(0,1) 70 De esta manera se recogen los resultados de la operación pooling de una manera rápida y compacta gracias a la concatenación. La última modificación realizada en el código es que guarde en un archivo de texto el vector theta si la precisión alcanzada con la red supera el 97% para no tener que entrenarla de nuevo: Figura 6-7: Almacenamiento de parámetros si la precisión supera el 97% 6.4 Resultados obtenidos Una vez explicadas las diferencias existentes con respecto al lenguaje de programación MATLAB y las soluciones adoptadas para solventar dichos problemas, el resto del código se implementa de forma análoga a MATLAB, pero con la sintaxis de Python. Si se abre un terminal en Python y uno se dirige al directorio donde está contenido el código, al ejecutar python run_train.py, éste se ejecutará. Comenzará entrenando a la red durante unas 700 iteraciones y, tras esto, calculará la precisión de la red, así como el tiempo invertido, con el añadido de que, si la precisión está por encima del 97%, guardará los parámetros en un archivo de texto. Al ejecutar los pasos anteriores, se obtiene: Figura 6-8: Resultados obtenido en Python En la figura anterior se puede observar que la precisión obtenida es del 97.25%, la cual está por encima de lo obtenido en MATLAB y moviéndose en el rango de 97 y 97.4%, por lo que puede decirse que la implementación obtenida en Python ha sido un éxito. Con respecto al tiempo invertido que es aproximadamente de 18 minutos, es superior al tiempo que tarda MATLAB en entrenar a la red, sin embargo, hay que tener varias cosas en cuenta. La primera de ellas es que este tiempo se ha conseguido después de optimizar el código de la forma descrita en el apartado 6.3, ya que al principio el tiempo era del orden de 30 minutos. Lo segundo es que se está ejecutando sobre una máquina virtual, por lo que, si se instalara Ubuntu en el mismo ordenador y se ejecutara el código, el tiempo de ejecución sería bastante menor ya que aprovecharía más los recursos del ordenador. Por último, Python también ofrece la posibilidad de usar la GPU para acelerar los cálculos (recurso que no se está utilizando) de manera que podría disminuir el tiempo de ejecución y, si se combina ésto con el uso de Ubuntu como Sistema Operativo principal, el tiempo disminuiría aún más. Por tanto, si se realizaran las mejoras descritas anteriormente, el tiempo de ejecución estaría a la par que el de MATLAB y la diferencia entre ambos tanto en precisión como en duración sería mínima. Gracias a esto, las únicas diferencias entre ambos a la hora de realizar una CNN residirían en la propia sintaxis de cada lenguaje y la aplicación para la que se va a implementar, siendo el resto de elementos análogos. 71 7 CONCLUSIONES Tal y como se ha podido observar a lo largo del presente proyecto, las redes neuronales resultan muy interesantes en cuanto a problemas de clasificación se refieren, en concreto las CNN para el reconocimiento de imágenes. La posibilidad de obtener una precisión superior al 97% y de poder añadir clases o imágenes a clasificar sin tener que modificar solo algunos parámetros de la red, hacen que las redes neuronales sean una herramienta muy poderosa y útil. Además, se han probado diferentes tipos de redes tanto en MATLAB como Python, y en todas ellas, independientemente del lenguaje de programación, han dado unos resultados más que válidos. Si bien la programación en MATLAB resulta más intuitiva que en Python debido al completo álgebra matricial que presenta y la disposición de las dimensiones de las matrices, la implementación en Python conociendo estas diferencias es igual de sencilla. El único inconveniente de las redes neuronales, sobre todo cuando el número de capas comienza a incrementarse, es la necesidad de utilizar equipos relativamente potentes (al menos a la hora de entrenar a la red), ya que la ejecución del código en un portátil de hace aproximadamente 4 años era demasiado lenta (alrededor de 2 horas). Sin embargo, ésto es algo que todas aquellas personas que usan redes neuronales conocen y aceptan. Con este proyecto se pretendía brindar una introducción a las redes neuronales, su funcionamiento y las matemáticas detrás de ellas. Este campo está en constante crecimiento e investigación, pero gracias a este proyecto, el lector podrá tener una idea general de cómo se han llevado a cabo las redes más importantes a día de hoy, sobre todo aquellas relacionadas con el reconocimiento de imágenes, y asienta las bases a la hora de crear redes más complejas. 72 73 8 LÍNEAS DE TRABAJO FUTURAS Debido a la gran variedad de redes neuronales existentes a día de hoy, existen muchas líneas de investigación que podrían servir como continuación a este proyecto. A continuación, se presentan las más interesantes: • Creación de una red neuronal que, en vez de reconocer imágenes, sea capaz de reconocer patrones de sonido para, por ejemplo, reconocer cual es la persona que está hablando por una grabación. • Creación de una CNN cuyo set de imágenes sea de unas dimensiones más grandes, de manera que esto fuerce el uso de más capas en la red. Esto se podría complementar con, por ejemplo, tomar las imágenes a mano mediante una cámara (de señales de tráfico, por ejemplo) y que un determinado sistema actúe en función de que detecte la red. • Usar diferentes metodologías para entrenar a la red, ya que en este proyecto siempre se usó el aprendizaje supervisado. Sin embargo, existen otros tipos de aprendizaje como el aprendizaje no supervisado o el aprendizaje por refuerzo. Este último es muy interesante, pues la red va aprendiendo a medida que prueba diferentes opciones (como probar que acción en un videojuego es más beneficiosa en una situación y repetirlo en situaciones análogas). 80 81 APÉNDICE II: CÓDIGOS DE REGRESIÓN LOGÍSTICA Código 1 Script ex1b_logreg.m Descripción Carga los datos, ejecuta una regresión logística y muestra los resultados clear all; clc; close all; addpath ../common addpath ../common/minFunc_2012/minFunc addpath ../common/minFunc_2012/minFunc/compiled % Load the MNIST data for this exercise. % train.X and test.X will contain the training and testing images. % Each matrix has size [n,m] where: % m is the number of examples. % n is the number of pixels in each image. % train.y and test.y will contain the corresponding labels (0 or 1). binary_digits = true; [train,test] = ex1_load_mnist(binary_digits); % Add row of 1s to the dataset to act as an intercept term. train.X = [ones(1,size(train.X,2)); train.X]; test.X = [ones(1,size(test.X,2)); test.X]; % Training set dimensions m=size(train.X,2); n=size(train.X,1); % Train logistic regression classifier using minFunc options = struct('MaxIter', 100); % First, we initialize theta to some small random values. theta = rand(n,1)*0.001; % % Implement batch logistic regression in logistic_regression_vec.m using % MATLAB's vectorization features to speed up your code. Compare the running % tic; theta=minFunc(@logistic_regression_vec, theta, options, train.X, train.y); fprintf('Optimization took %f seconds.\n', toc); % Print out training accuracy. tic; accuracy = binary_classifier_accuracy(theta,train.X,train.y); fprintf('Training accuracy: %2.1f%%\n', 100*accuracy); % Print out accuracy on the test set. accuracy = binary_classifier_accuracy(theta,test.X,test.y); fprintf('Test accuracy: %2.1f%%\n', 100*accuracy); 82 Código 2 Script logistic_regression_vec.m Descripción Calcula la función coste y el gradiente de la regresión logística function [f,g] = logistic_regression_vec(theta, X,y) % % Arguments: % theta - A column vector containing the parameter values to optimize. % X - The examples stored in a matrix. % X(i,j) is the i'th coordinate of the j'th example. % y - The label for each example. y(j) is the j'th example's label. % m=size(X,2); n=size(X,1); % initialize objective value and gradient. f = 0; g = zeros(size(theta)); % % Compute the logistic regression objective function and gradient % using vectorized code. (It will be just a few lines of code!) % Store the objective function value in 'f', and the gradient in 'g'. % y_hat=logsig(theta'*X); f=-sum(y.*log(y_hat)+(1-y).*log(1-y_hat)); g=X*(y_hat-y)'; end 83 APÉNDICE III: CÓDIGOS DE REGRESIÓN SOFTMAX Código 1 Script ex1c_logreg.m Descripción Carga los datos, ejecuta una regresión softmax y muestra los resultados clear all; clc; addpath ../common addpath ../common/minFunc_2012/minFunc addpath ../common/minFunc_2012/minFunc/compiled % Load the MNIST data for this exercise. % train.X and test.X will contain the training and testing images. % Each matrix has size [n,m] where: % m is the number of examples. % n is the number of pixels in each image. % train.y and test.y will contain the corresponding labels (0 to 9). binary_digits = false; num_classes = 10; [train,test] = ex1_load_mnist(binary_digits); % Add row of 1s to the dataset to act as an intercept term. train.X = [ones(1,size(train.X,2)); train.X]; test.X = [ones(1,size(test.X,2)); test.X]; train.y = train.y+1; % make labels 1-based. test.y = test.y+1; % make labels 1-based. % Training set info m=size(train.X,2); n=size(train.X,1); % Train softmax classifier using minFunc options = struct('MaxIter', 200); % Initialize theta. We use a matrix where each column corresponds to a class, % and each row is a classifier coefficient for that class. % Inside minFunc, theta will be stretched out into a long vector (theta(:)). % We only use num_classes-1 columns, since the last column is always assumed 0. theta = rand(n,num_classes-1)*0.001; % Call minFunc with the softmax_regression_vec.m file as objective. % % Implement batch softmax regression in the softmax_regression_vec.m % file using a vectorized implementation. % tic; theta(:)=minFunc(@softmax_regression_vec, theta(:), options, train.X, train.y); fprintf('Optimization took %f seconds.\n', toc); theta=[theta, zeros(n,1)]; % expand theta to include the last class. % Print out training accuracy. tic; accuracy = multi_classifier_accuracy(theta,train.X,train.y); fprintf('Training accuracy: %2.1f%%\n', 100*accuracy); % Print out test accuracy. accuracy = multi_classifier_accuracy(theta,test.X,test.y); fprintf('Test accuracy: %2.1f%%\n', 100*accuracy); 84 Código 2 Script softmax_regression_vec.m Descripción Calcula la función coste y el gradiente de la regresión softmax function [f,g] = softmax_regression(theta, X,y) % % Arguments: % theta - A vector containing the parameter values to optimize. % In minFunc, theta is reshaped to a long vector. So we need to % resize it to an n-by-(num_classes-1) matrix. % Recall that we assume theta(:,num_classes) = 0. % % X - The examples stored in a matrix. % X(i,j) is the i'th coordinate of the j'th example. % y - The label for each example. y(j) is the j'th example's label. % m=size(X,2); n=size(X,1); % theta is a vector; need to reshape to n x num_classes. theta=reshape(theta, n, []); num_classes=size(theta,2)+1; % initialize objective value and gradient. f = 0; g = zeros(size(theta)); % % Compute the softmax objective function and gradient using vectorized code. % Store the objective function value in 'f', and the gradient in 'g'. % Before returning g, make sure you form it back into a vector with g=g(:); % y_hat=exp(theta'*X); % It’s the numerator y_hat=[y_hat; ones(1,size(y_hat,2))]; % Last row set to 1 to have 10 classes y_hat_sum=sum(y_hat,1); % It’s the denominator y_hat_sum=y_hat_sum-1; % As last row doesn’t count, substract it from the total sum coc=bsxfun(@rdivide,y_hat,y_hat_sum); A=log(coc); ind=sub2ind(size(A),y,1:size(y_hat,2)); A(end,:)=0; % Last row set to 0 because the prob that y belongs to 10 is null f=-sum(A(ind)); B=zeros(size(A)); B(ind)=1; g=-X*(B-coc)'; g=g(:,1:end-1); % Last row doesn’t count g=g(:); % make gradient a vector for minFunc end 85 APÉNDICE IV: CÓDIGOS DE LA RED NEURONAL MULTICAPA Código 1 Script run_train.m Descripción Carga los datos, crea una red neuronal multicapa, la entrena y muestra los resultados clear all; clc; tic; % runs training procedure for supervised multilayer network % softmax output layer with cross entropy loss function %% setup environment % experiment information % a struct containing network layer sizes etc ei = []; % add common directory to your path for % minfunc and mnist data helpers addpath ../common; addpath(genpath('../common/minFunc_2012/minFunc')); %% load mnist data [data_train, labels_train, data_test, labels_test] = load_preprocess_mnist(); %% populate ei with the network architecture to train % ei is a structure you can use to store hyperparameters of the network % the architecture specified below should produce 100% training accuracy % You should be able to try different network architectures by changing ei % only (no changes to the objective function code) % dimension of input features ei.input_dim = 784; % number of output classes ei.output_dim = 10; % sizes of all hidden layers and the output layer ei.layer_sizes = [256, ei.output_dim]; % scaling parameter for l2 weight regularization penalty ei.lambda = 0; % which type of activation function to use in hidden layers % feel free to implement support for only the logistic sigmoid function ei.activation_fun = 'logistic'; %% setup random initial weights stack = initialize_weights(ei); params = stack2params(stack); %% setup minfunc options options = []; options.display = 'iter'; options.maxFunEvals = 1e6; options.Method = 'lbfgs'; %% run training [opt_params,opt_value,exitflag,output] = minFunc(@supervised_dnn_cost,... params,options,ei, data_train, labels_train); %% compute accuracy on the test and train set [~, ~, pred] = supervised_dnn_cost( opt_params, ei, data_test, [], true); 86 [~,pred] = max(pred); acc_test = mean(pred'==labels_test); fprintf('test accuracy: %f\n', acc_test*100); [~, ~, pred] = supervised_dnn_cost( opt_params, ei, data_train, [], true); [~,pred] = max(pred); acc_train = mean(pred'==labels_train); fprintf('train accuracy: %f\n', acc_train*100); toc; Código 2 Script load_preprocess_mnist.m Descripción Carga los datos del set MNIST y hace las transformaciones pertinentes function [data_train, labels_train, data_test, labels_test] = load_preprocess_mnist() %% TODO ensure this is consistent with common loaders % assumes relative paths to the common directory % assumes common directory on paty for access to load functions % adds 1 to the labels to make them 1-indexed data_train = loadMNISTImages('../common/train-images-idx3-ubyte'); labels_train = loadMNISTLabels(['../common/train-labels-idx1-ubyte']); labels_train = labels_train + 1; data_test = loadMNISTImages('../common/t10k-images-idx3-ubyte'); labels_test = loadMNISTLabels(['../common/t10k-labels-idx1-ubyte']); labels_test = labels_test + 1; Código 3 Script initialize_weights.m Descripción Inicializa los pesos y bias, es decir, los parámetros de la red function [ stack ] = initialize_weights( ei ) %INITIALIZE_WEIGHTS Random weight structures for a network architecture % eI describes a network via the fields layerSizes, inputDim, and outputDim % % This uses Xavier's weight initialization tricks for better backprop % See: X. Glorot, Y. Bengio. Understanding the difficulty of training % deep feedforward neural networks. AISTATS 2010. %% initialize hidden layers stack = cell(1, numel(ei.layer_sizes)); for l = 1 : numel(ei.layer_sizes) if l > 1 prev_size = ei.layer_sizes(l-1); else prev_size = ei.input_dim; end; cur_size = ei.layer_sizes(l); % Xaxier's scaling factor 87 s = sqrt(6) / sqrt(prev_size + cur_size); stack{l}.W = rand(cur_size, prev_size)*2*s - s; stack{l}.b = zeros(cur_size, 1); end Código 4 Script stack2params.m Descripción Convierte los parámetros de la forma pesos y bias a vector function [params] = stack2params(stack) % Converts a "stack" structure into a flattened parameter vector and also % stores the network configuration. This is useful when working with % optimization toolboxes such as minFunc. % % [params, netconfig] = stack2params(stack) % % stack - the stack structure, where stack{1}.w = weights of first layer % stack{1}.b = weights of first layer % stack{2}.w = weights of second layer % stack{2}.b = weights of second layer % ... etc. % This is a non-standard version of the code to support conv nets % it allows higher layers to have window sizes >= 1 of the previous layer % If using a gpu pass inParams as your gpu datatype % Setup the compressed param vector params = []; for d = 1:numel(stack) % This can be optimized. But since our stacks are relatively short, it % is okay params = [params ; stack{d}.W(:) ; stack{d}.b(:) ]; % Check that stack is of the correct form assert(size(stack{d}.W, 1) == size(stack{d}.b, 1), ... ['The bias should be a *column* vector of ' ... int2str(size(stack{d}.W, 1)) 'x1']); % no layer size constrain with conv nets if d < numel(stack) assert(mod(size(stack{d+1}.W, 2), size(stack{d}.W, 1)) == 0, ... ['The adjacent layers L' int2str(d) ' and L' int2str(d+1) ... ' should have matching sizes.']); end end end 88 Código 5 Script params2stack.m Descripción Convierte los parámetros de la formavectorial a pesos y bias function stack = params2stack(params, ei) % Converts a flattened parameter vector into a nice "stack" structure % for us to work with. This is useful when you're building multilayer % networks. % % stack = params2stack(params, netconfig) % % params - flattened parameter vector % ei - auxiliary variable containing % the configuration of the network % % Map the params (a vector into a stack of weights) depth = numel(ei.layer_sizes); stack = cell(depth,1); % the size of the previous layer prev_size = ei.input_dim; % mark current position in parameter vector cur_pos = 1; for d = 1:depth % Create layer d stack{d} = struct; hidden = ei.layer_sizes(d); % Extract weights wlen = double(hidden * prev_size); stack{d}.W = reshape(params(cur_pos:cur_pos+wlen-1), hidden, prev_size); cur_pos = cur_pos+wlen; % Extract bias blen = hidden; stack{d}.b = reshape(params(cur_pos:cur_pos+blen-1), hidden, 1); cur_pos = cur_pos+blen; % Set previous layer size prev_size = hidden; end end 89 Código 6 Script supervised_dnn_cost.m Descripción Calcula la probabilidad de pertenecer a cada clase, la función coste y el gradiente function [ cost, grad, pred_prob] = supervised_dnn_cost( theta, ei, data, labels, pred_only) %SPNETCOSTSLAVE Slave cost function for simple phone net % Does all the work of cost / gradient computation % Returns cost broken into cross-entropy, weight norm, and prox reg % components (ceCost, wCost, pCost) %% default values po = false; if exist('pred_only','var') po = pred_only; end; %% reshape into network numHidden = numel(ei.layer_sizes) - 1; numSamples = size(data, 2); hAct = cell(numHidden+1, 1); gradStack = cell(numHidden+1, 1); stack = params2stack(theta, ei); %% forward prop for l = 1 : numHidden if l > 1 hAct{l} = stack{l}.W * hAct{l - 1} + repmat(stack{l}.b, 1,numSamples); else hAct{l} = stack{l}.W * data + repmat(stack{l}.b, 1, numSamples); end hAct{l} = sigmoid(hAct{l}); end l = numHidden+1; y_hat = stack{l}.W * hAct{l - 1} + repmat(stack{l}.b, 1, numSamples); y_hat = exp(y_hat); hAct{l} = bsxfun(@rdivide, y_hat, sum(y_hat, 1)); [pred_prob pred_labels] = max(hAct{l}); pred_prob = hAct{l}; %% return here if only predictions desired. if po cost = -1; ceCost = -1; wCost = -1; numCorrect = -1; grad = []; return; end; %% compute cost y_hat = log(hAct{numHidden+1}); index = sub2ind(size(y_hat), labels', 1:numSamples); ceCost = -sum(y_hat(index))./(size(data,2)); %% compute gradients using backpropagation % Cross entroy gradient targets = zeros(size(hAct{numHidden+1})); % numLabels * numSamples targets(index) = 1; gradInput = hAct{numHidden+1} - targets; for l = numHidden+1 : -1 : 1 if l > numHidden gradFunc = ones(size(gradInput)); else gradFunc = hAct{l} .* (1 - hAct{l}); end gradOutput = gradInput .* gradFunc; if l > 1 gradStack{l}.W = gradOutput * hAct{l-1}'./(size(data,2)); else 96 % Backpropagating the error w.r.t the softmax layer is as usual. To % backpropagate through the pooling layer, you will need to upsample the % error with respect to the pooling layer for each filter and each image. % Use the kron function and a matrix of ones to do this upsampling % quickly. targets=zeros(size(probs)); targets(ind)=1; errorsOutputSoftmax=(probs-targets)/numImages; % DIVIDED BY numImages??? errorsInputSoftmax=Wd'*errorsOutputSoftmax; errorsOutputPool=reshape(errorsInputSoftmax,outputDim,outputDim,numFilters,numImages); errorsInputPool=zeros(convDim,convDim,numFilters,numImages); for imageNum=1:numImages for filterNum=1:numFilters errorsInputPool(:,:,filterNum,imageNum)=(kron(errorsOutputPool(:,:,filterNum,imageNum),ones(poolD im)))/(poolDim^2); end end errorsInputConv=errorsInputPool.*activations.*(1-activations); %%====================================================================== %% STEP 1d: Gradient Calculation % After backpropagating the errors above, we can use them to calculate the % gradient with respect to all the parameters. The gradient w.r.t the % softmax layer is calculated as usual. To calculate the gradient w.r.t. % a filter in the convolutional layer, convolve the backpropagated error % for that filter with each image and aggregate over images. Wd_grad=errorsOutputSoftmax*(activationsPooled'); bd_grad=sum(errorsOutputSoftmax,2); for imageNum=1:numImages for filterNum=1:numFilters e=errorsInputConv(:,:,filterNum,imageNum); errorsInputConv(:,:,filterNum,imageNum)=rot90(e,2); end end for filterNum=1:numFilters for imageNum=1:numImages Wc_grad(:,:,filterNum)=Wc_grad(:,:,filterNum)+conv2(images(:,:,imageNum),errorsInputConv(:,:,filt erNum,imageNum),'valid'); end end for filterNum=1:numFilters e=errorsInputConv(:,:,filterNum,:); bc_grad(filterNum)=sum(e(:)); end %% Unroll gradient into grad vector for minFunc grad = [Wc_grad(:) ; Wd_grad(:) ; bc_grad(:) ; bd_grad(:)]; end 97 Código 4 Script computeNumericalGradient.m Descripción Calcula el gradiente de una segunda forma para verificar que el proporcionado por la función anterior es correcto function numgrad = computeNumericalGradient(J, theta) % numgrad = computeNumericalGradient(J, theta) % theta: a vector of parameters % J: a function that outputs a real-number. Calling y = J(theta) will return the % function value at theta. % Initialize numgrad with zeros numgrad = zeros(size(theta)); %% Implement numerical gradient checking, and return the result in numgrad. % You should write code so that numgrad(i) is (the numerical approximation to) the % partial derivative of J with respect to the i-th input argument, evaluated at theta. % I.e., numgrad(i) should be the (approximately) the partial derivative of J with % respect to theta(i). % % You will probably want to compute the elements of numgrad one at a time. epsilon = 1e-4; for i =1:length(numgrad) oldT = theta(i); theta(i)=oldT+epsilon; pos = J(theta); theta(i)=oldT-epsilon; neg = J(theta); numgrad(i) = (pos-neg)/(2*epsilon); theta(i)=oldT; if mod(i,100)==0 fprintf('Done with %d\n',i); end; end; end 98 Código 5 Script minFuncSGD.m Descripción Entrena la CNN utilizando el algoritmo SGD con Momentum function [opttheta] = minFuncSGD(funObj,theta,data,labels,... options) % Runs stochastic gradient descent with momentum to optimize the % parameters for the given objective. % % Parameters: % funObj - function handle which accepts as input theta, % data, labels and returns cost and gradient w.r.t % to theta. % theta - unrolled parameter vector % data - stores data in m x n x numExamples tensor % labels - corresponding labels in numExamples x 1 vector % options - struct to store specific options for optimization % % Returns: % opttheta - optimized parameter vector % % Options (* required) % epochs* - number of epochs through data % alpha* - initial learning rate % minibatch* - size of minibatch % momentum - momentum constant, defualts to 0.9 %%====================================================================== %% Setup assert(all(isfield(options,{'epochs','alpha','minibatch'})),... 'Some options not defined'); if ~isfield(options,'momentum') options.momentum = 0.9; end; epochs = options.epochs; alpha = options.alpha; minibatch = options.minibatch; m = length(labels); % training set size % Setup for momentum mom = 0.5; momIncrease = 20; velocity = zeros(size(theta)); %%====================================================================== %% SGD loop it = 0; for e = 1:epochs % randomly permute indices of data for quick minibatch sampling rp = randperm(m); for s=1:minibatch:(m-minibatch+1) it = it + 1; % increase momentum after momIncrease iterations if it == momIncrease mom = options.momentum; end; % get next randomly selected minibatch mb_data = data(:,:,rp(s:s+minibatch-1)); mb_labels = labels(rp(s:s+minibatch-1)); % evaluate the objective function on the next minibatch [cost grad] = funObj(theta,mb_data,mb_labels); 99 % Add in the weighted velocity vector to the % gradient evaluated above scaled by the learning rate. % Then update the current weights theta according to the % sgd update rule velocity=mom*velocity+alpha*grad; theta=theta-velocity; fprintf('Epoch %d: Cost on iteration %d is %f\n',e,it,cost); end; % aneal learning rate by factor of two after each epoch alpha = alpha/2.0; end; opttheta = theta; end Código 6 Script cnnParamsToStack.m Descripción Transforma los parámetros de la red en forma vector a pesos y bias function [Wc, Wd, bc, bd] = cnnParamsToStack(theta,imageDim,filterDim,... numFilters,poolDim,numClasses) % Converts unrolled parameters for a single layer convolutional neural % network followed by a softmax layer into structured weight % tensors/matrices and corresponding biases % % Parameters: % theta - unrolled parameter vectore % imageDim - height/width of image % filterDim - dimension of convolutional filter % numFilters - number of convolutional filters % poolDim - dimension of pooling area % numClasses - number of classes to predict % % % Returns: % Wc - filterDim x filterDim x numFilters parameter matrix % Wd - numClasses x hiddenSize parameter matrix, hiddenSize is % calculated as numFilters*((imageDim-filterDim+1)/poolDim)^2 % bc - bias for convolution layer of size numFilters x 1 % bd - bias for dense layer of size hiddenSize x 1 outDim = (imageDim - filterDim + 1)/poolDim; hiddenSize = outDim^2*numFilters; %% Reshape theta indS = 1; indE = filterDim^2*numFilters; Wc = reshape(theta(indS:indE),filterDim,filterDim,numFilters); indS = indE+1; indE = indE+hiddenSize*numClasses; Wd = reshape(theta(indS:indE),numClasses,hiddenSize); indS = indE+1; indE = indE+numFilters; bc = theta(indS:indE); bd = theta(indE+1:end); end 100 101 APÉNDICE VI: CÓDIGOS EN PYTHON TRAS OPTIMIZACIÓN Código 1 Script cnnTrain.py Descripción Carga los datos, crea una red neuronal convolucional, la entrena y muestra los resultados from __future__ import division import os import numpy as np from pdb import set_trace as pb import cnnInitParams import cnnCost import computeNumericalGradient import minFuncSGD import time from numpy import linalg as LA import struct t=time.time() """ Convolutional Neural Network """ # Instructions # ------------ # # This file contains code that helps you get started in building a single # layer convolutional network. In this execirse, you will only # need to modify cnnCost.py and minFuncSGD.py. You will not need to # modify this file """ ====================================================================== """ """ STEP 0: Initialize Parameters and Load Data """ # Here we initialize some parameters used for the exercise # Configurtion imageDim = 28 # Dimension of the input images numClasses = 10 # Number of classes (MNIST images fall into 10 classes) filterDim = 9 # Filter size for convolutional layer numFilters = 20 # Numbre of filters for convolutional layer poolDim = 2 # Pooling dimension, (should divide imageDim-filterDim+1) # Load MNIST train path="." fname_img = os.path.join(path,'train-images-idx3-ubyte') fname_lbl = os.path.join(path,'train-labels-idx1-ubyte') with open(fname_img, 'rb') as fimg: magic, num, rows, cols = struct.unpack(">IIII",fimg.read(16)) img = np.fromfile(fimg, dtype=np.uint8).reshape(-1,rows,cols) images = img/255 with open(fname_lbl, 'rb') as flbl: magic, num = struct.unpack(">II",flbl.read(8)) labels = np.fromfile(flbl, dtype=np.int8) # Initialize Parameters theta=cnnInitParams.cnnInitParams(imageDim,filterDim,numFilters,poolDim,numClasses) """ ====================================================================== """ """ STEP 1: Implement ConvNet Objective """ 102 # Implement the funcion cnnCost.py """ ====================================================================== """ """ STEP 2: Gradient Check """ # Use the file computeNumericalGradient.py to check the gradient # calculation for you cnnCost.py function. You may need to add the # appropiate path or copy the file to this directory DEBUG=True # Set this to true to check the gradient if DEBUG: # To speed up gradient checking, we will use a reduced network and # a debugging data set print "Starting Gradient Check" db_numFilters = 2 db_filterDim = 9 db_poolDim = 5 db_images = images[0:10][:][:] db_labels = labels[0:10] db_theta = cnnInitParams.cnnInitParams(imageDim,db_filterDim,db_numFilters,db_poolDim,numClasses) cost, grad, preds = cnnCost.cnnCost(db_theta,db_images,db_labels,numClasses,db_filterDim,db_numFilters,db_poolDim,Fal se) # Check gradients numGrad = computeNumericalGradient.computeNumericalGradient(lambda x: cnnCost.cnnCost(x, db_images,db_labels,numClasses,db_filterDim,db_numFilters,db_poolDim,False),db_theta) diff = LA.norm(numGrad-grad)/LA.norm(numGrad+grad) # Should be small. In our implementation, these values are usually # less than 1e-9. print "Difference is",diff assert(diff < 1e-9),"Difference too large. Check your gradient computation again" elapsed=time.time()-t print "Elapsed time in checking is",elapsed,"seconds" """ ================================================================== """ """ STEP 3: Learn Parameters """ # Implement minFuncSGD.py, then train the model print "Starting Training" options = np.zeros((4,1)) options[0] = 3 options[1] = 256 options[2] = 1e-1 options[3] = 0.95 opttheta = minFuncSGD.minFuncSGD(lambda x, y, z: cnnCost.cnnCost(x,y,z,numClasses,filterDim,numFilters,poolDim,False),theta,images,labels,options) print "The sum of theta is:",np.sum(opttheta) """ =================================================================== """ """ STEP 4: Test """ # Test the performance of the trained model using the MNIST test set. Your # accuracy should be above 97 after 3 epochs of training print "Testing trained network" """testImages, testLabels = mndata.load_testing() testImages = np.reshape(testImages,(-1,imageDim,imageDim)) testImages = testImages/255""" fname_img = os.path.join(path,'t10k-images-idx3-ubyte') fname_lbl = os.path.join(path,'t10k-labels-idx1-ubyte') with open(fname_img, 'rb') as fimg: magic, num, rows, cols = struct.unpack(">IIII",fimg.read(16)) img = np.fromfile(fimg, dtype=np.uint8).reshape(-1,rows,cols) testImages = img/255 with open(fname_lbl, 'rb') as flbl: magic, num = struct.unpack(">II",flbl.read(8)) 103 testLabels = np.fromfile(flbl, dtype=np.int8) _, cost, preds = cnnCost.cnnCost(opttheta,testImages,testLabels,numClasses,filterDim,numFilters,poolDim,True) acc = np.sum(preds==testLabels)/(np.shape(testLabels)[0]) # Accuracy should be around 97.4% after 3 epochs print "Accuracy is",acc*100,"%" print "The number of successfully predicted images is",np.sum(preds==testLabels) elapsed=time.time()-t print "Elapsed time is",elapsed/60,"minutes" if acc>0.97: print "Saving optimized theta in file savedata.txt" np.savetxt('savedata.txt',opttheta) print "Done" Código 2 Script cnnInitParams.py Descripción Inicializa los parámetros de la red import math import numpy as np def cnnInitParams(imageDim,filterDim,numFilters,poolDim,numClasses): # Initialize parameters for a single layer convolutional neural # network followed by a softmax layer # # Parameters: # imageDim - height/width of image # filterDim - dimension of convolutional filter # numFilters - number of convolutional filters # poolDim - dimension of pooling area # numClasses - number of classes to predict # # # Returns: # theta - unrolled parametervector with initialized weights """ Initialize parameters randomly based on layer sizes """ assert (filterDim < imageDim),"filterDim must be less than imageDim" Wc = 1e-1*np.random.normal(0,1,(numFilters,filterDim,filterDim)) outDim = imageDim - filterDim + 1 # dimension of convolved image # Assume outDim is multiple of poolDim assert (outDim % poolDim == 0),"poolDim must divide outDim" outDim = outDim/poolDim hiddenSize = outDim**2*numFilters # We'll choose weights uniformly from the interval [-r, r] r=math.sqrt(6)/math.sqrt(numClasses+hiddenSize+1) Wd=np.random.uniform(0,1,(numClasses,hiddenSize))*2*r-r bc=np.zeros((numFilters,1)) bd=np.zeros((numClasses,1)) # Convert weights and bias gradients to the vector form. # This step will "unroll" (flatten and concatenate together) all # your parameters into a vector, which can then be used with minFunc Wc=np.reshape(Wc,(-1,1)) 104 Wd=np.reshape(Wd,(-1,1)) bc=np.reshape(bc,(-1,1)) bd=np.reshape(bd,(-1,1)) theta=np.append(Wc,Wd) theta=np.append(theta,bc) theta=np.append(theta,bd) theta=np.reshape(theta,(-1,1)) return theta Código 3 Script cnnCost.py Descripción Calcula la probabilidad de pertenecer a cada clase, la función coste y el gradiente from __future__ import division import numpy as np import scipy from scipy import signal import sigmoid import cnnParamsToStack from pdb import set_trace as pb from numpy import linalg as LA import sys def cnnCost(theta,images,labels,numClasses,filterDim,numFilters,poolDim,pred): # Calcualte cost and gradient for a single layer convolutional # neural network followed by a softmax layer with cross entropy # objective. # # Parameters: # theta - unrolled parameter vector # images - stores images in imageDim x imageDim x numImges # array # numClasses - number of classes to predict # filterDim - dimension of convolutional filter # numFilters - number of convolutional filters # poolDim - dimension of pooling area # pred - boolean only forward propagate and return # predictions # # # Returns: # cost - cross entropy cost # grad - gradient with respect to theta (if pred==False) # preds - list of predictions for each example (if pred==True) imageDim = np.shape(images)[1] # height/width of image numImages = np.shape(images)[0] # number of images """ Reshape parameters and setup gradient matrices """ # Wc is filterDim x filterDim x numFilters parameter matrix # bc is the corresponding bias # Wd is numClasses x hiddenSize parameter matrix where hiddenSize # is the number of output units from the convolutional layer # bd is corresponding bias Wc, Wd, bc, bd = cnnParamsToStack.cnnParamsToStack(theta,imageDim,filterDim,numFilters,poolDim,numClasses) # Same sizes as Wc,Wd,bc,bd. Used to hold gradient w.r.t above params. Wc_grad = np.zeros(np.shape(Wc)) 105 Wd_grad = np.zeros(np.shape(Wd)) bc_grad = np.zeros(np.shape(bc)) bd_grad = np.zeros(np.shape(bd)) """ ====================================================================== """ """ STEP 1a: Forward Propagation """ # In this step you will forward propagate the input through the # convolutional and subsampling (mean pooling) layers. You will then use # the responses from the convolution and pooling layer as the input to a # standard softmax layer. """ Convolutional Layer """ # For each image and each filter, convolve the image with the filter, add # the bias and apply the sigmoid nonlinearity. Then subsample the # convolved activations with mean pooling. Store the results of the # convolution in activations and the results of the pooling in # activationsPooled. You will need to save the convolved activations for # backpropagation. convDim = imageDim-filterDim+1 # dimension of convolved output outputDim = (convDim)//poolDim # dimension of subsampled output # convDim x convDim x numFilters x numImages tensor for storing activations activations = np.zeros((numImages,numFilters,convDim,convDim)) # outputDim x outputDim x numFilters x numImages tensor for storing # subsampled activations activationsPooled = np.zeros((numImages,numFilters,outputDim,outputDim)) """ YOUR CODE HERE """ W = np.ones((poolDim,poolDim)) ind=range(convDim) ind=ind[0:None:poolDim] for imageNum in xrange(numImages): for filterNum in xrange(numFilters): convolvedImage = np.zeros((convDim,convDim)) filter = Wc[filterNum][:][:] filter = np.rot90(np.rot90(filter)) im = images[:][:][imageNum] convolvedImage = scipy.signal.convolve2d(im,filter,'valid') convolvedImage = convolvedImage+np.tile(bc[filterNum],np.shape(convolvedImage)) convolvedImage = sigmoid.sigmoid(convolvedImage) activations[imageNum][filterNum][:][:] = convolvedImage pooledImage = np.zeros((outputDim,outputDim)) filter = W im = convolvedImage convolved = scipy.signal.convolve2d(im,filter,'valid')/(poolDim**2) pooledImage=convolved[ind].swapaxes(0,1)[ind].swapaxes(0,1) activationsPooled[imageNum][filterNum][:][:] = pooledImage # Reshape activations into 2-d matrix, hiddenSize x numImages, # for Softmax layer activationsPooledAfter=np.zeros((outputDim**2*numFilters,numImages)) activationsPooledAfter=activationsPooled.swapaxes(0,2).swapaxes(1,3).swapaxes(2,3).reshap e((-1,numImages),order='F') """ Softmax Layer """ # Forward propagate the pooled activations calculated above into a # standard softmax layer. For your convenience we have reshaped # activationPooled into a hiddenSize x numImages matrix. Store the # results in probs. # numClasses x numImages for storing probability that each image belongs to # each class. probs = np.zeros((numClasses,numImages)) softmax=np.dot(Wd,activationsPooledAfter)+np.tile(bd,(1,numImages)) softmax=np.exp(softmax) probs=softmax/np.sum(softmax,axis=0) """ ====================================================================== """ """ STEP 1b: Calculate Cost """ # In this step you will use the labels given as input and the probs # calculate above to evaluate the cross entropy objective. Store your # results in cost. 112