Implementación de técnicas para análisis cluster robusto en torno a subespacios afines

Grado en Estad ´ ıstica Trabajo Fin de Grado Implementaci´on de t´ecnicas para An´alisis Cluster robusto en torno a subespacios afines Autor D. Jes´us Fern´andez Iglesias Tutor D. Luis ´ Angel Garc´ıa Escudero a 2 Resumen La constante generaci´on de conjuntos de datos masivos que se produce en la actualidad ha provocado que el desarrollo de t´ecnicas de aprendizaje autom´atico capaces de extraer conocimiento ´util de dicha informaci´on sea un campo del conocimiento en auge y en constante desarrollo. En muchos de estos problemas, las observaciones no tienen asociadas ning´un tipo de etiqueta, categor´ıa o clase, ´unicamente se dispone de los propios datos. Por tanto, la b´usqueda de patrones ocultos en los mismos se torna una tarea fundamental. Con ese fin, el paradigma de aprendizaje no supervisado ofrece una amplia gama de procedimientos que permiten el estudio y agrupaci´on de objetos en base a sus similitudes. En la intersecci´on de la incesante creaci´on de conjuntos de datos enormes y el aprendizaje no supervisado surge la necesidad de desarrollar e implementar procedimientos computacionalmente eficientes para poder aplicar estas t´ecnicas de aprendizaje no supervisado y en particular de la aplicaci´on de t´ecnicas de an´alisis cluster. En este trabajo se estudia y se desarrollan versiones computacionalmente eficientes de un procedimiento de an´alisis cluster robusto entorno a subespacios afines. Un enfoque robusto al an´alisis cluster evita que unas pocas observaciones at´ıpicas pueden condicionar de manera muy negativa la detecci´on correcta de clusters. La metodolog´ıa desarrollada examina varias opciones de implementaci´on, explorando el enfoque secuencial, el paralelizado y el h´ıbrido sacando partido a varios lenguajes de programaci´on, adem´as de realizar los correspondientes an´alisis de eficiencia computacional para determinar qu´e versi´on es la m´as adecuada. Adem´as, un ejemplo de aplicaci´on real del procedimiento desarrollado es mostrado en el ´ambito de la segmentaci´on de im´agenes. 3 a 4 Abstract The constant generation of massive data sets nowadays has made the development of machine learning techniques capable of extracting useful knowledge from such information a booming and constantly developing field of knowledge. In many of these problems, the observations do not have any kind of label, category or class associated with them; only the data itself is available. Therefore, the search for hidden patterns in the data becomes a fundamental task. To that end, the unsupervised learning paradigm offers a wide range of procedures that allow the study and clustering of objects based on their similarities. At the intersection of the incessant creation of huge datasets and unsupervised learning arises the need to develop and implement computationally efficient procedures to be able to apply these unsupervised learning techniques and in particular the application of cluster analysis techniques. In this paper we study and develop computationally efficient versions of a robust cluster analysis procedure around affine subspaces. A robust approach to cluster analysis avoids that a few outlier observations can condition in a very negative way the correct detection of clusters. The developed methodology examines several implementation options, exploring sequential, parallelized and hybrid approaches taking advantage of several programming languages, as well as performing the corresponding computational efficiency analysis to determine which version is the most suitable. In addition, an example of real application of the developed procedure is shown in the field of image segmentation. 5 a 6 Agradecimientos A los profesores del Grado en Estad´ıstica que me han acompa˜nado durante estos a˜nos por transmitirme la pasi´on y el conocimiento necesarios para formar mi futuro profesional. A Luis ´ Angel, por la constante dedicaci´on y ayuda que han puesto en este proyecto y resolver todas las cuestiones que le he ido planteando de manera r´apida y eficaz. A mis compa˜neros de carrera por haber compartido esta traves´ıa a mi lado, superando juntos todas las dificultades surgidas. Y por ´ultimo, a mi familia, en espacial a mis padres, por inculcarme los valores del sacrificio y el esfuerzo tan necesarios, y ser los pilares que me han sostenido durante toda esta etapa. a 7 a 8 ´ INDICE GENERAL ´ Indice general 1. Introducci´on 14 1.1. Contextualizaci´on.................................... 14 1.2. Contenidosatratar................................... 15 1.3. Resumen global de resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.4. Herramientasutilizadas................................. 16 2. An´alisis Cluster y Componentes Principales 18 2.1. M´etodos de Clustering jer´arquicos........................... 18 2.2. M´etodos de Clustering no jer´arquicos o particionales . . . . . . . . . . . . . . . . 19 2.2.1. K-medias .................................... 19 2.3. An´alisis en Componentes Principales . . . . . . . . . . . . . . . . . . . . . . . . . 20 3. Robustez 24 3.1. Mediaymediana .................................... 25 3.2. K-mediasrecortadas .................................. 28 3.3. Aplicaci´on de los recortes al ACP . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4. Clustering robusto entorno a subespacios lineales 36 4.1. Agrupaciones robustas entorno a subespacios afines d-dimensionales . . . . . . . . 37 4.1.1. Caso I: problema te´orico o poblacional . . . . . . . . . . . . . . . . . . . . 37 4.1.2. Caso II: problema emp´ırico o muestral . . . . . . . . . . . . . . . . . . . . 38 4.2. Descripci´on del algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5. Implementaci´on 41 5.1. Mejorasa˜nadidas .................................... 41 5.1.1. Subespacios afines de dimensiones diferentes . . . . . . . . . . . . . . . . . 41 5.1.2. Estimador de m´ınimo determinante de la covarianza y distancia robusta de Mahalanobis................................... 45 5.1.3. Funcionalidad gr´afica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.1.3.1. 2D................................... 50 5.1.3.2. 3D................................... 51 5.2. Pseudoc´odigo ...................................... 54 5.3. Computaci´onparalela ................................. 58 5.3.1. B´usqueda de subespacios afines inicial . . . . . . . . . . . . . . . . . . . . 59 5.3.2. Iteraci´on completa de las inicializaciones m´as prometedoras . . . . . . . . . 62 5.3.3. Flujo general del programa . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.4. Integraci´on de RyC++ ................................. 66 5.4.1. B´usqueda de subespacios afines inicial . . . . . . . . . . . . . . . . . . . . 68 9 1.3 Resumen global de resultados El ´ultimo cap´ıtulo de esta memoria est´a dedicado a exponer las conclusiones del trabajo y a plantear l´ıneas futuras de investigaci´on. 1.3. Resumen global de resultados En este trabajo se ha logrado desarrollar una versi´on computacionalmente muy eficiente del algoritmo que agrupa observaciones entorno a subespacios afines de manera robusta. Dicha versi´on adapta un enfoque h´ıbrido, donde el c´odigo es desarrollado en dos lenguajes de programaci´on que interact´uan e intercambian informaci´on durante la ejecuci´on del procedimiento. Los lenguajes que cooperan son RyC++, siendo el primero el encargado de llevar a cabo el flujo general del programa y el segundo el que realiza las partes computacionalmente m´as costosas del algoritmo. La versi´on h´ıbrida mejora con creces tanto a una versi´on que aprovecha el poder de la computaci´on paralela tambi´en desarrollada y detallada en este trabajo como a una funci´on disponible en el CRAN de R. Los experimentos desarrollados y los resultados obtenidos son expuestos en el Cap´ıtulo 5. Aprovechando la eficiencia de la implementaci´on desarrollada, el procedimiento se aplica a la segmentaci´on de im´agenes, las cuales han de ser aplanadas constituyendo en s´ı mismas un conjunto de datos inmenso para tama˜nos de im´agenes comunes. Se aprecia como la aplicaci´on del algoritmo con distintas configuraciones genera m´ascaras de segmentaci´on diferentes, remarcando o suavizando partes de intensidad m´as d´ebil de la imagen, las cu´ales pueden ser muy ´utiles en multitud de dominios de aplicaci´on. 1.4. Herramientas utilizadas A continuaci´on se muestra la relaci´on de herramientas utilizadas durante la realizaci´on del proyecto: R: es un lenguaje de programaci´on de distribuci´on libre ampliamente usada por la comunidad estad´ıstica y en ciencia de datos, est´a orientado al desarrollo de software estad´ıstico y al an´alisis de datos. En este trabajo, buena parte de los gr´aficos se han creado con esta herramienta, adem´as del desarrollo ´ıntegro de dos versiones del procedimiento implementado en el Cap´ıtulo 5, una secuencial y una paralelizada. C++: es un lenguaje de programaci´on generalista que extiende al lenguaje Cpara permitir la manipulaci´on de objetos. En este trabajo se ha utilizado para implementar parte de la versi´on h´ıbrida (combinando dos lenguajes) del an´alisis cluster robusto entorno a subespacios afines. Se trata de uno de los lenguajes de programaci´on m´as eficientes que existen. 16 1 INTRODUCCI´ ON RStudio: es un entorno de desarrollo integrado para el lenguaje de programaci´on R. Contiene herramientas para la gesti´on, el trazado y la depuraci´on del espacio de trabajo. Los c´odigos RyC++ han sido desarrollados en este IDE (Integrated Development Environment). TeXStudio: es un editor de L A T EXde c´odigo abierto y multiplataforma. La decisi´on de utilizar L A T EXpara la redacci´on de esta memoria es por su caracter´ıstica de producir documentos con alta calidad tipogr´afica. 17 2. An´alisis Cluster y Componentes Principales A fecha de realizaci´on de este documento, se generan, aproximadamente, 2.5 quintillones de gigabytes de informaci´on al d´ıa. Si queremos obtener conocimiento ´util a partir de esta informaci´on mediante la aplicaci´on de t´ecnicas de ML (Machine Learning) nos encontraremos que la gran mayor´ıa de las observaciones no est´an etiquetadas en ninguna clase o no se dispone de una variable respuesta medida en ellas. Es por ello que el aprendizaje no supervisado es cada vez m´as usado, cobrando especial importancia en ´areas como la detecci´on de anomal´ıas, fraudes y defectos, visualizaci´on, reducci´on de la dimensionalidad, sistemas de recomendaci´on... Las t´ecnicas com´unmente aplicadas en estos problemas pueden agruparse principalmente en dos familias: M´etodos param´etricos: en este caso se asume una distribuci´on de probabilidad subyacente a los datos, es decir, los individuos provienen de una poblaci´on que sigue una distribuci´on probabil´ıstica con un conjunto de par´ametros fijos. El aprendizaje no supervisado param´etrico requiere de la construcci´on de modelos, t´ıpicamente de mezcla gaussiana (Gaussian Mixture Models), y del uso de algoritmos EM (Expectation-Maximization). M´etodos no param´etricos: al contrario que en los m´etodos param´etricos, no se requiere hacer suposiciones acerca de la distribuci´on subyacente a los datos, por tanto, en ocasiones, se denominan distribution-free methods. Dentro de la inmensidad de m´etodos y aproximaciones que se encuentran dentro del aprendizaje no supervisado, en este trabajo se pondr´a el foco en el clustering. Su objetivo consiste en agrupar observaciones, compuestas por un conjunto de atributos o variables, ´unicamente atendiendo a criterios basados en sus diferencias y semejanzas. Los conjuntos obtenidos tras realizar el agrupamiento se denominan clusters, conglomerados o grupos. El principal inconveniente que nos encontramos es la dificultad para definir lo que es un cluster. La imprecisi´on en dicha definici´on ha tra´ıdo consigo gran cantidad de m´etodos. 2.1. M´etodos de Clustering jer´arquicos En esta familia de m´etodos, la pertenencia a un cluster en un nivel espec´ıfico de jerarqu´ıa condiciona su posible pertenencia en niveles superiores. La manera en la que se construye la jerarqu´ıa genera dos tipos de t´ecnicas de clustering: Aglomerativos: su construcci´on es de abajo hacia arriba. Se parte de clusters formados por una ´unica observaci´on. Iterativamente, se van agrupando para formar clusters m´as grandes. 18 2 AN´ ALISIS CLUSTER Y COMPONENTES PRINCIPALES Divisivos: su construcci´on es de arriba hacia abajo. Se parte de un cluster formado por la totalidad de los individuos. Iterativamente se va dividiendo dando lugar a conglomerados m´as peque˜nos seg´un se va descendiendo en la jerarqu´ıa. 2.2. M´etodos de Clustering no jer´arquicos o particionales Esta familia de m´etodos obtienen una ´unica partici´on del espacio muestral. Las t´ecnicas de agrupamiento pertenecientes a esta categor´ıa son las m´as famosas y utilizadas. El m´aximo exponente es el algoritmo de las k-medias, aunque otros como DBSCAN tambi´en tienen su importancia. Respecto al n´umero de grupos en los que se desea fragmentar el conjunto de datos, la estrategia a seguir sigue dos l´ıneas diferenciadas. Algunos de los algoritmos, como es el caso de k-medias, necesitan como condici´on previa a su uso que se especifique el n´umero de clusters. Otros, como es el caso de DBSCAN, calculan autom´aticamente el n´umero de agrupamientos, sin necesidad de que sean especificados previamente. Ambas aproximaciones tienen sus ventajas e inconvenientes, ofreciendo variedad y diversidad de soluciones para multitud de problemas a los que son aplicables este tipo de t´ecnicas. En particular, en este trabajo, se va a exponer con m´as detalle el algoritmo de las k-medias, ya que guarda relaci´on con el algoritmo que se va a implementar en el Cap´ıtulo 5. Concretamente, se trata de un caso particular, donde las dimensiones de los ksubespacios afines entorno a los que agrupar son 0 (es decir, kpuntos o centroides) y no recorta ninguna observaci´on (ninguna observaci´on es tratada como un punto at´ıpico). 2.2.1. K-medias El algoritmo de las k-medias es un m´etodo de an´alisis cluster particional que trata de segmentar el conjunto de datos X={x1, ..., xn}en kgrupos disjuntos de tal manera que cada observaci´on sea asignada a un ´unico grupo (Figura 2.1). Los pasos de los que consta el algoritmo de las k-medias son los siguientes: 1. Especificar el n´umero de clusters ky seleccionar kobservaciones aleatorias del conjunto de datos. Estos kpuntos son los denominados centroides iniciales. 2. Calcular la distancia eucl´ıdea entre cada observaci´on pertenciente a Xy cada uno de los k centroides. 3. Asignar cada observaci´on al cluster del que el centroide m´as cercano es representante. 4. Calcular los nuevos centroides como la media de todas las observaciones pertenecientes al mismo grupo. 5. Repetir los pasos 2, 3 y 4 hasta que los centroides no var´ıen en 2 iteraciones consecutivas. 19 2.3 An´alisis en Componentes Principales Figura 2.1: Agrupamiento realizados por k-medias (k= 3). Sin embargo, este algoritmo presenta diversos problemas. Entre ellos destacan la sensibilidad a la elecci´on inicial de los centroides, la necesidad de prefijar el n´umero, k, de clusters, el hecho de que ´unicamente se encontrar´an clusters convexos y la falta de robustez frente a observaciones at´ıpicas. Como ya se ha mencionado, el algoritmo de las k-medias es un caso particular del algoritmo que se va a detallar en este trabajo, con la configuraci´on de par´ametros (d1, ..., dk) = (0, ..., 0) (representando djla dimensi´on del subespacio af´ın aproximante j-´esimo) y α= 0 (representando la proporci´on de observaciones que ser´an recortadas), donde kes el n´umero de clusters. 2.3. An´alisis en Componentes Principales Hay muchos problemas candidatos a la aplicaci´on de modelos de aprendizaje autom´atico en los que es necesario tratar con una cantidad enorme de variables medidas. La aplicaci´on directa de algunos de estos algoritmos a conjuntos de datos con esta caracter´ıstica puede causar numerosos problemas. Por un lado, al disponer de muchas variables explicativas, el n´umero de par´ametros del modelo ser´a, por lo general, consecuentemente muy elevado. Para que las estimaciones de dichos par´ametros sean suficientemente precisas se requiere una gran cantidad de observaciones, lo cual muchas veces no es posible. Por otra parte, es m´as que sabido el peligro que conlleva el incluir demasiadas variables ex- 20 2 AN´ ALISIS CLUSTER Y COMPONENTES PRINCIPALES plicativas en un modelo, puede darse el problema del sobreajuste. Recordar que el sobreajuste aparece cuando se permite una gran complejidad del modelo ajustado, y ´este parece ajustar de manera ´optima los datos con los que se han estimado los par´ametros. Sin embargo, la capacidad de generalizaci´on de un modelo con sobreajuste es mucho menor que lo que hace indicar su funcionamiento sobre los datos iniciales. Otro problema encontrado pr´acticamente a diario en la estad´ıstica es el de la multicolinealidad. La presencia de un n´umero excesivo de variables explicativas facilita que muchas de ellas tengan una gran correlaci´on entre s´ı, es decir, multicolinealidad. Por tanto, se puede decir que la presencia de un n´umero muy elevado de variables explicativas act´ua como catalizador de la multicolinealidad. Adem´as, los requerimientos de memoria y capacidad de c´omputo son m´as ambiciosos a medida que aumenta el tama˜no de los conjuntos de datos, ya sea debido a un aumento de los casos, de las variables explicativas, o de ambos. Por estos motivos, en muchos contextos y dominios de aplicaci´on, es necesaria la utilizaci´on de t´ecnicas de reducci´on de la dimensionalidad como paso previo al uso de otros m´etodos. Algunos clasificadores del paradigma de aprendizaje por ensembles como Random Forest, los filtros de baja varianza o los filtros de alta correlaci´on pueden utilizarse con el fin de alcanzar este objetivo. Sin embargo, la t´ecnica m´as conocida y utilizada, y la que se va a tratar con detalle en este trabajo, es el An´alisis en Componentes Principales (ACP). Las componentes principales de un conjunto de observaciones representadas en Rpson un conjunto de pvectores directores ortogonales entre s´ı (Figura 2.2). La d-´esima componente principal vendr´ıa dada por el d-´esimo vector que mejor recoge la variabilidad de la muestra, siendo ortogonal a los d−1 vectores anteriores. Por tanto, el primer vector ser´a aquel que represente la direcci´on que recoja la m´axima varianza de la nube de puntos. El fundamento te´orico b´asico que reside en este m´etodo es que la variabilidad total del conjunto de datos es invariante frente a la rotaci´on de ejes ortogonales y existen rotaciones que hacen que gran parte de la informaci´on contenida en nuestros datos se encuentre disponible en la proyecci´on de dichos datos en el subespacio generado por unos pocos vectores directores (Figura 2.3). De manera m´as intuitiva, el An´alisis en Componentes Principales puede verse como ajustar un elipsoide p-dimensional a los datos, donde los ejes del elipsoide representan las pcomponentes principales. Por tanto, si uno de los ejes del elipsoide es peque˜no, entonces la proporci´on de varianza recogida por ese eje tambi´en ser´a peque˜na. 21 2.3 An´alisis en Componentes Principales Figura 2.2: Representaci´on esquem´atica de las componentes principales z1yz2de una nube de puntos el´ıptica en IR2. Figura 2.3: Rotaci´on de los datos en el sistema de coordenadas de las componentes principales z1 yz2. Cada componente principal se obtiene por combinaci´on lineal de las variables originales. Sea (X1, X2, ..., Xp) un vector p-dimensional con las variables explicativas, entonces la primera componente principal Z1es la combinaci´on lineal de dichas variables que recoja la mayor cantidad de varianza: Z1=φ11X1+φ21X2+... +φp1Xp donde p X j=1 φ2 j1= 1 debido a que la combinaci´on lineal es normalizada. De forma an´aloga: 22 2 AN´ ALISIS CLUSTER Y COMPONENTES PRINCIPALES Zd=φ1dX1+φ2dX2+... +φpdXp ser´ıa la combinaci´on lineal de dichas variable que recoge la d-´esima mayor variabilidad. Los t´erminos φij se conocen como loadings, y son los pesos en las variables que caracterizan a la componente principal. Cuando es necesario interpretar el significado de una componente principal estos valores son muy ´utiles, ya que permiten conocer el peso que tiene cada variable independiente en la componente, ayudando a determinar que tipo de informaci´on recoge la misma. Para calcular el valor de los loadings, y por tanto la componente principal, es necesario resolver un problema de optimizaci´on para maximizar la variabilidad recogida por la componente. En la pr´actica, esta operaci´on es equivalente a obtener los autovectores de la matriz de varianzascovarianzas muestral del conjunto de datos ordenados seg´un los autovalores. Se calcula la matriz de varianzas-covarianzas del conjunto de datos X ∈ Rn×pcomo: S=1 n−1 n X i=1 (xi−¯x)(xi−¯x)T donde xies una fila de la matriz de datos Xy ¯xcontiene las medias aritm´eticas por columnas de las pvariables. Entonces, la descomposici´on en autovalores y autovectores de la matriz de varianzas-covarianzas puede hacerse de la siguiente manera: SV =Vλ donde Ves una matriz cuyas columnas son los autovectores de S, y λuna matriz diagonal con los autovalores. Entonces: S=VλV−1 Por tanto, la primera columna de Vasociada al autovalor m´as alto ser´ıa el autovector que nos marca la direcci´on de la primera componente principal. La columna asociada al segundo autovalor m´as alto marca la direcci´on de la segunda componente principal, y as´ı hasta completar las p componentes principales. Al igual que suced´ıa con el algoritmo de las k-medias, el An´alisis en Componentes Principales es otro caso particular del algoritmo que se trata e implementa en este trabajo. Las componentes principales se obtienen, como caso particular simple, del enfoque adoptado cuando k=1. Es decir, si no se busca ninguna estructura de clusters y todas las observaciones est´an en un mismos cluster. 23 3. Robustez Como ya se ha comentado al comienzo del Cap´ıtulo 2, diariamente se genera una cascada inmensa de datos. Una parte de esos datos se puede modelizar razonablemente bien por unas pocas distribuciones de probabilidad muy conocidas, permitiendo trabajar con ellos y explotarlos mediante la aplicaci´on de m´etodos que aprovechan las suposiciones a nivel poblacional realizadas sobre la distribuci´on subyacente que ha generado los datos. Por tanto, en esta situaci´on, los procedimientos desarrollados y estudiados en la estad´ıstica cl´asica son perfectamente v´alidos. Sin embargo, en aplicaciones de an´alisis reales de datos es muy cuestionable el poder justificar el cumplimiento de las hip´otesis sobre las distribuciones probabil´ısticas en las que se basan muchos de estos m´etodos. Algunos conjuntos de datos presentan ruido, es decir, las observaciones presentan una variabilidad que no puede ser explicada por los modelos que se tratan de ajustar. En general, es imposible evitar la presencia de ruido en un conjunto de datos, ya que se extraen de entornos de producci´on donde esta caracter´ıstica es innata. La presencia de un ruido o contaminaci´on excesivos puede da˜nar irreversiblemente un conjunto de datos, y hacerlo inutilizable. Otros conjuntos de datos se caracterizan por la presencia de valores at´ıpicos o outliers, observaciones significativamente diferentes del resto. Por ejemplo, estas observaciones at´ıpicas aparecen cuando nos enfrentamos a datos generados desde distribuciones con colas pesadas (heavy tailed) o incluso cuando aparecen puntos que, sin ser extremos, se separan del comportamiento mayoritario del resto de observaciones. La presencia de valores at´ıpicos puede deberse o bien a variabilidad intr´ınseca de la muestra o bien a un error experimental, es decir, un error de medici´on, de transcripci´on... En este ´ultimo caso, y si se comprueba el error, se suele optar por eliminar el outlier, siempre y cuando ese valor no pueda ser corregido. Los valores at´ıpicos pueden causar serios problemas en el an´alisis estad´ıstico de los datos. De ah´ı, surge la necesidad del desarrollo de t´ecnicas estad´ısticas que sean robustas frente a la presencia de outliers. Un m´etodo no robusto frente a outliers ofrecer´a un resultado diferente en funci´on de si los datos que toma como entrada presentan o no este tipo de puntos. La soluci´on calculada cuando hay presencia de valores at´ıpicos no suele ser adecuada. Es lo que se busca precisamente con el desarrollo de este tipo de t´ecnicas, un m´etodo robusto frente a outliers ofrecer´a un resultado similar, y razonable la mayor´ıa de las veces, tanto si los datos que toma como entrada presentan valores at´ıpicos como si no. Por tanto, se puede concluir que un m´etodo estad´ıstico robusto es un procedimiento que mantiene un buen funcionamiento para datos provenientes de una amplia gama de distribuciones de probabilidad, si el modelo que ha generado los datos es pr´oximo a las hip´otesis supuestas ideales, aunque est´as hip´otesis no se cumplan exactamente. Los m´etodos no robustos pueden tener un funcionamiento muy malo incluso con muestras de datos que han sido generadas por modelos muy 24 3 ROBUSTEZ parecidos al modelo asumido pero que no coinciden con dicho modelo. 3.1. Media y mediana Una primera y sencilla aproximaci´on para comprender qu´e diferencias hay entre estad´ısticos que carecen o presentan robustez estad´ıstica es la diferencia entre la media y la mediana. Sea x1, x2, ..., x40 una muestra proveniente de una distribuci´on normal multivariante tomando valores en R3. Para simular una situaci´on en la que varios casos de un conjunto de datos presentan ruido, a las 8 ´ultimas observaciones se les a˜nadir´a en cada una de las 3 componentes una cantidad aleatoria de ruido proveniente de una distribuci´on uniforme continua (Figura 3.1). x1, x2, ..., x32 ∼N   0 0 0 ,  1 0 0 0 1 0 0 0 1   x33, x34, ..., x40 ∼N   0 0 0 ,  1 0 0 0 1 0 0 0 1  + (U(15,20), U(15,20), U(15,20))0 Al generar los datos siempre supondremos que todas las distribuciones que aparezcan han sido generadas de forma independiente entre s´ı. Figura 3.1: Representaci´on de los datos x1, . . . , x40 25 3.3 Aplicaci´on de los recortes al ACP Figura 3.5: Agrupaciones realizadas por las k-medias recortadas (k= 2, α= 0.05). Bajo el primero de los grupos se engloban las 950 observaciones provenientes de la distribuci´on normal multivariante con vector de medias (0,0,0), mientras que en el segundo grupo se encuentran las 950 observaciones provenientes de la distribuci´on normal multivariante con vector de medias (5 2,5 2,5 2). Las 100 observaciones extra´ıdas aleatoriamente de la ´ultima distribuci´on normal multivariante, el 5 % del total, son dejadas sin asignar a ning´un grupo (graficadas en color negro), por lo que el algoritmo consigue aislarlas y que no intervengan en la determinaci´on de los centroides de los clusters. Este ejemplo sirve para mostrar como el procedimiento de las k-medias recortadas es m´as robusto a las observaciones at´ıpicas que el algoritmo de k-medias cl´asico. Las k-medias recortadas han probado contar con buenas propiedades de robustez como se puede ver en [2]. 3.3. Aplicaci´on de los recortes al ACP Al igual que sucede con el algoritmo de las k-medias, el An´alisis en Componentes Principales es otro m´etodo susceptible de ofrecer malos resultados cuando los datos sobre los que es aplicado presentan valores extremos. Para la determinaci´on de la direcci´on de m´axima varianza de un conjunto de observaciones todos los puntos influyen de la misma manera, y unas pocas observaciones alejadas de la nube de puntos principal pueden acabar desviando dicha direcci´on de la que cabr´ıa esperar ´unicamente teniendo en cuenta a la nube de puntos principal. Para ilustrar este fen´omeno, se dise˜na un experimento generando datos artificialmente al igual 32 3 ROBUSTEZ que se hizo en el caso de la k-medias. Por ejemplo, sup´ongase que se tienen 2 submuestras, cada una de ellas proveniente de una distribuci´on normal multivariante tal y como se especifica a continuaci´on: x1, x2, ..., x950 ∼N   0 0 0 ,  10 0 0 01 20 0 0 0 1 20    x951, x952, ..., x1000 ∼N   10 10 10 ,  10 0 0 0 10 0 0 0 10   El primer conjunto est´a compuesto por una muestra de 950 observaciones, y el segundo conjunto est´a compuesto por una muestra de 50 observaciones. Los 2 conjuntos se unen en un mismo conjunto de datos, teniendo as´ı 1000 puntos en R3provenientes de 2 distribuciones distintas. La representaci´on gr´afica del conjunto de datos se aprecia en la Figura 3.6. Figura 3.6: Representaci´on del conjunto de datos x1, . . . , x1000. Cuando se realiza un ACP sobre el conjunto de datos, se observa que el vector que marca la direcci´on de m´axima varianza de los datos es ~v =(0.8022, 0.4205, 0.4239), representado con una l´ınea azul en la Figura 3.7. Analizando gr´aficamente la direcci´on de m´axima varianza se aprecia como, lejos de recoger ´unicamente la informaci´on de la nube de puntos principal formada 33 3.3 Aplicaci´on de los recortes al ACP por las 950 observaciones provenientes de la primera distribuci´on normal, las 50 observaciones extra´ıdas de una normal con vector de medias (10,10,10) son determinantes a la hora de calcular el vector director. La recta cuya direcci´on es marcada por el vector est´a claramente inclinada hacia los 50 puntos provenientes de la segunda distribuci´on normal, con lo que se concluye que estas observaciones est´an influyendo de manera clara en la determinaci´on de la direcci´on de m´axima varianza. Figura 3.7: Direcci´on de m´axima varianza encontrada por ACP. Queda claro que es necesario disponer de un procedimiento de c´alculo de las componentes principales que sea resistente a las observaciones at´ıpicas, y que en este ejemplo sea capaz de determinar la direcci´on de m´axima dispersi´on de la nube principal de puntos, sin que observaciones secundarias intervengan en dicha decisi´on. Para ello, se puede generalizar el algoritmo de las kmedias recortadas de tal manera que, en vez de buscar kgrupos, ´unicamente se busque uno para determinar la direcci´on de m´axima varianza. Las agrupaciones se realizar´an entorno a un subespacio af´ın de dimensi´on 1, una recta, en vez de un subespacio af´ın de dimensi´on 0 como ser´ıa un centroide. Este procedimiento de An´alisis en Componentes Principales recortadas es un caso particular del algoritmo que se detallar´a en el Cap´ıtulo 4, al igual que las k-medias recortadas. Concretamente, en el caso de las k-medias recortadas la configuraci´on de par´ametros ser´ıa (d1, ..., dk) = (0, ..., 0) (ksubespacios afines de dimensi´on 0 donde kes el n´umero de grupos) y α > 0 (la proporci´on de observaciones a recortar), mientras que en el caso del ACP con recortes la relaci´on de par´ametros ser´ıa d= 1 y α > 0. Si aplicamos el procedimiento con el par´ametro de recorte α= 5 % al conjunto de datos 34 3 ROBUSTEZ generado a partir de 2 distribuciones normales multivariantes donde el ACP no ofrec´ıa una soluci´on ´optima nos encontramos con que ahora la direcci´on de m´axima varianza encontrada si es la esperada, como se aprecia en la Figura 3.8. La direcci´on de m´axima variabilidad del conjunto de puntos encontrada es mucho m´as razonable que la propuesta por el la t´ecnica cl´asica de An´alisis en Componentes Principales, debido a que los 50 puntos provenientes de la distribuci´on normal multivariante con vector de medias (10,10,10) no son utilizados para determinar la recta que mejor resume la varianza de la muestra. Al fijar el par´ametro de recorte α= 0.05 se consigue que esas 50 observaciones no modifiquen la direcci´on de m´axima varianza que encontrar´ıa el ACP aplic´andose ´unicamente sobre los 950 puntos provenientes de la distribuci´on normal con vector de medias (0,0,0), y que representa el 95 % del conjunto de datos. Frente al vector que marca la direcci´on de m´axima varianza de los datos ~v =(0.8022, 0.4205, 0.4239) encontrado originalmente, el procedimiento robusto encuentra el vector ~v =(0.9999, -0.0037, 0.0008), la variaci´on y la mejor´ıa del ajuste es evidente. Figura 3.8: Direcci´on de m´axima varianza encontrada por ACP con recortes. 35 4. Clustering robusto entorno a subespacios lineales Bajo los m´etodos de an´alisis cluster no jer´arquicos es com´un que resida el concepto de formar grupos o conglomerados en torno a puntos centrales, los cu´ales representan el comportamiento t´ıpico de las observaciones pertenecientes al grupo del que dichos puntos son representantes. El m´aximo exponente de esta familia de m´etodos es el algoritmo de las k-medias, basado en el criterio de los m´ınimos cuadrados. Este criterio presenta mayoritariamente dos problemas. El primero de ellos es la falta de robustez estad´ıstica y que ya hemos comentado en secciones anteriores. El segundo de los problemas que presentan esta familia de m´etodos es que tienen tendencia a encontrar clusters esf´ericos debido a que realizan la suposici´on de que el comportamiento de un conglomerado de puntos puede ser resumido por una observaci´on central. Sin embargo, en ocasiones la presencia de agrupaciones en un conjunto de datos es debida a la existencia de relaciones entre las variables medidas que no tienen por qu´e responder a patrones esf´ericos. Por ejemplo, dada una muestra de observaciones, se pueden encontrar diferentes estructuras lineales como l´ıneas, planos o hiperplanos alrededor de los cu´ales las observaciones pueden ser agrupadas de una manera natural, sin forzar a crear clusters inducidos por criterios que favorecen el agrupamiento en estructuras esf´ericas. Los primeros intentos de realizar an´alisis cluster entorno a subespacios lineales vinieron por parte de Hosmer (1984), Lenstra et al. (1982) y Sp¨ath (1982), ajustando una mixtura de dos modelos de regresi´on lineal simple. Para el problema de estructuras lineales bidimensionales, aproximaciones alternativas fueron presentadas por Murtagh y Raftery (1984) y Phillips y Rosenfeld (1988). DeSarbo y Cron (1988) presentaron una soluci´on para un n´umero general de dimensiones y un n´umero arbitrario de grupos. Van Aelst et al. (2006) [3] abordaron el problema del agrupamiento lineal usando una aproximaci´on basada en regresi´on ortogonal. Dicho planteamiento obtuvo muy buenos resultados en diversos problemas donde no hab´ıa presencia de observaciones at´ıpicas. Bajo esta aproximaci´on y buscando un ´unica agrupaci´on, el problema se reduce al An´alisis en Componentes Principales cl´asico, que como ya se ha mostrado en esta memoria tambi´en sufre serios problemas de robustez. En muchas de las potenciales ´areas de aplicaci´on como la visi´on artificial, el reconocimiento de patrones o la tomograf´ıa la presencia de ruido en los conjuntos de datos es muy frecuente, por lo que urg´ıa mejorar las aproximaciones al agrupamiento entorno a estructuras lineales para que fuesen robustas frente al ruido y los outliers. 36 4 CLUSTERING ROBUSTO ENTORNO A SUBESPACIOS LINEALES 4.1. Agrupaciones robustas entorno a subespacios afines d-dimensionales Garc´ıa-Escudero et al. (2009) [4] introducen una metodolog´ıa que busca agrupamientos entorno a estructuras lineales en presencia de observaciones at´ıpicas, dotando de robustez al algoritmo de agrupamiento lineal basado en regresi´on ortogonal en [3]. El m´etodo est´a basado en la idea de impartial trimming (Gordaliza (1991) [5], Cuesta-Albertos et al. (1997) [1]), es decir, son los propios datos los que deciden cu´ales han de ser las observaciones que tienen que ser eliminadas. El m´etodo permite realizar agrupaciones entorno a subespacios afines generales d-dimensionales, no ´unicamente entorno a l´ıneas rectas. Adem´as, como ya se ha comentado, la robustez estad´ıstica se incluye de una manera muy natural utilizando el criterio de los m´ınimos cuadrados recortados. Dada una muestra x1, . . . , xnde observaciones en Rp, 0 ≤α < 1 (la proporci´on de observaciones que ser´an recortadas), d(la dimensi´on de los subespacios afines siendo 0 ≤d<p) y k∈N (el n´umero de grupos que se pretenden buscar), se busca la soluci´on al siguiente problema de minimizaci´on: m´ın Y⊂{x1,...,xn},#Y=[n(1−α)] m´ın {h1,...,hk}∈Ad 1 [n(1 −α)] X xi∈Y m´ın j=1,...,k{kxi−Prhj(xi)k2}!(1) donde Ad:= {h⊂Rp,hes un subespacio af´ın d-dimensional}yP rh(·) denota la proyecci´on ortogonal en h. Cualquier soluci´on H0={h0 1, . . . , h0 k}del problema induce a una partici´on de las observaciones no recortadas en kagrupaciones entorno a subespacios afines, de tal manera que la agrupaci´on Cjest´a formada por todas las observaciones no recortadas que est´an m´as cerca de h0 jque de cualquiera de los otros k−1 subespacios en H0. 4.1.1. Caso I: problema te´orico o poblacional Asumiendo que {x1, . . . , xn}es el resultado de una muestra aleatoria X1, . . . , Xnde una distribuci´on de probabilidad P, el problema emp´ırico o muestral en (1) admite un hom´ologo te´orico o poblacional. Sea Puna distribuci´on de probabilidad continua en Rp,α∈(0,1) y k∈N. Para cada H={h1, . . . , hk} ⊂ Ady cada conjunto Atal que P(A) = 1 −α, se mide la k-variaci´on entorno a Hdado Acomo: VA(H) := 1 1−αZA d(x, H)2dP(x) con d(x, H) = minj=1,...,kkx−P rhj(x)k. Entonces, se obtiene la k-variaci´on dado Ahaciendo una minimizaci´on en H: VA:= ´ınf H⊂Ad,#H=k{VA(H)} 37 4.2 Descripci´on del algoritmo y finalmente se obtiene la k-variaci´on recortada con recorte αminimizando en A: Vk,α := ´ınf A:P(A)=1−α(VA) Resolviendo este doble problema de minimizaci´on, se consigue un conjunto ´optimo A0yksubes- pacios afines ´optimos H0={h0 1, . . . , h0 k}tal que VA0(H0) = Vk,α Dado que se asume que no existe ninguna variable explicativa que ha de ser resumida por otro conjunto de variables se utilizan distancias ortogonales para medir las discrepancias entre los puntos y los hiperplanos. Esto es, se utiliza el criterio de la regresi´on ortogonal y no el de la regresi´on cl´asica. Adem´as, estas distancias no est´an escaladas, por lo que se asume igualdad de varianza a lo largo de todos los subespacios. Esta metodolog´ıa fue extendida en [6] al caso en que s´ı existe una variable respuesta privilegiada y se usan los residuales t´ıpicamente utilizados en An´alisis de Regresi´on. 4.1.2. Caso II: problema emp´ırico o muestral Siendo {Xn}nuna secuencia de vectores aleatorios independientes e igualmente distribuidos extra´ıdos de la distribuci´on P, la distribuci´on emp´ırica se define como: Pn(A) = 1 n n X i=1 IA(Xi) El problema original (1) tiene un desarrollo id´entico al mostrado en el caso poblacional, pero sustituyendo la distribuci´on desconocida Ppor la distribuci´on emp´ırica Pn. Fijados X1= x1, . . . , Xn=xn, el resultado de existencia de soluciones probada en el caso de distribuciones te´oricas Ppuede ser utilizada para derivar f´acilmente la existencia de soluciones en el caso emp´ırico. Existen formas ´optimas de dividir {x1, . . . , xn}en kgrupos de tal manera que el n´umero total de elementos en el total de los grupos sea [n(1 −α)]. Para cada una de las particiones, los k subespacios afines ´optimos se obtienen recurriendo a la regresi´on ortogonal de las observaciones en cada grupo. La referencia [4] muestra una prueba de la convergencia de las soluciones muestrales a las soluciones del problema poblacional, la convergencia entre subespacios afines ha de ser vista como la convergencia de distancias al origen y la posibilidad de elegir secuencias convergentes de los vectores generadores de los subespacios (spanning vectors). 4.2. Descripci´on del algoritmo El c´alculo de los ksubespacios afines α-recortados ´optimos en el problema emp´ırico tiene una alta complejidad computacional, debido a que para alcanzar ese ´optimo es necesario realizar una b´usqueda en un espacio combinatorio de subconjuntos, un espacio que generalmente es demasiado 38 4 CLUSTERING ROBUSTO ENTORNO A SUBESPACIOS LINEALES amplio. Por tanto, un algoritmo que encuentre la soluci´on ´optima siempre realizando una b´usqueda de fuerza bruta no es factible, y se abre el camino para un m´etodo que encuentre una aproximaci´on adecuada al ´optimo. El algoritmo que se va a describir y posteriormente implementar es una adaptaci´on del propuesto para realizar el c´alculo de las k-medias recortadas por Garc´ıa-Escudero et al. (2003) [7], y debe ser visto como una combinaci´on del algoritmo de las k-medias cl´asico y del algoritmo FASTMCD en Rousseeuw y Van Driessen (1999) [8] para calcular el estimador de m´ınimo determinante de la covarianza (MCD). En las k-medias recortadas y en FAST-MCD un paso de concentraci´on es requerido donde se seleccionan las observaciones con menores distancias eucl´ıdeas a sus respectivos centros, en este algoritmo se reservan las [n(1 −α)] observaciones con menor distancia ortogonal al subespacio m´as cercano de entre los ksubespacios afines calculados en la iteraci´on anterior. Tras este paso, se obtienen knuevos subespacios afines resolviendo kproblemas de regresi´on ortogonal. Dada una muestra {x1, . . . , xn}de observaciones ∈Rp, 0 ≤α < 1 (la proporci´on de observaciones que ser´an recortadas), d(la dimensi´on de los subespacios afines siendo 0 ≤d < p), Ad:= {h⊂Rp,hes un subespacio af´ın d-dimensional}yk∈N(el n´umero de grupos que se pretenden buscar), el algoritmo consta de los siguientes pasos: 1. Realizar un escalado de las variables para evitar problemas de precisi´on num´erica. El m´etodo de escalado elegido es robusto, dividiendo cada variable entre su MAD (Median Absolute Deviation). 2. Seleccionar ksubespacios afines iniciales en Ad. Por ejemplo, seleccionar aleatoriamente (d+ 1) ×kobservaciones del conjunto de datos y usarlos para obtener ksubespacios afines, donde cada uno es determinado por d+ 1 observaciones. El c´alculo del subespacio af´ın viene dado por el punto medio xj 0de los d+ 1 puntos y una matriz Uj 0cuyas columnas son los d autovectores unitarios asociados a los autovalores distintos de 0 de la matriz de varianzascovarianzas muestral de las d+ 1 observaciones. 3. Paso de concentraci´on. Sea H={h1, . . . , hk} ⊂ Adlos ksubespacios afines calculados en la iteraci´on anterior. a) Calcular las distancias di=d(xi, H), i = 1, . . . , n, entre cada observaci´on y el subespacio af´ın m´as pr´oximo de entre los ksubespacios afines obtenidos en la iteraci´on anterior. Determinar el conjunto Cque consiste en las [n(1 −α)] observaciones con menores distancias di, donde: d2 i= ´ınf j=1,...,kk{I−Uj 0(Uj 0)0}(xi−xj 0)k2 39 4.2 Descripci´on del algoritmo b) Hacer una partici´on del conjunto Cen C={C1, . . . , Ck}donde los puntos en Cjson aquellas observaciones que est´an m´as cercanas a hjque a cualquiera de los dem´as subespacios afines hlcon l6=j, es decir: Cj:= {xi∈C:d(xi, hj)2=di} c) Sea mjla media muestral de las observaciones en CjyUj 1una matriz cuyas columnas son los dautovectores asociados a los dmayores autovalores de la matriz de varianzas-covarianzas muestral de las observaciones en Cj. Los ksubespacios afines H={h1, . . . , hk}para la siguiente iteraci´on ser´an los ksubespacios afines que pasan por m1, . . . , mky son generados por los vectores disponibles en las columnas de U1 1, . . . , Uk 1respectivamente. 4. Repetir el paso de concentraci´on varias veces. Posteriormente, calcular el valor de la siguiente funci´on de coste: 1 [n(1 −α)] k X j=1 X xi∈Cj d2 i 5. Repetir los pasos 2, 3 y 4 varias veces, de tal manera que se inicialice el algoritmo con subespacios afines iniciales diferentes en cada comienzo. Guardar las soluciones que obtengan menores valores en la funci´on de coste, iterar completamente esas soluciones y elegir la ´optima en base al criterio de minimizar la funci´on de coste. Cuando el par´ametro de recorte αtoma el valor 0, entonces el algoritmo expuesto se reduce al agrupamiento lineal en Van Aelst et al. (2006) [3]. Cuando se busca agrupar las obervaciones entorno a un ´unico subespacio k= 1 se est´a realizando un An´alisis en Componentes Principales recortado en el caso de que α > 0 o un An´alisis en Componentes Principales cl´asico en el caso de que α= 0. Si se busca agrupar las observaciones entorno a ksubespacios de dimensi´on 0, k centroides, el m´etodo se transforma en las k-medias recortadas de la secci´on 3.2 en el caso de que α > 0 y en las k-medias cl´asicas en el caso de que α= 0. 40 5 IMPLEMENTACI´ ON 5. Implementaci´on En este cap´ıtulo se va a explicar el desarrollo de dos versiones del m´etodo presentado en el Cap´ıtulo 4 con el objetivo de optimizar la computaci´on de los subespacios afines. Una versi´on utiliza el poder de la computaci´on paralela en R, y otra la eficiencia del lenguaje C++. Una versi´on secuencial en Rtambi´en ha sido desarrollada pero no ser´a detallada al ser m´as sencilla y carecer de menos valor que las otras dos versiones comentadas. Adem´as, se han a˜nadido nuevos pasos y modificaciones para mejorar el algoritmo, obteniendo una versi´on m´as completa del mismo. Estas mejoras ser´an explicadas como paso previo al pseudoc´odigo desarrollado. 5.1. Mejoras a˜nadidas 5.1.1. Subespacios afines de dimensiones diferentes La primera de los mejoras consideradas es poder agrupar las observaciones entorno a subespacios afines de diferentes dimensiones. Se rompe con la limitaci´on de que todas las estructuras lineales que se busquen tengan que residir en el mismo subespacio vectorial. En la versi´on original del m´etodo, se buscaban H={h1, h2, . . . , hk} ⊂ Adsubespacios afines, donde h1∈Rp, h2∈Rp, . . . , hk∈Rpcon p < d. Ahora, adem´as de seguir contemplando el caso anterior, se permite buscar H={h1, h2, . . . , hk} ⊂ Adsubespacios afines en Rpcon dimensiones d1, . . . , dk(no necesariamente coincidentes) y dj< p para j= 1, . . . , k. Para ilustrar una situaci´on donde sea ´util esta generalizaci´on de agrupar entorno a subespacios afines de diferentes dimensiones, se dise˜na un experimento generando datos artificialmente provenientes de 2 distribuciones normales multivariantes y una mixtura de la distribuci´on normal y la distribuci´on uniforme, tal y como se especifica a continuaci´on: x1, x2, ..., x1000 ∼N   −2 −2 −2 ,  1 50 0 01 50 0 0 1 5    x1001, x1002, ..., x2000 ∼N   5 5 5 ,  10 0 0 01 20 0 0 0 1 20    xi,1∼N12,1 10,xi,2∼U(5,15) , xi,3∼U(−5,15) , i= 2001,...,6000 41 5.1 Mejoras a˜nadidas Figura 5.8: Elipses de tolerancia cl´asica y robusta sobre un mismo conjunto de datos. Aplicando esta t´ecnica al algoritmo desarrollado en esta secci´on, una vez se hayan estimado los par´ametros de localizaci´on y escala de los clusters, estos se utilizar´an para calcular la distancia con las observaciones mediante la distancia de Mahalanobis, que permite medir la distancia entre un punto y una distribuci´on o nube de puntos, relativa al tama˜no de la nube. Por tanto, para el m´etodo aqu´ı descrito, se calcular´an de manera robusta mediante el determinante de m´ınima covarianza los par´ametros de localizaci´on y escala de las knubes de puntos que representan los kclusters: ˆµCkyˆ ΣCk. Posteriormente, para cada observaci´on, se calcular´an sus distancias de Mahalanobis con las knubes de puntos formadas por los kclusters: dMah-Rob(xi, Ck) = q(xi−ˆµCk)0ˆ Σ−1 Ck(xi−ˆµCk) Por tanto, siendo x1, . . . , xn∈Rpel conjunto de observaciones y C={C1, . . . , Ck}los clusters encontrados por el m´etodo, la asignaci´on final en clusters vendr´a dada por: Cj:= {xi∈C:dMah-Rob(xi, Cj) = m´ın l=1,...,k dMah-Rob(xi, Cl)} Realizando este paso final, se consigue que cuando los subespacios lineales encontrados por el algoritmo provoquen la intersecci´on de dos o m´as subespacios afines entre s´ı, los resultados que se obten´ıan anteriormente sean corregidos, y las observaciones que se encuentran cerca de la intersecci´on ahora sean asignadas a la nube de puntos que conforme un cluster m´as cercana. En la Figura 5.9 se ven los resultados obtenidos tras aplicar la reasignaci´on final en base al MCD y la distancia robusta de Mahalanobis. Ahora, en la intersecci´on entre los dos subespacios afines, el comportamiento es el esperado, siendo las observaciones fronterizas asignadas a la nube de puntos m´as cercana (el plano) y no al hiperplano m´as pr´oximo siguiendo el criterio de la distancia ortogonal. 48 5 IMPLEMENTACI´ ON Figura 5.9: Agrupaciones realizadas (d=(1,2) y α= 0) tras reasignar los grupos mediante MCD y Mahalanobis. La funci´on desarrollada en Rpara realizar la funcionalidad explicada es la siguiente: allocateMahanalobis<- function (n,d, grupos ){ robustMahalanobis<-matrix(rep (0 ,n*length(d)), nrow = n, ncol =length(d)) for(i in 1: length(d)){ mcd<- covMcd(puntos[grupos==i,]) for(j in 1:n){ robustMahalanobis[j,i]<-mahalanobis( puntos [j ,] , mcd $ center,mcd $ cov) } } gruposRobustos<- apply (robustMahalanobis ,1,which .min) gruposRobustos [ grupos ==1000] <-1000 return(gruposRobustos) } 5.1.3. Funcionalidad gr´afica Como parte de la funcionalidad a˜nadida a la implementaci´on del algoritmo de clustering entorno a subespacios lineales se ha dise˜nado un sistema para representar gr´aficamente las agrupaciones realizadas y los subespacios ajustados en dimensiones bajas. La funcionalidad gr´afica est´a disponible cuando las observaciones del conjunto de datos se encuentran en R2yR3. A continuaci´on se detallan las dos opciones. 49 5.1 Mejoras a˜nadidas 5.1.3.1 2D Cuando se pretenden agrupar las observaciones utilizando ´unicamente la informaci´on proporcionada por dos variables los gr´aficos que se generan son, como es l´ogico, en dos dimensiones. En este caso, los espacios en los que agrupar las observaciones s´olo pueden ser de dimensi´on 0 y 1, y ser´an representados mediante puntos y rectas respectivamente. Un ejemplo de representaci´on gr´afica cuando los subespacios afines entorno a los que agrupar son de dimensi´on 0 se encuentra en la Figura 5.10. Se aprecia como el color en el que se encuentran las observaciones representa su grupo, teniendo el mismo color todos los individuos que pertenecen a un mismo cluster. En el centro de los clusters, del mismo color, se sit´ua un punto m´as grande que el resto, representando el baricentro o centroide de la agrupaci´on. Figura 5.10: Representaci´on bidimensional de las agrupaciones realizadas entorno a subespacios afines de dimensi´on 0. Cuando los subespacios afines son de dimensi´on unitaria ´estos quedan definidos por una recta (Figura 5.11). La recta, al igual que los centroides en el caso de subespacios de dimensi´on 0, est´a graficada en el mismo tono que los puntos que pertenecen al grupo que representa. Figura 5.11: Representaci´on bidimensional de las agrupaciones realizadas entorno a subespacios afines de dimensi´on 1. 50 5 IMPLEMENTACI´ ON Los puntos negros y con forma de aspa presentes en las representaciones gr´aficas se corresponden con las observaciones recortadas. Recordemos que dichas observaciones han sido tratadas como puntos at´ıpicos y no han influido en la determinaci´on de los subespacios afines. 5.1.3.2 3D Figura 5.12: Representaciones tridimensionales de las agrupaciones realizadas entorno a subespacios afines de dimensi´on 0. Si se pretende realizar el an´alisis cluster robusto en un conjunto de datos en R3, entonces se ofrecen dos posibilidades de representaci´on, una m´as sencilla y una m´as estilizada, cada una de las 51 5.1 Mejoras a˜nadidas cuales tiene sus ventajas y desventajas. Ambas representaciones son interactivas y permiten las rotaciones y el aumento sobre una zona concreta de la figura. Para ambos estilos de representaciones se ha utilizado el paquete contribuido rgl. Respecto a la representaci´on m´as sencilla, la principal ventaja que tiene es que conlleva menor tiempo para su generaci´on que la variante m´as estilizada y no se requieren especiales prestaciones en cuanto al hardware para su generaci´on. Como defecto se puede mencionar su aspecto menos llamativo. Figura 5.13: Representaciones tridimensionales de las agrupaciones realizadas entorno a subespacios afines de dimensi´on 1. 52 5 IMPLEMENTACI´ ON Respecto a la representaci´on m´as estilizada, su ´unica ventaja y no por ello menos importante es facilitar su visualizaci´on. Como h´andicaps se encuentran el trabajo nada desde˜nable que ha de realizar la tarjeta gr´afica de la computadora para generarlo e interactuar con ´el cuando los conjuntos de datos son de un tama˜no considerable. Las posibles dimensiones de los subespacios entorno a los que agrupar cuando los datos son tridimensionales son 1, 2 y 3. En los dos primeros casos, el modo de representar tanto las observaciones como los subespacios es id´entico al caso bidimensional, mostrando los centroides de los subespacios de dimensi´on 0 como puntos y los subespacios de dimensi´on 1 como rectas (Figuras 5.12 y 5.13). Figura 5.14: Representaciones tridimensionales de las agrupaciones realizadas entorno a subespacios afines de dimensi´on 2. 53 5.2 Pseudoc´odigo Cuando los uno o varios subespacios afines entorno a los que realizar clustering son de dimen- si´on 2, entonces dichas estructuras son representadas mediante la rotaci´on en los ejes ortogonales del subespacio de dos rectas (Figura 5.14). Dichas rectas iniciales se corresponden con las dos direcciones de m´axima varianza de los puntos que pertenecen al subespacio y que contengan al punto medio de dicha agrupaci´on. Una vez se disponen de esas dos rectas, ´estas se rotan como si se tratase de una peonza para generar una estructura con forma de plano que representa al subespacio af´ın de dimensi´on de 2. Igualmente, el color del plano es id´entico al color de las observaciones que pertenecen a su grupo. Al igual que en el caso bidimensional, y para ambos formatos de representaciones, los puntos recortados son representados en negro, esta vez sin forma de aspa. 5.2. Pseudoc´odigo Algorithm 1 Robust Linear Clustering Input:X={x1, . . . , xn}: conjunto de datos, d: vector con las dimensiones de los subespacios, α∈[0,1]: par´ametro de recorte, nstarts: n´umero de inicializaciones Output:clusters: grupo al que pertenece cada observaci´on, fObj: valor de la funci´on objetivo 1: function RobustLinearSubspacesInit(X,d,α) 2: k←tama˜no de d 3: n←n´umero de observaciones en X 4: p←dimensi´on de las observaciones en X 5: D ← matriz de dimensiones n×k 6: for i∈ {1, . . . , k}do 7: if d[i] == 0 then 8: c←observaci´on aleatoria de X 9: for l∈ {1, . . . , n}do 10: D[l, i]←qPj∈{1,...,p}(Xl,j −cj)2 11: end for 12: else 13: C←matriz de dimensiones (d[i] + 1) ×p 14: for s∈ {1, . . . , d [i]+1}do 15: Cs←observaci´on aleatoria de X 16: end for 17: ¯x←media de las observaciones en C 18: U ← matriz de dimensiones p×d[i] con los d[i] mayores autovectores como columnas 19: for l∈ {1, . . . , n}do 54 5 IMPLEMENTACI´ ON 20: D[l, i] = k{I− U(U)0}(Xl−¯x)k2 21: end for 22: end if 23: end for 24: clusters, minimos ←vectores n-dimensionales 25: for l∈ {1, . . . , n}do 26: minimos [l]←m´ıni=1,...,k D[l, i] 27: clusters [l]←ital que minimos [l] == m´ıni=1,...,k D[l, i] 28: end for 29: lim ← {n·αmayores valores en minimos} 30: clusters [minimos ∈lim]←1000 31: iter ←1 32: clustersant ←vector n-dimensional 33: while iter < 10 AND clusters 6=clustersant do 34: clustersant ←clusters 35: for i∈ {1, . . . , k}do 36: if d[i] == 0 then 37: ¯x←media de las observaciones donde clusters == i 38: for l∈ {1, . . . , n}do 39: D[l, i]←qPj∈{1,...,p}(Xl,j −cj)2 40: end for 41: else 42: ¯x←media de las observaciones donde clusters == i 43: U ← matriz de dimensiones p×d[i] con los d[i] mayores autovectores de la matriz de covarianzas de las observaciones donde clusters == icomo columnas 44: for l∈ {1, . . . , n}do 45: D[l, i] = k{I− U(U)0}(Xl−¯x)k2 46: end for 47: end if 48: end for 49: for l∈ {1, . . . , n}do 50: minimos [l]←m´ıni=1,...,k D[l, i] 51: clusters [l]←ital que minimos [l] == m´ıni=1,...,k D[l, i] 52: end for 53: lim ← {n·αmayores valores en minimos} 54: clusters [minimos ∈lim]←1000 55: iter ←iter + 1 56: end while 57: costfunction ←Pl∈{1,...,n},minimosl6=1000 minimos2 l [n[1−α]] 55 5.2 Pseudoc´odigo 58: return (clusters, costfunction) 59: end function 60: function RobustLinearSubspaces(X,d,α,clustersant,costfunction) 61: k←tama˜no de d 62: n←n´umero de observaciones en X 63: p←dimensi´on de las observaciones en X 64: D ← matriz de dimensiones n×k 65: clusters ←clustersant 66: clustersant ←vector n-dimensional 67: flag ←0 68: while flag < 100 AND clusters 6=clustersant do 69: clustersant ←clusters 70: for i∈ {1, . . . , k}do 71: if d[i] == 0 then 72: ¯x←media de las observaciones donde clusters == i 73: for l∈ {1, . . . , n}do 74: D[l, i]←qPj∈{1,...,p}(Xl,j −cj)2 75: end for 76: else 77: ¯x←media de las observaciones donde clusters == i 78: U ← matriz de dimensiones p×d[i] con los d[i] mayores autovectores de la matriz de covarianzas de las observaciones donde clusters == icomo columnas 79: for l∈ {1, . . . , n}do 80: D[l, i] = k{I− U(U)0}(Xl−¯x)k 81: end for 82: end if 83: end for 84: for l∈ {1, . . . , n}do 85: minimos [l]←m´ıni=1,...,k D[l, i] 86: clusters [l]←ital que minimos [l] == m´ıni=1,...,k D[l, i] 87: end for 88: lim ← {n·αmayores valores en minimos} 89: clusters [minimos ∈lim]←1000 90: flag ←flag + 1 91: end while 92: costfunction ←Pl∈{1,...,n},minimosl6=1000 minimos2 l [n[1−α]] 93: return (clusters, costfunction) 94: end function 56 5 IMPLEMENTACI´ ON 95: function MCD Mahalanobis(X,d,clusters) 96: k←tama˜no de d 97: n←n´umero de observaciones en X 98: p←dimensi´on de las observaciones en X 99: ∆←matriz de dimensiones n×k 100: for i∈ {1, . . . , k}do 101: ˆµ←media de las observaciones donde clusters == iestimada por MCD. 102: ˆ Σ←matriz de covarianzas de las observaciones donde clusters == i estimada por MCD. 103: for l∈ {1, . . . , n}do 104: ∆l,i ←q(Xl−ˆµ)0ˆ Σ−1(Xl−ˆµ) 105: end for 106: end for 107: for l∈ {1, . . . , n}do 108: minimos [l]←m´ıni=1,...,k D[l, i] 109: clustersnuevos [l]←ital que minimos [l] == m´ıni=1,...,k D[l, i] 110: end for 111: clustersnuevos = 1000 donde clusters == 1000 112: return (clustersnuevos) 113: end function 114: Escalar robustamente Xdividiendo entre su MAD. 115: gr ←matriz de dimensiones nstarts ×n 116: fObjs ←vector nstarts-dimensional 117: for init ∈ {1, . . . , nstarts}do 118: grinit, fObjsinit ←ROBUSTLINEARSUBSPACESINIT(X, d, α) 119: end for 120: gr10, fObjs10 ←agrupaciones y valores de las funciones objetivo de los 10 mejores inicios en base a fObjsinit 121: for init ∈ {1,...,10}do 122: grMininit, fObjsMininit ←ROBUSTLINEARSUBSPACES(X, d, α,gr10init, fObjs10init) 123: end for 124: agrupacionOpt, fObjOpt ←agrupaciones y valor de la funci´on objetivo del menor valor en fObjsMin 125: clusters ←MCD MAHALANOBIS(X,d,agrupacionOpt) 57 5.3 Computaci´on paralela Figura 5.15: Diagrama de flujo de la versi´on paralela del c´odigo. El esqueleto de la funci´on principal del programa que tiene que ser llamada para realizar el procedimiento es el siguiente: trimksubspaces<- function ( puntos ,d, alpha ,n_starts,plot=NULL , style = 1){ . . . return(list(funcion_coste = optimo , clusters = gruposRobustos )) } El primer paso que ha de hacerse es escalar robustamente las variables de los datos que se introduzcan al algoritmo, dividi´endolas entre su desviaci´on mediana absoluta. for(i in 1: dim( puntos ) [2]) { puntos [,i] <-puntos [,i]/mad( puntos [,i ]) } Una vez se disponen de los datos escalados y la funci´on a paralelizar adaptada (secci´on 5.2.1.1), hay que preparar el entorno de computaci´on paralela. Para ello, hay que seguir los siguientes pasos: 64 5 IMPLEMENTACI´ ON 1. Detectar el n´umero de n´ucleos disponibles en la m´aquina donde est´a ejecutando el c´odigo y crear un cluster con el n´umero de n´ucleos deseado. En este caso, todos los n´ucleos que haya menos 1 para no sobrecargar en exceso la m´aquina y no comprometer la ejecuci´on de otros procesos coexistentes en segundo plano. 2. Definir una secuencia {1, . . . , n}donde nes el n´umero de veces que se quiere ejecutar la funci´on deseada. En nuestro caso, nrepresenta el n´umero de inicializaciones aleatorias del algoritmo. 3. Utilizar una funci´on de la familia par*apply a la que se le indica el cluster de servidores creado, la secuencia, la funci´on a ejecutar y los argumentos que tiene la funci´on (si los tiene). No hay que olvidar que la funci´on tiene que estar adaptada para soportar el par´ametro dummy que representa la secuencia. En nuestro caso, la funci´on elegida es parLapply, que devolver´a los resultados de las nejecuciones de rlc ini en una lista. 4. Destruir el cluster de n´ucleos creado. if (!require(’parallel ’)) { install.packages (’parallel ’) library(’parallel ’) } cl <- makeCluster ( detectCores () -1) inicializacion<- seq (1 ,n_starts) results <- parLapply (cl , inicializacion , rlc _ini , puntos = puntos , d=d , alpha = alpha ) stopCluster (cl ) Una vez la funci´on parLapply ha ejecutado las veces indicadas por el usuario la funci´on rlc ini, se seleccionan las 10 inicializaciones m´as prometedoras minimizando el valor de la funci´on objetivo. Puede darse el caso de que no haya 10 inicializaciones aleatorias, ya sea o bien porque el usuario ha introducido como argumento n starts un n´umero menor que 10 o bien porque hay muchas inicializaciones que han generado una excepci´on durante su ejecuci´on por la configuraci´on de los subespacios, la naturaleza de los datos y el par´ametro de recorte. En ese caso, se seleccionan todas las disponibles. Una vez se han iterado completamente las mejores inicializaciones en la funci´on rlc se elige la que tiene un menor valor en la funci´on objetivo, y se realiza una ´ultima asignaci´on de grupos en base al MCD y la distancia de Mahalanobis robusta, tal y como se detalla en la secci´on 5.1.2. gruposRobustos<- allocateMahanalobis(dim ( puntos ) [1] ,d, grupos ) 65 5.4 Integraci´on de RyC++ El agrupamiento que devuelve la funci´on allocateMahanalobis es el agrupamiento final que devuelve el algoritmo. El valor de la funci´on de coste devuelto no es modificado por las posibles reasignaciones que realice la funci´on allocateMahanalobis. Por ´ultimo, si se ha indicado que se desea una salida gr´afica, se llama a una funci´on que se encarga de realizarla. if(!is.null(plot)){ if(plot) plottrim3d ( gruposRobustos ,d , dim ( puntos ) [2] , style ) } El cuerpo de la funci´on que realiza la funcionalidad gr´afica no es mostrada en esta memoria al considerarse que no es lo suficientemente compleja. 5.4. Integraci´on de RyC++ En ocasiones, Rno es suficientemente r´apido. Por mucho que se trate de optimizar todas y cada una de las l´ıneas de c´odigo disponibles en un programa, ´este es demasiado lento. Para conseguir una versi´on a´un m´as r´apida que la desarrollada mediante la paralelizaci´on del c´odigo escrito en Ry mejorar de manera dr´astica el rendimiento del procedimiento se desarrollar´an las partes m´as pesadas del algoritmo en otro lenguaje de programaci´on, concretamente C++. C++ es un lenguaje de programaci´on derivado del lenguaje C cuyo objetivo fue extender el lenguaje C a la orientaci´on a objetos. Una de las caracter´ısticas m´as atractivas de este lenguaje es su eficiencia computacional, motivo por el que se va a utilizar en este trabajo para implementar las porciones m´as costosas del procedimiento. Debido a que algunas partes del algoritmo no son computacionalmente pesadas y su implementaci´on es mucho m´as sencilla en Rque en C++, se crear´a un c´odigo en formato h´ıbrido, parte del algoritmo ser´a escrito en Ry parte en C++. En la Figura 5.16 se aprecia el pipeline de la ejecuci´on del procedimiento en funci´on de si la computaci´on es llevada a cabo por Ro por C++. Las partes m´as costosas del algoritmo, es decir la b´usqueda inicial de subespacios afines realizada en varias ocasiones y el ajuste completo de las b´usquedas iniciales m´as prometedoras, van a ser implementadas en C++. El resto del c´odigo, encargado de iniciar la ejecuci´on, seleccionar los inicios m´as prometedores y seleccionar finalmente el mejor resultado, as´ı como la etapa de Mahalanobis y la funcionalidad gr´afica, ser´a desarrollado en R. 66 5 IMPLEMENTACI´ ON Figura 5.16: Diagrama de flujo de la versi´on que integra RyC++. Para integrar c´odigo C++ en el proyecto de Rse hace uso del paquete Rcpp [11]. if (!require(’Rcpp ’)) { install.packages (’Rcpp ’) library(’Rcpp ’) } Rcpp ofrece una API (Application Programming Interface) que permite intercambiar objetos de R(incluyendo objetos de tipo S3 y S4) entre RyC++, es decir, act´ua como wrapper entre ambas plataformas. Permite definir funciones en C++ e intercalarlas dentro de un fichero con c´odigo R, o bien crear un fichero con c´odigo fuente en C++, compilarle y disponer de las funciones creadas para su uso como si se tratasen de funciones de R. Esta ´ultima opci´on es la escogida en este trabajo. Otro paquete de gran utilidad que se va a utilizar es RcppArmadillo [12], una biblioteca de ´algebra lineal con multitud de funciones implementadas en C++ necesarias para la implementaci´on del algoritmo. Por ejemplo, dispone de una funci´on para realizar un An´alisis en Componentes Principales, etapa necesaria en el procedimiento a desarrollar. El paquete ha sido desarrollado para proporcionar una sintaxis relativamente similar a la disponible en Matlab. 67 5.4 Integraci´on de RyC++ if (!require(’ RcppArmadillo ’)) { install.packages (’RcppArmadillo ’) library(’ RcppArmadillo ’) } 5.4.1. B´usqueda de subespacios afines inicial Como ya se explic´o en la secci´on anterior, la b´usqueda de subespacios afines inicial es la parte m´as costosa del procedimiento, y por ello es candidata a programarse en un lenguaje de programaci´on muy eficiente como lo es C++. El no determinismo que viene inducido por la selecci´on incial aleatoria de subespacios afines provoca que sea necesario un alto n´umero de incializaciones para garantizar una soluci´on razonalble. La funci´on, al estar escrita en C++, ha de residir en un fichero con extensi´on .cpp, para que posteriormente el compilador pueda compilarla correctamente. Recordar que esta parte se corresponde con la funci´on ROBUSTLINEARSUBSPACESINIT del pseudoc´odigo desarrollado. // [[ Rcpp :: export ]] Rcpp :: List rlc (Rcpp :: NumericMatrix puntos , Rcpp :: NumericVector d, double alpha ) { . . . return result; } En la cabecera de la funci´on se encuentra la orden // [[Rcpp::export]], la cu´al es necesaria para disponer de la funci´on en la sesi´on de Rdesde la que se enlace al fichero .cpp en el que se defina la funci´on. Si en ese mismo fichero se define una funci´on sin esa cabecera, entonces no ser´a accesible desde la sesi´on de R. Esto ´ultimo es ´util, y en este trabajo se ha realizado, en el caso de funciones que sirven de ayuda para realizar una tarea pero que ´unicamente tienen que ser llamadas desde otro punto del fichero C++. En este caso no es necesario encerrar la funci´on dentro de un bloque tryCatch para evitar que un error en una inicializaci´on afecte a las dem´as. El wrapper que proporciona Rcpp se encarga de gestionar dichos errores y no son propagados hasta el c´odigo Rque gestiona las llamadas a la funci´on escrita en C++. En el interior de la funci´on lo primero que ha de hacerse es declarar las variables que ser´an usadas dentro de la misma. int k = d. length (); 68 5 IMPLEMENTACI´ ON int n = puntos . nrow () ; int p = puntos . ncol () ; int n_atip = floor ( alpha *n); Rcpp :: NumericVector grupos (n); Rcpp :: NumericMatrix x(k,p); Rcpp :: NumericMatrix di(n,k); int numRand; Rcpp :: NumericVector minimos (n); Rcpp :: NumericVector minimos2 (n); double lim; Rcpp :: NumericVector temporal ; int asignados = 0; double funcion_coste = 0; int iter_ext = 0; Rcpp :: NumericVector gruposAnt (n); Posteriormente, se inicializan los subespacios afines mediante la extracci´on aleatoria de observaciones del conjunto de datos y se computa la distancia ortogonal de cada observaci´on a cada uno de los k clusters. for(int i = 0; i < k; i ++) { if(d[i] == 0){ numRand = rand () %n; x. row (i) = puntos . row ( numRand ); for(int j = 0; j < n; j ++) { di (j,i) = sqrt ( Rcpp :: sum ( pow ( puntos . row (j)-x. row (i) ,2) )); } }else { Rcpp :: NumericMatrix ci (d[i ]+1 , p); for(int j = 0; j < (d[i ]+1) ; j++) { numRand = rand () %n; ci .row (j) = puntos . row ( numRand ); } for(int j = 0; j < p; j ++) { x(i,j) = Rcpp :: mean (ci.column (j)); } Rcpp :: NumericMatrix U(d[i],p); Rcpp :: NumericMatrix eigen_vec (p ,p); eigen_vec = localpca (ci); for(int j = 0; j < d[i]; j ++) { U. row (j) = eigen_vec . column (j); } Rcpp :: NumericMatrix productoMat (p,p); arma :: mat arma_identidad = arma :: eye (p,p); arma :: mat arma_U = Rcpp ::as <arma ::mat >( U); arma :: mat arma_productoMat = arma_identidad - arma_U .t() * arma_U ; 69 5.4 Integraci´on de RyC++ arma :: mat arma_puntos = Rcpp ::as < arma :: mat >( puntos ); arma :: mat arma_x = Rcpp ::as <arma ::mat >( x); for(int j = 0; j < n; j ++) { di(j,i) = arma :: norm ( arma_productoMat *( arma_puntos .row (j)-arma_x .row(i)). t());; } } } Una vez realizada esta etapa, hay que asignar cada observaci´on al subespacio m´as cercano. Aquellas α·nobservaciones con distancias m´as altas son categorizadas como puntos at´ıpicos, asign´andolas al grupo 1000. for(int j = 0; j < n; j ++) { minimos [j] = Rcpp :: min (di. row (j)); minimos2 [j] = Rcpp :: min (di. row (j)); temporal = di.row (j); grupos [j] = std :: min_element ( temporal . begin () , temporal . end () ) - temporal . begin () ; } std :: sort ( minimos . begin () , minimos . end () ); if(n_atip >0) { lim = minimos [n - n_atip ]; }else { lim = 100000.0; } for(int j = 0; j < n; j ++) { if( minimos2 [j] > lim ){ grupos [j] = 1000; asignados += 1; } } for(int j = 0; j < n; j ++) { if(minimos2 [j] == lim && asignados < n_atip ){ grupos [j] = 1000; asignados += 1; } } La ´ultima etapa es realizar el paso de concentraci´on 9 veces, obteniendo una asignaci´on en clusters de las observaciones. Se mantiene la decisi´on se interrumpir el bucle y la ejecuci´on de la funci´on si en dos iteraciones consecutivas la asignaci´on de las observaciones a los diferentes grupos no var´ıa, significando que en las posteriores iteraciones tampoco variar´a. 70 5 IMPLEMENTACI´ ON while ( iter_ext < 10 && ! compara ( grupos , gruposAnt )){ funcion_coste = 0; asignados = 0; for(int ind = 0; ind < n; ind ++){ gruposAnt [ ind ] = grupos [ ind ]; } for(int i = 0; i < k; i ++) { x. row (i) = media ( puntos ,grupos ,i); if(d[i] == 0){ for(int j = 0; j < n; j ++) { di (j,i) = sqrt ( Rcpp :: sum ( pow ( puntos . row (j)-x. row (i) ,2) )); } }else { Rcpp :: NumericMatrix U(d[i],p); Rcpp :: NumericMatrix eigen_vec (p ,p); Rcpp :: NumericMatrix puntos_grupo = selecciona_puntos (puntos ,grupos ,i); eigen_vec = localpca ( puntos_grupo ); for(int j = 0; j < d[i]; j ++) { U. row (j) = eigen_vec . column (j); } Rcpp :: NumericMatrix productoMat (p,p); arma :: mat arma_identidad = arma :: eye (p,p); arma :: mat arma_U = Rcpp ::as <arma ::mat >( U); arma :: mat arma_productoMat = arma_identidad - arma_U .t() * arma_U ; arma :: mat arma_puntos = Rcpp ::as < arma :: mat >( puntos ); arma :: mat arma_x = Rcpp ::as <arma ::mat >( x); for(int j = 0; j < n; j ++) { di(j,i) = arma :: norm ( arma_productoMat *( arma_puntos .row (j)-arma_x .row(i) ).t());; } } } for(int j = 0; j < n; j ++) { minimos [j] = Rcpp :: min (di. row (j)); minimos2 [j] = Rcpp :: min (di. row (j)); temporal = di.row (j); grupos [j] = std :: min_element ( temporal . begin () , temporal . end () ) - temporal . begin () ; } std :: sort ( minimos . begin () , minimos . end () ); if(n_atip >0) { lim = minimos [n - n_atip ]; }else { lim = 100000.0; } for(int j = 0; j < n; j ++) { if( minimos2 [j] > lim ){ grupos [j] = 1000; asignados += 1; 71 5.4 Integraci´on de RyC++ } } for(int j = 0; j < n; j ++) { if( minimos2 [j] == lim && asignados < n_atip ){ grupos [j] = 1000; asignados += 1; } if( grupos [j] != 1000) { funcion_coste += minimos2 [j]; } } funcion_coste /= (n - n_atip ) iter_ext += 1; } La asignaci´on de las observaciones a los clusters y el valor de la funci´on de coste son introducidos en una lista para ser devueltos como resultado de la ejecuci´on de la funci´on. Rcpp :: List result = Rcpp :: List :: create (grupos , funcion_coste ); 5.4.2. Iteraci´on completa de las inicializaciones m´as prometedoras En la secci´on anterior se explic´o que, en la mayor´ıa de ocasiones, la iteraci´on completa de las inicializaciones m´as prometedoras no era un proceso lo suficientemente costoso desde el plano computacional como para paralelizarlo entre los distintos n´ucleos del procesador. Esto suced´ıa porque al crear el cluster de cores se necesitaba una cantidad razonablemente alta de tiempo que luego no se ve´ıa amortizado al paralelizar la ejecuci´on de la funci´on (en la mayor´ıa de los casos). Por el contrario, esta problem´atica no acontece al programar la funci´on en un script C++ que, posteriormente, se integra en Rmediante un wrapper. Las llamadas a las funciones escritas en C++ desde una sesi´on de Rson inmediatas, y es palmaria la mejor´ıa en tiempos de procesamiento que requiere una misma pieza de c´odigo escrita en C++ frente a su versi´on en R. Por esta raz´on, esta etapa es implementada en el script de C++. Al ser una funci´on que necesita ser llamada desde la sesi´on de R, hay que incluir en su cabecera la orden // [[Rcpp::export]]. // [[ Rcpp :: export ]] Rcpp :: List rlc2 (Rcpp :: NumericMatrix puntos , Rcpp :: NumericVector d, double alpha , Rcpp :: NumericVector grupos ) { . . . return result; } 72 5 IMPLEMENTACI´ ON La funci´on rlc2 itera la asignaci´on a subespacios afines inicial que es pasada en el argumento grupos hasta que o bien los grupos a los que se asigna cada observaci´on no var´ıan en 2 iteraciones consecutivas o bien se alcanza un n´umero m´aximo de iteraciones. El n´umero m´aximo de iteraciones es fijado en 100 al igual que en el caso de la versi´on paralelizada del c´odigo, permitiendo una posterior comparaci´on de tiempos. Al igual que en el caso de la funci´on anterior, lo primero consiste en declarar las variables que se van a utilizar dentro de la funci´on. int k = d. length (); int n = puntos . nrow () ; int p = puntos . ncol () ; int n_atip = floor ( alpha *n); Rcpp :: NumericMatrix x(k,p); Rcpp :: NumericMatrix di(n,k); Rcpp :: NumericVector minimos (n); Rcpp :: NumericVector minimos2 (n); double lim; Rcpp :: NumericVector temporal ; Rcpp :: NumericVector gruposAnt (n); int asignados = 0; double funcion_coste = 0; int flag = 0; Tras esta etapa previa de definici´on de variables, se procede a, en un bucle while, ajustar en profundidad la asignaci´on a subespacios afines pasada como argumento de la funci´on. while ( flag < 100 && ! compara ( grupos , gruposAnt )){ funcion_coste = 0; asignados = 0; for(int ind = 0; ind < n; ind ++){ gruposAnt [ ind ] = grupos [ ind ]; } for(int i = 0; i < k; i ++) { x. row (i) = media ( puntos ,grupos ,i); if(d[i] == 0){ for(int j = 0; j < n; j ++) { di (j,i) = sqrt ( Rcpp :: sum ( pow ( puntos . row (j)-x. row (i) ,2) )); } }else { Rcpp :: NumericMatrix U(d[i],p); Rcpp :: NumericMatrix eigen_vec (p ,p); Rcpp :: NumericMatrix puntos_grupo = selecciona_puntos (puntos ,grupos ,i); eigen_vec = localpca ( puntos_grupo ); for(int j = 0; j < d[i]; j ++) { U. row (j) = eigen_vec . column (j); } Rcpp :: NumericMatrix productoMat (p,p); 73 5.5 An´alisis comparativo de tiempos de ejecuci´on Figura 5.22: Regresi´on lineal a los tiempos medidos. 5.5.2. Funci´on trimkmeans vs funci´on trimksubspaces En esta segunda comparaci´on, se va a enfrentar la versi´on del algoritmo implementado m´as eficiente (R&C++) con una funci´on de un paquete disponible en el CRAN (Comprehensive R Archive Network) de R. Esta funci´on, llamada trimkmeans, pertenece al paquete trimcluster [13], desarrollado por C.Hennig, y permite ajustar el m´etodo de las k-medias recortadas a un conjunto de datos. Como la funci´on trimkmeans ´unicamente permite realizar agrupaciones entorno a centroides, los experimentos dise˜nados para realizar la comparaci´on de tiempos tienen la caracter´ıstica de que los subespacios afines entorno a los que agrupar son de dimensi´on 0. Cualquier otra combinaci´on de dimensiones de los subespacios soportada por la funci´on trimksubspaces (funci´on implementada en este trabajo) no ser´ıa soportada por la funci´on del paquete trimcluster. El n´umero de experimentos dise˜nados es 3. En el primero experimento, computacionalmente el menos costoso, se generan 200 observaciones en R2con el fin de agrupar entorno a dos clusters de dimensi´on 0. El segundo consta de 3000 puntos en R2, tratando de agrupar las observaciones entorno a 2 centroides. Por ´ultimo, constituyendo el experimento m´as exigente, en el tercero se generan 50000 observaciones en R3. El n´umero de clusters en este ´ultimo caso ser´a de 5. Al igual que anterior, cada par funci´on-experimento es realizado 5 veces, recogiendo los tiempos de computaci´on empleados. El promedio de dichos tiempos seg´un funci´on y experimento puede apreciarse en el Cuadro 5.3. 80 5 IMPLEMENTACI´ ON Versi´on Exp. 1 Exp. 2 Exp. 3 trimkmeans 0.687765 26.92146 3650.621 trimksubspaces 0.07548022 1.929447 57.65721 Cuadro 5.3: Tiempo promedio en segundos de las funciones sobre los experimentos. Figura 5.23: Diagrama de cajas del tiempo en escala logar´ıtmica frente a la funci´on (Exp. 1). Figura 5.24: Diagrama de cajas del tiempo en escala logar´ıtmica frente a la funci´on (Exp. 2). La diferencias en el tiempo de computaci´on entre ambas versiones son evidentes. En los 3 experimentos realizados, la funci´on construida en este trabajo es respectivamente 91, 13 y 63 veces m´as r´apida que la funci´on disponible en CRAN, todo ello a pesar de que la funci´on creada es m´as potente y permite agrupar entorno a subespacios de cualquier dimensi´on, incluso mezclando dimensiones. 81 5.5 An´alisis comparativo de tiempos de ejecuci´on En los diagramas de cajas de las Figuras 5.23, 5.24 y 5.25 en los que se representa el logaritmo de los 5 tiempos de ejecuci´on frente a la funci´on utilizada por cada experimento se observa como la funci´on implementada en este trabajo tiene un rendimiento muy superior a la que podemos encontrar en los repositorios de R. Al igual que se hizo con la comparaci´on anterior, se ajusta una recta de regresi´on a los tiempos logar´ıtmicos de cada experimento seg´un la versi´on (Figura 5.26). La evoluci´on de las rectas es bastante diferente, apreci´andose que la que menos pendiente tiene es la ajustada con los 15 tiempos disponibles de la funci´on desarrollada en ese trabajo. Figura 5.25: Diagrama de cajas del tiempo en escala logar´ıtmica frente a la funci´on (Exp. 3). Figura 5.26: Regresi´on lineal a los tiempos medidos. Esto se confirma viendo la pendiente de cada recta ajustada en el Cuadro 5.4. Si los experimentos realizados y la respuesta de las versiones fuesen un reflejo de la casu´ıstica total de posibles 82 5 IMPLEMENTACI´ ON combianciones a darse, esto significar´ıa que, con conjuntos de datos de cada vez m´as tama˜no, lafunci´on desarrollada en este trabajo ser´ıa la que se impondr´ıa en cuanto a que tardar´ıa menos tiempo en ajustar los subespacios. El escaso tama˜no muestral de 5 tomas de tiempo por cada par funci´on-experimento no permite sacar conclusiones fiables. Versi´on Intercept Pendiente trimkmeans -4.870 4.288 trimksubspaces -5.944 3.325 Cuadro 5.4: Intercept y pendiente de las rectas de regresi´on ajustadas. 83 6. Aplicaci´on real en la segmentaci´on de im´agenes El an´alisis cluster tiene una amplia variedad de aplicaciones en el panorama actual cient´ıfico. Para ilustrar un ejemplo de aplicaci´on real del algoritmo desarrollado se va a realizar segmentaci´on de im´agenes, en concreto segmentaci´on basada en colores. El proceso de dividir una imagen en m´ultiples regiones o segmentos es conocido como segmentaci´on. El objetivo de esta t´ecnica es dividir una imagen en regiones que, al contener menos informaci´on, pueden ser m´as f´aciles de analizar. T´ıpicamente, el proceso de segmentaci´on de im´agenes se utiliza para localizar objetos y fronteras contenidos en una imagen [14]. Ahondando m´as en el tema, esto se logra etiquetando cada p´ıxel de una imagen en una categor´ıa. Es decir, la segmentaci´on de im´agenes puede considerarse un problema de clasificaci´on donde cada p´ıxel ha de ser etiquetado en una categor´ıa distinta, donde todos los p´ıxeles que compartan etiqueta significa que tienen una caracter´ıstica com´un: pertenecen a un mismo cuerpo o a un objeto de las mismas caracter´ısticas t´ıpicamente (Figura 6.1). Figura 6.1: Segmentaci´on de im´agenes. La segmentaci´on de im´agenes es una tarea vital en una amplia gama de ramas del conocimiento. La segmentaci´on de im´agenes m´edicas [15] o en el mundo del industria [16] tiene mucho inter´es y multitud de modelos en el estado del arte son desarrollados para estas tareas. Para automatizar el proceso de segmentaci´on sobre im´agenes complejas, lo m´as com´un es aprovechar la informaci´on proporcionada por los colores presentes en la imagen. P´ıxeles cercanos en color son asociados a un mismo objeto y se clasifican en una misma categor´ıa. Es en esta tarea, en la b´usqueda de relaciones entre p´ıxeles en base a la informaci´on proporcionada por su color, es donde el an´alisis cluster y en concreto el procedimiento desarrollado pueden sacar ventaja sobre el resto de t´ecnicas, al ser un paradigma excelente en la b´usqueda de patrones que relacionen elementos entre s´ı (p´ıxeles en este caso). La imagen que se va a utilizar para ilustrar la aplicaci´on del an´alisis cluster robusto entorno a 84 6 APLICACI´ ON REAL EN LA SEGMENTACI´ ON DE IM´ AGENES subespacios afines es la de la Figura 6.2. Figura 6.2: Imagen base a la que se va a aplicar clustering. La imagen es de 800x533 p´ıxeles y se encuentra en 3 canales: R (Red), G (Green) y B (Blue). La imagen se aplana, formando un conjunto de datos de 426400 observaciones en R3. Los 426400 puntos representan los 800x533 p´ıxeles, y por cada p´ıxel se recogen 3 valores: su intensidad en el canal R, su intensidad en el canal G y su intensidad en el canal B (Figura 6.3). Los valores de intensidad oscilan entre 0 y 1, correspondi´endose el 1 con el valor m´as alto posible (t´ıpicamente 255) y 0 el valor m´as bajo posible. Figura 6.3: Conversi´on imagen-conjunto de datos. 85 6.1 Clustering de una imagen en 3 grupos. 6.1. Clustering de una imagen en 3 grupos. El primer experimento que se va a realizar es el de generar las im´agenes de segmentaci´on ´unicamente clasificando los p´ıxeles en 3 posibles valores. Para ello, se utiliza el procedimiento de an´alisis cluster robusto entorno a subespacios afines desarrollado en su versi´on en C++ al ser la m´as eficiente y la que soporta una mayor carga computacional con diferencia. La decisi´on de usar esta versi´on es debido al tama˜no masivo del conjunto de datos, con pr´acticamente medio mill´on de observaciones. Las configuraciones de dimensiones de los subespacios afines que se utilizan del algoritmo son todas los posibles para realizar agrupamientos en R3, es decir: 1. Tres grupos entorno a subespacios de dimensi´on 0 (k-medias con k= 3). 2. Dos grupos entorno a subespacios de dimensi´on 0 y un grupo entorno a un subespacio de dimensi´on 1. 3. Dos grupos entorno a subespacios de dimensi´on 1 y un grupo entorno a un subespacio de dimensi´on 0. 4. Tres grupos entorno a subespacios de dimensi´on 1. Todas las configuraciones son usadas sin recortar ninguna observaci´on (α= 0 %), para no perder informaci´on de la imagen, y con 50 inicios aleatorios. Figura 6.4: Im´agenes de segmentaci´on obtenidas tras realizar los agrupamientos entorno a dimensiones: 0-0-0, 0-0-1, 0-1-1 y 1-1-1 (de izquierda a derecha y de arriba a abajo). 86 6 APLICACI´ ON REAL EN LA SEGMENTACI´ ON DE IM´ AGENES En la Figura 6.4 se aprecia como las im´agenes de segmentaci´on tienden a formar agrupaciones de p´ıxeles bajo los mismos colores m´as alargadas conforme se utiliza el algoritmo con subespacios de dimensi´on unitaria en detrimento de subespacios de dimensi´on 0. Estas agrupaciones pueden aportar un valor a˜nadido diferencial a lo que ya aporta el algoritmo de las k-medias al permitir asociar entre s´ı partes m´as d´ebiles de las im´agenes que pueden ser objeto de inter´es en multitud de dominios de aplicaci´on. 6.2. Clustering de una imagen en 6 grupos. El segundo experimento contin´ua la l´ınea del primero, pero en vez de clasificar los p´ıxeles en ´unicamente 3 categor´ıas se van a clasificar en 6. Igualmente se utiliza el procedimiento de an´alisis cluster robusto entorno a subespacios afines desarrollado en su versi´on en C++. Las configuraciones de dimensiones de los subespacios afines que se utilizan del algoritmo son todas los posibles para realizar agrupamientos en R3, es decir: 1. Seis grupos entorno a subespacios de dimensi´on cero (k-medias con k= 6). 2. Cinco grupos entorno a subespacios de dimensi´on cero y un grupo entorno a un subespacio de dimensi´on uno. 3. Cuatro grupos entorno a subespacios de dimensi´on cero y dos grupo entorno a subespacios de dimensi´on uno. 4. Tres grupos entorno a subespacios de dimensi´on cero y tres grupos entorno a subespacios de dimensi´on uno. 5. Cuatro grupos entorno a subespacios de dimensi´on uno y dos grupos entorno a subespacios de dimensi´on cero. 6. Cinco grupos entorno a subespacios de dimensi´on uno y un grupo entorno a un subespacio de dimensi´on cero. 7. Seis grupos entorno a subespacios de dimensi´on uno. Al igual que en el experimento anterior, todas las configuraciones son usadas sin recortar ninguna observaci´on (α= 0 %), para no perder informaci´on de la imagen, y con 50 inicios aleatorios. Observando las im´agenes de segmentaci´on generadas en la Figura 6.5, igualmente se observa que los grupos de p´ıxeles que comparten categor´ıa (color) tienden a ser m´as alargados conforme hay m´as subespacios afines de dimensi´on 1 entorno a los que agrupar. Como era de esperar, las im´agenes de segmentaci´on generadas por agrupaciones entorno a 6 subespacios son de mayor 87 6.3 Clustering de una imagen con ruido calidad que las generadas por agrupaciones entorno a 3 subespacios, generalizando y resumiendo de manera m´as ´optima la informaci´on presente en la imagen original. Figura 6.5: Im´agenes de segmentaci´on obtenidas tras realizar los agrupamientos entorno a dimensiones: 0-0-0-0-0-0, 0-0-0-0-0-1, 0-0-0-0-1-1, 0-0-0-1-1-1, 0-0-1-1-1-1, 0-1-1-1-1-1 y 1-1-1-1-1-1 (de izquierda a derecha y de arriba a abajo). 6.3. Clustering de una imagen con ruido El tercer y ´ultimo experimento que se va a realizar tiene como objetivo mostrar la robustez estad´ıstica del procedimiento desarrollado. Para ello, a la imagen de la Figura 6.2 se le va a a˜nadir una cantidad de ruido aleatorio sobre ciertos p´ıxeles. Posteriormente, se proceder´a a realizar clustering sobre la imagen con ruido con el algoritmo desarrollado, tanto aplicando recortes (α > 0) como sin ellos (α= 0), con el fin de determinar si en la versi´on en la que se aplican recortes las observaciones recortadas son aquellos p´ıxeles contaminados. Para ello, primero se muestra la m´ascara de segmentaci´on generada a partir de la Figura 6.2 88 6 APLICACI´ ON REAL EN LA SEGMENTACI´ ON DE IM´ AGENES con la configuraci´on del algoritmo d= (0,0,0,0) y α= 0 (Figura 6.6). Figura 6.6: Imagen de segmentaci´on obtenida tras realizar los agrupamientos con d= (0,0,0,0). Posteriormente, se contaminan un 2 % de los p´ıxeles de la imagen original, dando lugar a la imagen presentada en la Figura 6.7. Los p´ıxeles con ruido contienen en los canales R y B un valor aleatorio proveniente de un distribuci´on uniforme continua U(0,0.05), mientras que el canal G contiene un valor aleatorio extra´ıdo de una distribuci´on N(0.9,0.0004). Figura 6.7: Imagen base con un 2 % de p´ıxeles contaminados. Las im´agenes de segmentaci´on generadas (Figura 6.8) utilizando configuraciones del procedimiento con α= 0 y α= 0.02 muestran la robustez del m´etodo. En la primera configuraci´on, los p´ıxeles contaminados influyen en la determinaci´on de los subespacios afines y por tanto en 89