Análisis topológico de datos

Facultad de Matem´aticas Departamento de ´ Algebra An´alisis topol´ogico de datos Javier Perera Lago Doble Grado en Matem´aticas y Estad´ıstica A.A. 2019-2020 2 Abstract Topological data analysis is a branch of computational topology which uses algebra to obtain topological features from a data set. It has many applications in computer vision, shape description, time series analysis, biomedicine, drug design... The first step to learn topological information from data is to build a filtration of nested simplicial complexes estimating the global structure of the data set at different scales of precision. The second step is to study the evolution of the simplicial homology along this filtration, using a tool from algebraic topology called persistent homology. When the filtration is indexed by a single parameter, we can describe discrete and complete invariants for persistent homology with coefficients over a field, such as barcodes, persistent diagrams or persistent landscapes. Unfortunately, we don’t have a simple classification for persistent homology when the filtration is indexed by two or more parameters. The problem of classifying multiparameter persistent homology is being widely studied nowadays, and we present some proposals for it. One of the big problems of large sets is the presence of noisy, wrong or incomplete data. Although persistent homology is stable under small perturbation, we still need to give a statistical frame in order to separate insignificant features from topological signal and build some hypothesis tests applied to topological features. Another usual problem of large sets is the high dimensionality. Even if we have captured the main topological caracteristics of the data, we aren’t usually able to fully understand its structure. In order to reduce dimensionality and obtain a better visualization we present the Mapper algorithm, which groups the data using a filter function and summarizes the set as a graph or a low dimensional simplicial complex. 3 4 ´ Indice general 1. Complejos simpliciales 7 1.1. Primerasdefiniciones................................. 7 1.2. Recubrimientoynervio................................ 8 1.3. Complejo de ˇ Cech................................... 8 1.4. Complejo alpha .................................... 9 1.5. Complejo de Vietoris-Rips . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.6. Complejos witness .................................. 11 1.7. Filtraciones y multifiltraciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.8. Algoritmo Mapper .................................. 15 2. Homolog´ıa simplicial 17 2.1. Gruposdecadenas .................................. 17 2.2. Operadorborde.................................... 18 2.3. Gruposdehomolog´ıa................................. 18 2.4. Computaci´on de la homolog´ıa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3. Homolog´ıa persistente 25 3.1. Complejos y m´odulos de persistencia . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2. M´odulos de homolog´ıa persistente . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3. Computaci´on de la homolog´ıa persistente . . . . . . . . . . . . . . . . . . . . . . 28 3.4. C´odigos de barras y curvas de Betti . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.5. Diagramas de persistencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4. Homolog´ıa multipersistente 35 4.1. Computaci´on de la homolog´ıa multipersistente . . . . . . . . . . . . . . . . . . . 36 4.2. Invariantedelrango.................................. 37 4.3. Diagrama de persistencia multiparam´etrico . . . . . . . . . . . . . . . . . . . . . 38 5. Inferencia estad´ıstica aplicada a la homolog´ıa persistente 41 5.1. Paisajes de persistencia y siluetas . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.2. Teor´ıa de probabilidad sobre paisajes y siluetas . . . . . . . . . . . . . . . . . . 44 5.3. Bandas de confianza basadas en bootstrap ...................... 46 5.4. Bandas bootstrap de confianza para paisajes y siluetas . . . . . . . . . . . . . . . 46 5.5. Conjuntosdeconfianza................................ 48 5.6. An´alisis de m´axima persistencia . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 6. Aplicaciones del an´alisis topol´ogico de datos 51 6.1. Machine learning basado en homolog´ıa persistente . . . . . . . . . . . . . . . . . 51 6.2. An´alisisdeformas................................... 52 6.3. An´alisis de series temporales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5 6´ INDICE GENERAL 6.4. Algoritmo Mapper enbiomedicina.......................... 54 6.5. Software para el an´alisis topol´ogico de datos . . . . . . . . . . . . . . . . . . . . 55 A. Funciones de R 63 Cap´ıtulo 1 Complejos simpliciales El punto de partida en el an´alisis topol´ogico de datos es un conjunto de puntos S= {x1,··· , xN} ∈ Rd, aunque tambi´en pueden estar en un espacio m´etrico general. Siempre se asume uno de los siguientes escenarios: o bien el conjunto Sviene muestreado de una variedad topol´ogica no observable M, o bien es una muestra aleatoria simple de una distribuci´on de probabilidad Pcuyo soporte es M. El primer escenario asume que nuestros datos son correctos y est´an medidos con precisi´on pero podr´ıan ser m´as numerosos, mientras que el segundo asume que los datos tienen una componente aleatoria y pueden estar sujetos a error. En cualquier caso, nuestro objetivo es deducir las propiedades topol´ogicas de Ma partir del conjunto S. El primer paso para ello ser´a construir un complejo simplicial sobre los puntos de Sque reconstruya M de forma aproximada. 1.1. Primeras definiciones Definici´on 1.1. Un complejo simplicial abstracto es un conjunto no vac´ıo Ktal que todos sus elementos son conjuntos finitos y si σ∈K, y adem´as τ⊂σ, entonces se verifica que τ∈K. En particular, ∅ ∈ K. Cada elemento de un complejo simplicial se llama s´ımplice. En otras palabras, es una generalizaci´on del concepto de grafo donde se permiten relaciones m´as complejas que las aristas, es decir, con tres elementos o m´as [16]. Definici´on 1.2. Un s´ımplice con n+ 1 elementos tiene dimensi´on ny se llama un n-s´ımplice. En particular ∅tiene dimensi´on −1. La dimensi´on de un complejo simplicial es el m´aximo de las dimensiones de sus s´ımplices. Dados dos s´ımplices σ, τ ∈Ktales que τ⊂σ, se dir´a que τes una cara de σy que σes una cocara de τ. Un s´ımplice se dir´a maximal si no tiene cocaras propias en K. Definici´on 1.3. Un subcomplejo de Kes un subconjunto L⊂Kque satisface la definici´on de complejo simplicial. Un tipo especial de subcomplejo es el n-esqueleto de K, formado por todos los s´ımplices de Kcon dimensi´on no menor. Existe una versi´on geom´etrica de la definici´on de complejo simplicial, que nos har´a m´as f´acil las construcciones sobre S[49]. Definici´on 1.4. Un s´ımplice geom´etrico σde dimensi´on nes la envoltura convexa de {a0,··· , an} ⊂ Rdaf´ınmente independientes. La envoltura convexa de un subconjunto de {a0,··· , an}es una cara de σ. Un complejo simplicial geom´etrico Kes una colecci´on de s´ımplices geom´etricos en Rdtal que cualquier cara de σ∈Kest´a en Ky la intersecci´on de dos s´ımplices de Kes o bien vac´ıa o bien una cara com´un. La dimensi´on de un complejo simplicial geom´etrico es el m´aximo de las dimensiones de sus s´ımplices. 7 8CAP´ ITULO 1. COMPLEJOS SIMPLICIALES En esta versi´on geom´etrica, un s´ımplice de dimensi´on 0 es un v´ertice, un s´ımplice de dimensi´on 1 es un segmento o arista, uno de dimensi´on 2 es un tri´angulo, uno de dimensi´on 3 es un tetraedro y as´ı sucesivamente. Esta nomenclatura geom´etrica se aplica tambi´en en los complejos abstractos. 1.2. Recubrimiento y nervio Para construir algunos de los complejos simpliciales que definiremos m´as tarde necesitamos los conceptos de recubrimiento y nervio [16]. Definici´on 1.5. Dado un conjunto finito Sde puntos eucl´ıdeos, un recubrimiento de Ses una familia de conjuntos no vac´ıos U={Ui}i∈Ital que S⊂Si∈IUi. Si todos los Uison abiertos (cerrados) no vac´ıos, se dice que Ues un recubrimiento abierto (cerrado) de S. Definici´on 1.6. Dado un recubrimiento U={Ui}i∈I, su nervio es un complejo simplicial abstracto Nconstruido sobre el conjunto de ´ındices I, de forma que ∅ ∈ Ny cada subconjunto no vac´ıo de ´ındices J⊂Ipertenece a Nsi y solo si Tj∈JUj6=∅. En particular, el s´ımplice unitario ipertenece a Npara todo i∈I. Dicho de otra forma, el nervio asigna un v´ertice a cada conjunto del recubrimiento y despu´es crea una arista por cada pareja de conjuntos que se corten, un tri´angulo por cada terna de conjuntos que se corten y as´ı sucesivamente. Un motivo por el que esta construcci´on resulta tan ´util es la existencia del siguiente teorema, que garantiza una relaci´on de homotop´ıa entre el recubrimiento y su nervio [45]. Teorema 1.1 (del nervio).Sea Mun espacio topol´ogico inmerso en Rdy sea U={Ui}i∈Iun recubrimiento abierto numerable para M. Supongamos que para cada ∅ 6=J⊂Ise tiene que Tj∈JUjes contr´actil o vac´ıo. Entonces el nervio Nde Ues homot´opicamente equivalente a M. Recordemos que dos espacios topol´ogicos XeYse dicen homot´opicamente equivalentes si existen dos aplicaciones continuas f:X→Y,g:Y→Xtales que g◦fes homot´opica (se puede deformar continuamente) a IdXyf◦ges homot´opica a IdY. 1.3. Complejo de ˇ Cech El complejo de ˇ Cech es el complejo simplicial m´as elemental que se puede construir sobre S, y el que describe de una forma m´as natural c´omo este conjunto recubre M[16]. Definici´on 1.7. Dados Sy un cierto ε > 0, sea el recubrimiento abierto Uε={B(x, ε)|x∈S}. Definimos el complejo de ˇ Cech C(S, ε)como el nervio de Uε. Definimos tambi´en C(S, 0) = ∅ Los v´ertices de C(S, ε) son los puntos de S, y dado un subconjunto de ellos, generan un s´ımplice si la intersecci´on de sus respectivas bolas es no vac´ıa. Como estas bolas son convexas y contr´actiles, el complejo de ˇ Cech es homot´opicamente equivalente a Uεseg´un el teorema del nervio. Es claro que si 0 < δ < ε, entonces C(S, δ)⊂C(S, ε). Adem´as para un cierto M > 0 suficientemente grande, el complejo C(S, M) ser´a un (N−1)-s´ımplice. Esta construcci´on es muy ´util desde el punto de vista te´orico gracias al siguiente resultado [16]: 1.4. COMPLEJO ALPHA 9 (a) Recubrimiento Uεde Spara ε= 0.6 (b) Complejo de ˇ Cech C(S, 0.6) Figura 1.1: Construcci´on de un complejo de ˇ Cech a partir de su recubrimiento. N´otese que en la parte inferior derecha hay tres puntos alineados, que definen un tri´angulo en el complejo aunque no quede representado visualmente. Teorema 1.2. Sea Muna variedad Riemanniana compacta inmersa en Rd. Entonces existe un cierto e > 0tal que para todo ε≤ese tiene que C(M, ε)es homot´opicamente equivalente a M. Es m´as, existe un subconjunto finito V⊂Mtal que C(V, ε)tambi´en es homot´opicamente equivalente a M. Es decir, que bajo la hip´otesis de que los puntos de Sest´en bien distribuidos podemos conseguir un complejo homot´opicamente equivalente a M. No obstante, su tama˜no final desaconseja usarlo en la pr´actica. En el caso extremo en que se obtiene un (N−1)-s´ımplice, se obtienen un total de 2N s´ımplices de dimensiones entre −1 y N−1, que en general es mucho mayor que la dimensi´on del espacio donde se encuentran nuestros puntos. Por otra parte, su computaci´on no solo es compleja sino que adem´as no es exacta sino aproximada. Comprobar si una pareja de puntos define una arista es sencillo pues basta ver si su distancia es menor que 2ε, lo cual es f´acil de verificar usando la matriz de distancias de S. En cambio si queremos comprobar si tres o m´as v´ertices definen un s´ımplice, debemos probar que la intersecci´on de sus bolas con radio εtenga al menos un punto en com´un. Esto es equivalente a probar que existe una bola con radio εtal que contenga a dichos v´ertices. Lo cual nos obliga a calcular el problema de la m´ınima bola envolvente (o MEB por sus siglas en ingl´es de ”Minimum Enclosing Ball”) para cada posible s´ımplice. Los algoritmos que resuelven este problema en general no calculan el radio exacto de esta bola, sino una sobreestimaci´on suya [71]. Por tanto, corremos el riesgo de obtener un complejo mayor que el real. A continuaci´on describiremos nuevos complejos que traten de mantener las buenas propiedades topol´ogicas del complejo de ˇ Cech y a la vez resuelvan sus dificultades computacionales. 1.4. Complejo alpha Esta estructura es una variante del complejo de ˇ Cech donde limitamos su dimensi´on a la del espacio ambiente. En primer lugar tomamos el conjunto de puntos Sy construimos su diagrama de Voronoi. Definici´on 1.8. La regi´on de Voronoi de un punto x∈Ses el conjunto R(x) := {y∈ Rd|d(x, y)≤d(x0, y),∀x0∈S\ {x}}. Estas regiones definen un recubrimiento cerrado de S, y su nervio es el complejo de Delaunay D(S). Las regiones de Voronoi nos sirven para limitar la expansi´on de las bolas de radio ε, dando lugar a un subcomplejo del complejo de ˇ Cech llamado el complejo alpha [38]. 16 CAP´ ITULO 1. COMPLEJOS SIMPLICIALES Ep(y) = Px∈S||y−x||p N1/p .(1.3) En general los complejos simpliciales que se obtengan con este algoritmo no se usar´an para el an´alisis topol´ogico que se desarrollar´a en los pr´oximos cap´ıtulos, sino para el resumen y representaci´on visual de conjuntos de datos muy numerosos [16]. No obstante, es ampliamente usado por un gran n´umero de autores para ejemplos de aplicaci´on pr´actica [51] y ser´a en la ´ultima secci´on cuando veamos algunos de ellos. Cap´ıtulo 2 Homolog´ıa simplicial El grupo fundamental π1(X), que estudia los lazos de un espacio topol´ogico X, tiene una generalizaci´on en los grupos de homotop´ıa πn(X), que se definen en t´erminos de aplicaciones entre cubos n-dimensionales InyXy homotop´ıas entre estas aplicaciones. Aunque tienen un gran inter´es te´orico, en general estos grupos son dif´ıciles de computar. En esta secci´on presentamos una alternativa: los grupos de homolog´ıa. En particular estudiaremos la homolog´ıa de un complejo simplicial Kdado. Estos grupos son m´as f´acilmente computables que los grupos de homotop´ıa, pero tienen una definici´on menos natural y requieren resolver ciertas cuestiones t´ecnicas para su c´alculo [40]. Dado un complejo simplicial K, supondremos a lo largo de toda la secci´on que sus v´ertices est´an ordenados, y se escribir´an como x1, x2,··· , xN. 2.1. Grupos de cadenas Siguiendo la definici´on 1.1 de complejo simplicial, un s´ımplice σviene determinado por un conjunto de v´ertices {xi0, xi1,··· , xip} ∈ K, siendo irrelevante el orden en que se presenten. Por ejemplo, los conjuntos {x1, x2}y{x2, x1}definen el mismo s´ımplice. Para definir la homolog´ıa simplicial supondremos que los s´ımplices est´an ordenados, y para remarcar que el orden importa cambiamos las llaves por corchetes, escribiendo σ= [xi0, xi1,··· , xip]. Entonces, [x1, x2]y[x2, x1] son diferentes ordenaciones del mismo s´ımplice. Definici´on 2.1. Decimos que un s´ımplice ordenado [xi0, xi1,··· , xip]est´a bien ordenado si i0< i1<··· < ip, es decir, si sus v´ertices siguen la ordenaci´on que dimos al principio en K. Por otra parte, diremos que un s´ımplice ordenado est´a orientado positivamente si existe una permutaci´on par de sus v´ertices que lo convierta en un s´ımplice bien ordenado. Si no existe esta permutaci´on par, diremos que el s´ımplice ordenado est´a orientado negativamente. Siguiendo con el ejemplo anterior, [x1, x2] tendr´ıa orientaci´on positiva porque est´a bien ordenado mientras que [x2, x1] tendr´ıa orientaci´on negativa. Una vez definidos estos t´erminos, podemos definir los grupos de cadenas, que ser´an nuestro primer paso en la definici´on de la homolog´ıa [49]. Definici´on 2.2. Dado un complejo simplicial K, su grupo n-´esimo de cadenas Cn(K)es un grupo abeliano libre que tiene por generadores a sus n-s´ımplices bien ordenados, donde un reordenamiento de sus v´ertices es equivalente a un cambio de signo si la permutaci´on es impar y no afecta nada si la permutaci´on es par. En nuestro ejemplo, [x1, x2] = −[x2, x1]. Si un complejo simplicial no tiene n-s´ımplices, entonces no hay generadores y Cn(K) = {0}. A partir de K, podemos crear el grupo Cn(K) para todo n∈Z. 17 18 CAP´ ITULO 2. HOMOLOG´ IA SIMPLICIAL Definici´on 2.3. Un elemento c∈Cn(k)es una n-cadena, y se puede representar como c= Pici·σi, donde los σison n-s´ımplices bien ordenados de K, y los cison coeficientes en Z. En las pr´oximas secciones veremos c´omo se relacionan los grupos de cadenas entre s´ı dando lugar al concepto de homolog´ıa simplicial. 2.2. Operador borde Una vez definidos los grupos de cadenas, pasamos a estudiar c´omo se relacionan entre ellos mediante los operadores borde, que asignan a cada s´ımplice orientado sus caras orientadas [40]. A partir de ahora omitiremos el complejo simplicial Kcuando nombremos a sus grupos de cadenas. Definici´on 2.4. El operador borde n-´esimo es un homomorfismo de grupos ∂n:Cn→Cn−1 que opera sobre sus s´ımplices de la siguiente forma: ∂n[v0,··· , vn] = n X i=0 (−1)i[v0,··· ,ˆvi,··· , vn],(2.1) donde ˆviindica que el v´ertice vise ha eliminado del s´ımplice. ∂nse extiende al resto de cadenas compuestas por linealidad. Definici´on 2.5. Dada una n-cadena c, al elemento ∂n(c)se le llama borde de c. Una n-cadena con borde nulo se llama n-ciclo. Los operadores borde conectan los diferentes grupos de cadenas entre s´ı, dando lugar a un complejo de cadenas [40]. Definici´on 2.6. Dado un complejo simplicial K, su complejo de cadenas C∗es un conjunto {{Cn}n∈N,{∂n}n∈N}de grupos de cadenas y operadores borde formando el siguiente diagrama: C∗:··· → Cn+1 ∂n+1 −→ Cn ∂n −→ Cn−1→ ··· (2.2) A continuaci´on probamos una caracter´ıstica importante de los complejos de cadenas [49]. Proposici´on 2.1. ∂p−1◦∂p≡0∀p∈Z. Corolario 2.1. im ∂p⊂ker ∂p−1∀p= 1,··· , m + 1. Esta proposici´on implica que im ∂p⊂ker ∂p−1. Dicho en lenguaje natural, el borde de cualquier n-cadena es un ciclo de Cn−1. 2.3. Grupos de homolog´ıa El corolario 2.1 afirma que im ∂p⊂ker ∂p−1. En particular, im ∂pes un subgrupo de ker ∂p−1 y podemos dar la siguiente definici´on [40]: Definici´on 2.7. Sea un complejo simplicial Ky sea C∗su complejo de cadenas. Definimos el grupo de homolog´ıa n-´esimo de K(o de C∗) como el cociente: Hn=ker ∂n im ∂n+1 (2.3) 2.4. COMPUTACI ´ ON DE LA HOMOLOG´ IA 19 Recordemos que ker ∂nes el grupo de todos los n-ciclos de K, mientras que im ∂n+1 es el grupo de los n-ciclos que son borde de alguna (n+ 1)-cadena de K. Al hacer el cociente, Hn solo conserva los ciclos que no bordean ninguna cadena de dimensi´on mayor, es decir, conserva los ciclos que representan “agujeros vac´ıos” en K. Como los grupos de homolog´ıa son grupos abelianos con un n´umero finito de generadores, en particular son Z-m´odulos finitamente generados y les podemos aplicar el siguiente teorema [32] para dar una clasificaci´on suya. Teorema 2.1 (teorema fundamental de los R-m´odulos finitamente generados).Sea Run dominio de ideales principales (como Z). Todo R-m´odulo finitamente generado se descompone de forma ´unica salvo orden como la suma directa: β M i=1 R⊕ m M j=1 R/tjR, (2.4) donde βes un entero no negativo y los tjson elementos no nulos y no unidades de Rtales que tj|tj+1. El primer elemento de esta suma directa se llama subm´odulo libre y se caracteriza por su rango o n´umero de Betti β. Cuando hablemos de un grupo de homolog´ıa Hn, el n´umero de Betti se escribir´a βn. El segundo elemento de la suma se llama subm´odulo de torsi´on y se caracteriza por sus coeficientes de torsi´on t1,··· , tm(salvo producto por una unidad). En conjunto el R-m´odulo queda parametrizado por los m+ 1 elementos β, t1,··· , tm. En complejos simpliciales de dimensi´on 3 sin torsi´on en su grupo de homolog´ıa, los n´umeros de Betti tienen un significado intuitivo: β0es el n´umero de componentes conexas, β1es el n´umero de “t´uneles” y β2representa el n´umero de “vac´ıos” o “burbujas” [4]. Aunque nuestra definici´on de homolog´ıa se basa en grupos abelianos o Z-m´odulos, se puede extender poniendo coeficientes en cualquier dominio de ideales principales en lugar de Z. Por este motivo, de aqu´ı en adelante no hablaremos de grupos de cadenas o de homolog´ıa sino de m´odulos. En el caso particular en que usemos coeficientes en un cuerpo k, el subm´odulo de torsi´on desaparece, y Hnse reduce a un espacio vectorial de dimensi´on igual a su n´umero de Betti βn. Una propiedad muy importante de los m´odulos de homolog´ıa simplicial es que se mantienen invariantes bajo la acci´on de la homotop´ıa [40]. Proposici´on 2.2. Sean K1yK2dos complejos simpliciales geom´etricos homot´opicamente equivalentes. Entonces los R-m´odulos de homolog´ıa Hn(K1)yHn(K2)son isomorfos para todo ny para todo anillo R. Entonces dado Sy un εlos complejos C(S, ε) y A(S, ε) no son en general iguales pero al ser homot´opicamente equivalentes sus m´odulos de homolog´ıa s´ı lo son. Como veremos en la pr´oxima secci´on, la elecci´on de un anillo u otro puede simplificar o complicar la computaci´on de las distintas homolog´ıas para un complejo K, pero modifica su estructura as´ı que debemos ser cautelosos si no queremos perder informaci´on relevante. 2.4. Computaci´on de la homolog´ıa Una vez definida la homolog´ıa simplicial, procedemos a su computaci´on efectiva. Para ello debemos calcular el n´ucleo y la imagen de cada operador borde ∂n. Como son operadores lineales entre m´odulos, se pueden representar a trav´es de una matriz Mnen t´erminos de ciertas 20 CAP´ ITULO 2. HOMOLOG´ IA SIMPLICIAL (a) Complejo simplicial K         bc bd be cd ce a00000 b−1−1−1 0 0 c100−1−1 d01010 e00101         (b) M1         bcd bc 1 bd −1 be 0 cd 1 ce 0         (c) M2 Figura 2.1: Usaremos este complejo Kpara ilustrar los c´alculos descritos en esta secci´on. Partimos de las matrices M1yM2respecto de las bases naturales de C0,C1yC2. Por comodidad cambiamos la notaci´on con corchetes por la yuxtaposici´on para remarcar el orden de los v´ertices. bases. Nuestra elecci´on natural ser´a la de tomar los n-s´ımplices orientados de Kcomo base de Cn. Calcular el n´ucleo y la imagen de ∂nequivale a calcular el espacio nulo y el espacio de columnas de la matriz Mnrespectivamente. La herramienta principal que usaremos para calcular la homolog´ıa es la forma normal de Smith [5]. Teorema 2.2. Sea Auna matriz con entradas en el anillo R. Existen dos matrices PyQ productos de de matrices elementales tales que A0=Q−1AP, donde A0tiene la forma:      a1...0 ar 0 0      ,(2.5) ya1|a2|···|ar. En estas condiciones, se dice que A0es la forma normal de Smith de A. La forma m´as sencilla de leer la estructura de Hnse basa en expresar la matriz del homomorfismo ∂n+1 con respecto a la base natural de Cn+1 y una base de ker ∂n. Este cambio de base es trivial gracias al siguiente resultado: Lema 2.1 (del cambio de base).Sea Mk+1 la matriz de ∂k+1 con respecto a las bases naturales de Cn+1 yCny sea M0 kla forma normal de Smith de Mk. Se puede obtener una matriz de ∂k+1 con respecto a la base natural de Ck+1 y una base de ker ∂ksimplemente eliminando de Mk+1 las filas que corresponden con columnas pivotales (es decir, no nulas) en M0 k. Demostraci´on. Dado que ∂k◦∂k+1 ≡0, entonces Mk·Mk+1 = 0. Para transformar Mken M0 ksolo se usan operaciones elementales por filas y columnas de los tipos 1 y 3. Para que el producto de las dos matrices siga siendo nulo durante la reducci´on de Mk, una operaci´on en las columnas de Mkdebe ir acompa˜nada de una operaci´on por filas en Mk+1. Las operaciones por filas en Mkno afectan a la matriz Mk+1. Una transposici´on de columnas en Mkse corresponde con una transposici´on de sus filas respectivas en Mk+1. Si reemplazamos (columna i) por (columna i) +q(columna j) en Mk, debemos reemplazar (fila j) por (fila j) −q(fila i) en Mk+1. Tras acabar la reducci´on, se ve que si una columna de Mkha sido anulada, su correspondiente fila de Mk+1 no ha sido alterada en ning´un momento. Entonces tenemos: 2.4. COMPUTACI ´ ON DE LA HOMOLOG´ IA 21 0 = Mk·Mk+1 =M0 k·˜ Mk+1 = (T|0) ·A B=T·A+ 0 ·B=T·A, donde Tes una submatriz de ˜ Mkcon sus columnas pivotales, y 0 es una submatriz nula. Las filas de Bson aquellas de Mk+1 que no han sufrido alteraci´on en el proceso, y se corresponden con una base de ker ∂k. Dado que T·A= 0, y todas las columnas de Tson linealmente independientes, A= 0 y podemos eliminar sus filas de Mk+1. Este lema y su prueba se basan en el art´ıculo [75], que est´a referido a la homolog´ıa persistente. En su versi´on original M0 kno es la forma normal de Smith sino la forma escalonada por columnas de Mk. Nuestra prueba es esencialmente la misma, pero tambi´en considera el efecto de las operaciones por filas en Mk. Por tanto hemos conseguido una mejora en el tiempo de computaci´on reduciendo en buena parte el tama˜no de las matrices implicadas. Una vez que la matriz de ∂n+1 est´a escrita en funci´on de las bases de Cn+1 y ker ∂n, la llevamos a su forma normal de Smith mediante operaciones elementales por filas y columnas, como en (2.5). Sean c1,··· , cslos ciclos correspondientes a las filas de esta matriz, o equivalentemente la base de ker ∂n. Entonces, el R-m´odulo n-´esimo de homolog´ıa es: Hn=ker ∂n im ∂n+1 ∼ =hc1,··· , cr,··· , csi ha1c1,··· , arcri∼ = s−r M i=1 R⊕ r M j=1 R/ajR(2.6) Si ajes unidad de Rpara alg´un j, el cociente R/ajRse anula y no se incluye en la suma directa. Entonces, tenemos que βn=s−ry los elementos torsionarios de Hnson los ajno unidades [32]. Con estas herramientas, podemos dar un algoritmo que nos ayude a calcular todos los grupos de homolog´ıa en cadena [75]. Algoritmo 2.1 (c´alculo en cadena de todos los m´odulos de homolog´ıa).Sea Kun complejo simplicial de dimensi´on m, y sean M1,··· , Mmlas matrices de los operadores borde ∂1,··· , ∂m respecto de las bases naturales de C0, C1,··· , Cm. Reducimos M1a su forma normal de Smith. Extraemos la estructura de H0como en la ecuaci´on (2.6), recordando que ker ∂0=C0. Tomamos M2y eliminamos las filas que se correspondan con columnas pivotales de M1 en virtud del lema 2.1. Reducimos la matriz M2truncada a su forma normal de Smith y extraemos la estructura de H1como en (2.6). Repetimos los dos pasos anteriores con las matrices M3,··· , Mm, calculando los m´odulos H2,··· , Hm−1. Dado que im ∂m+1 = 0, tenemos que Hm= ker ∂m. El resto de m´odulos de homolog´ıa son triviales. Si solo quisi´eramos calcular la estructura de Hnpara un ´unico nsin tener en cuenta qu´e ciclos generan ker ∂ny im ∂n+1, bastar´ıa con llevar Mn+1 a su forma normal de Smith sin truncarla previamente. Es m´as, como los s´ımplices de dimensi´on superior a n+ 1 no juegan ning´un papel en el c´alculo de Hn, no necesitar´ıamos calcular todo el complejo simplicial sino solo su (n+ 1)-esqueleto [40]. 22 CAP´ ITULO 2. HOMOLOG´ IA SIMPLICIAL             cd ce +bc +bc bc∗bd∗be∗−bd −be c1 0 0 0 0 d0 1 0 0 0 e0 0 1 0 0 b+c+d+e0 0 0 0 0 a0 0 0 0 0             (a) Forma normal de Smith de M1   bcd cd +bc −bd 1 ce +bc −be 0  (b) M2truncada H0∼ =Z5/Z3=Z2 H1∼ =Z2/Z=Z H2∼ =0 (c) Grupos de homolog´ıa de K Figura 2.2: C´alculo de los grupos de homolog´ıa de Kusando la forma normal de Smith. V´ease que las columnas pivotales de ˜ M1(marcadas con asterisco) se han eliminado de M2, quedando ´esta de forma autom´atica en forma normal de Smith. Cuando R=Z, el algoritmo de reducci´on a la forma normal de Smith tiene un orden superior al polin´omico. Aunque existen t´ecnicas m´as sofisticadas para obtener la reducci´on con menos operaciones [61, 62], en la pr´actica se suele preferir la reducci´on cl´asica. Si en cambio los coeficientes se toman dentro de un cuerpo k, no es necesario llevar las matrices a su forma normal de Smith, sino a su forma escalonada por columnas. En ese caso, s´olo se deben calcular los n´umeros de Betti y ´estos se extraen f´acilmente a partir de los rangos de las matrices [66]. La computaci´on es por tanto mucho m´as veloz cuando fijamos los coeficientes en un cuerpo. Veamos con un peque˜no ejemplo c´omo afecta al resultado final la elecci´on del anillo de coeficientes. Sea Kun complejo simplicial homeomorfo a la botella de Klein. Los grupos de homolog´ıa de Ksobre Zson H0=Z,H1=Z⊕Z/Z2 y H2={0}. Tenemos entonces un 1-ciclo de torsi´on ctal que cno bordea ninguna 2-cadena, pero 2csi lo hace. Si tomamos coeficientes en el cuerpo Z/Z2, el ciclo 2cse considerar´ıa nulo dentro de C1(K) y no aportar´ıa nada a im ∂2. En ese caso, los n´umeros de Betti ser´ıan β0= 1, β1= 2, β2= 1. Si en cambio tomamos los coeficientes en Z/Zppara cualquier p≥2 primo, o en otro cuerpo como QoR, el subespacio generado por 2ces el mismo que el generado por c, y entonces el ciclo ces parte del denominador de Hn. En ese caso, los n´umeros de Betti ser´ıan β0= 1, β1= 1, β2= 0. En otras palabras, al tomar los coeficientes sobre un cuerpo no podemos observar directamente la torsi´on de un complejo pero podemos deducir su existencia comparando los resultados obtenidos sobre cuerpos diferentes [75]. Cuando nuestro complejo Kest´a inmerso en R3, sus Z-m´odulos o grupos de homolog´ıa no tienen torsi´on y se puede usar una t´ecnica diferente para calcular sus n´umeros de Betti, llamada el algoritmo incremental [30], que no se apoya en c´alculos sobre matrices sino en teor´ıa de grafos. Algoritmo 2.2 (algoritmo incremental para los n´umeros de Betti).Sea Kun complejo simplicial inmerso en R3. Ordenamos los s´ımplices de Kseg´un su dimensi´on, quedando K={σ1,··· , σn}. Definimos βj= 0 para j= 0,1,2,3 2.4. COMPUTACI ´ ON DE LA HOMOLOG´ IA 23 Para todo i= 1,··· , n k= dim σi Si σipertenece a un k-ciclo de Ki={σ1,··· , σi} entonces βk=βk+ 1 y si no, βk−1=βk−1−1 En base a este algoritmo, se ve que cada vez que entra un nuevo s´ımplice al complejo pasa una de dos cosas: o bien crea una nueva clase de homolog´ıa k-´esima o bien destruye una clase (k−1)-´esima. N´otese que todos los v´ertices forman 0-ciclos al entrar en el complejo. Para detectar 1-ciclos y 2-ciclos se consideran el grafo de aristas y el grafo dual de los tri´angulos respectivamente y se aplican algoritmos eficientes de la teor´ıa de grafos [27]. Por otra parte, no es posible que se forme ning´un 3-ciclo porque entonces el conjunto de sus 3-s´ımplices ser´ıa homeomorfo a S4, que no puede estar inmerso en R3. Entonces, cada tetraedro nuevo siempre destruye una clase de homolog´ıa de dimensi´on 2. 24 CAP´ ITULO 2. HOMOLOG´ IA SIMPLICIAL Cap´ıtulo 3 Homolog´ıa persistente Sea una filtraci´on de complejos simpliciales como en la definici´on 1.15. Cada uno de los complejos anidados que la componen es una aproximaci´on de la topolog´ıa de Sa un cierto nivel de precisi´on. En la secci´on anterior hemos definido los m´odulos de homolog´ıa para estudiar la topolog´ıa de Smediante herramientas algebraicas que puedan calcularse computacionalmente. Esta herramienta opera sobre cada complejo de forma individual, pero nosotros no sabemos a priori cu´al de estas diversas aproximaciones es la m´as acertada. Para resolver este problema nuestro enfoque consistir´a en analizar la homolog´ıa del complejo en todos los niveles de la filtraci´on [33] para despu´es resumir en una sola estructura algebraica la evoluci´on de este proceso. Esta estructura que definiremos se llamar´a el m´odulo de homolog´ıa persistente. A lo largo de este cap´ıtulo supondremos que los m´odulos de cadenas y los m´odulos de homolog´ıa tienen sus coeficientes sobre un anillo gen´erico R. 3.1. Complejos y m´odulos de persistencia En una filtraci´on K, cada complejo simplicial Kucuenta con su propio complejo de cadenas Cu ∗. Los complejos de persistencia nos permitir´an reunir todas estos complejos en uno solo [75]. Definici´on 3.1. Sea K={Ku}u∈N({Ku}u∈R) una filtraci´on de tipo finito. Su complejo de persistencia Ces la familia de complejos de cadenas {Cu ∗}u∈N({Cu ∗}u∈R) con coeficientes en un anillo R, junto con los homomorfismos de complejos de cadenas inducidos por la inclusi´on iu,v :Cu ∗→Cv ∗por cada par u≤v. Los complejos de persistencia de filtraciones con par´ametro natural se pueden expresar m´as naturalmente en forma de diagrama como: C:C0 ∗ i0,1 −→ C1 ∗ i1,2 −→ C2 ∗ i2,3 −→ ··· (3.1) Si expandimos los complejos de cadenas, mostrando todos sus m´odulos de cadenas, el diagrama se puede ver como: ∂3↓∂3↓∂3↓ C0 2 i0,1 −→ C1 2 i1,2 −→ C2 2 i2,3 −→ ··· ∂2↓∂2↓∂2↓ C0 1 i0,1 −→ C1 1 i1,2 −→ C2 1 i2,3 −→ ··· ∂1↓∂1↓∂1↓ C0 0 i0,1 −→ C1 0 i1,2 −→ C2 0 i2,3 −→ ··· (3.2) 25 32 CAP´ ITULO 3. HOMOLOG´ IA PERSISTENTE (a) C´odigo de barras (b) Curvas de Betti Figura 3.4: Representaci´on visual de los c´odigos de barras y curvas de Betti de la filtraci´on de la figura 3.1. Definici´on 3.10. Sea Bn={[bi, di)}iel c´odigo de barras del m´odulo de homolog´ıa persistente PHn. Por cada intervalo [bi, di)∈B, dibujamos el punto (bi, di)en el plano real extendido [0,∞]2. Definimos el diagrama de persistencia Dncomo Dn={(bi, di)|[bi, di)∈Bn} ∪ d+, donde d+es la diagonal d+={(x, x)|x∈R+}. Por construcci´on, todos los puntos no diagonales del diagrama de persistencia se encuentran en el primer cuadrante del plano real extendido, y por encima de d+. En esta representaci´on, los puntos m´as alejados de la diagonal se corresponden con las clases de homolog´ıa m´as persistentes, mientras que los puntos cercanos a ella se corresponden con clases de corta duraci´on. Estos diagramas son muy f´aciles de representar y visualizar, y por ello es habitual dibujar los diagramas de varias dimensiones diferentes simult´aneamente en una misma gr´afica, como se puede ver en la figura 3.5a. A continuaci´on presentaremos un resultado que relaciona los diagramas de persistencia con los n´umeros persistentes de Betti, d´andole a ´estos una representaci´on visual [33]. Definici´on 3.11. Dado un punto (b, d)de un diagrama de persistencia D, definimos Tb,d como el tri´angulo dado por las desigualdades: Tb,d ={(x, y)∈[−∞,∞]2|x≥b, y < d, y ≥x}(3.8) Lema 3.2 (de los tri´angulos).Sea el m´odulo de homolog´ıa persistente n-´esima PHny sea T el multiconjunto de todos los tri´angulos Tb,d de su diagrama de persistencia Dn. El n´umero de tri´angulos de Tque se superponen sobre el punto (l, p)coincide con βl,p n. Gracias a este lema podemos decir que el diagrama de persistencia no s´olo caracteriza el nacimiento y muerte de todas las clases de homolog´ıa, sino que adem´as aporta informaci´on sobre el rango de todos los homomorfismos ιl,p dentro del m´odulo PHn. Rec´ıprocamente, conociendo el valor de cada βl,p npodemos reconstruir todos los tri´angulos Tb,d y por tanto el diagrama Dn. Es decir, ambos invariantes son equivalentes y como adem´as son equivalentes al c´odigo de barras, son completos. Por ello existen algunos autores que definen la homolog´ıa persistente n-´esima como la imagen de todos los homomorfismos ιl,p :Hl n→Hp n[4]. Por ´ultimo daremos un resultado sobre la estabilidad de los diagramas de persistencia frente a peque˜nas perturbaciones en el conjunto de puntos S. Para medir la perturbaci´on de un diagrama de persistencia, antes debemos dar una medida de similaridad entre ellos [48, 25]. 3.5. DIAGRAMAS DE PERSISTENCIA 33 (a) Diagramas de persistencia (b) Lema de los tri´angulos Figura 3.5: Diagrama de persistencia asociado a nuestra filtraci´on. A la derecha vemos c´omo es el lema de los tri´angulos sobre el diagrama de dimensi´on 0. V´ease que T0,1es un tri´angulo doble porque el c´odigo de barras contiene dos veces al intervalo [0,1) Definici´on 3.12. Sean dos diagramas de persistencia DyD0. La distancia pde Wasserstein entre DyD0es: Wp(D, D0) = ´ınf φ∈Φ"X x∈D||x−φ(x)||p ∞#1/p ,(3.9) donde Φes el conjunto de todas las biyecciones entre DyD0. La distancia bottleneck entre DyD0es: dB(D, D0) = ´ınf φ∈Φsup x∈D||x−φ(x)||∞(3.10) En ambos casos, si no existe el ´ınfimo, se dice que la distancia entre DyD0es infinita. La distancia bottleneck dBes el l´ımite de Wpcuando ptiende a infinito, y por ello algunos autores la denotan como W∞. A la hora de calcular estas biyecciones, hay que tener en cuenta que d+tambi´en es parte de D, que algunos puntos pueden ser m´ultiples y otros pueden estar en la recta del infinito. En base a esta distancia bottleneck , podemos dar el siguiente resultado [25]: Teorema 3.3 (de estabilidad para diagramas de persistencia).Sea Xun espacio triangulable y sean f, g :X→Rdos funciones de Morse. Supongamos que se construyen sendas filtraciones por conjuntos de subnivel de fygy sean Dn(f)yDn(g)sus diagramas de persistencia de dimensi´on nasociados. Entonces los diagramas de persistencia satisfacen: dB(Dn(f), Dn(g)) ≤ ||f−g||∞(3.11) Es decir, que una peque˜na perturbaci´on en la funci´on de filtro induce un cambio peque˜no en la estructura de nuestros m´odulos de homolog´ıa. Mientras que una peque˜na perturbaci´on en el conjunto de puntos Spuede cambiar sustancialmente la homolog´ıa simplicial de un complejo, la homolog´ıa persistente se mantiene estable. Este es otro motivo por el cual es preferible estudiar una base de datos mediante filtraciones de complejos simpliciales antes que hacerlo a un ´unico nivel de precisi´on. 34 CAP´ ITULO 3. HOMOLOG´ IA PERSISTENTE Cap´ıtulo 4 Homolog´ıa multipersistente En el cap´ıtulo anterior hemos definido la homolog´ıa persistente de una filtraci´on de complejos simpliciales, hemos caracterizado su estructura algebraica y hemos encontrado una parametrizaci´on discreta de la misma. En esta secci´on trataremos de seguir un desarrollo similar para el caso de filtraciones multiparam´etricas o multifiltraciones. La definici´on 3.1 de complejo de persistencia de una filtraci´on, la definici´on 3.2 de m´odulo de persistencia y la definici´on 3.3 de m´odulo de homolog´ıa persistente se extienden f´acilmente al caso multiparam´etrico cambiando u∈N(u∈R) por u∈Np(u∈Rp) y recordando que u≤vsignifica que ui≤vipara todo i= 1,··· , p. Igualmente, los conceptos de anillo graduado y m´odulo graduado se extienden de forma natural teniendo en cuenta que la suma directa que los define pasa a estar indexada sobre Np(Rp). Dado un anillo no graduado R, el anillo de polinomios con pvariables R[x1,··· , xp] es un anillo p-graduado donde los monomios son los elementos homog´eneos. Cuando Res un cuerpo, es habitual denotar Ap=R[x1,··· , xp]. En analog´ıa con el caso uniparam´etrico, tambi´en contamos con un teorema de correspondencia similar al teorema 3.1 [19]: Teorema 4.1 (de correspondencia).Existe una equivalencia αentre la categor´ıa de los m´odulos de multipersistencia de tipo finito sobre Ry la categor´ıa de los m´odulos graduados finitamente generados sobre el anillo de polinomios R[x1,··· , xp], dada por: α(M) = M u Mu,(4.1) donde el homomorfismo fu,v se corresponde con la multiplicaci´on por el monomio xv−u= xv1−u1 1···xvp−up p. Al contrario que en el caso uniparam´etrico, R[x1,··· , xp] no es un dominio de ideales principales ni siquiera cuando Res un cuerpo, y por tanto este teorema de correspondencia no nos conduce a un teorema de estructura similar al teorema 3.2. Tampoco podremos definir un invariante similar al c´odigo de barras o al diagrama de persistencia que caracterice completamente y de forma ´unica los m´odulos de multipersistencia [19]. S´ı que existe una clasificaci´on completa para los m´odulos de multipersistencia, que depende de dos invariantes discretos y un tercer invariante que en general es continuo [19]. Esta clasificaci´on no nos resulta ´util porque en la pr´actica estar´ıa sujeta a la precisi´on que seamos capaces de dar dentro del rango continuo de opciones. Por este motivo la clasificaci´on completa se descartar´a y tendremos que renunciar a las buenas propiedades que tienen los invariantes del caso uniparam´etrico. El problema de clasificar los m´odulos graduados sobre Apest´a siendo ampliamente estudiado en la actualidad, y existen diferentes propuestas para atajarlo. Dedicaremos las restantes secciones de este cap´ıtulo para hacer una breve introducci´on a algunas de ellas. 35 36 CAP´ ITULO 4. HOMOLOG´ IA MULTIPERSISTENTE (a) Etapa (0,1) (b) Etapa (1,1) (c) Etapa (2,1) (d) Etapa (0,0) (e) Etapa (1,0) (f) Etapa (2,0) Figura 4.1: Ejemplo de bifiltraci´on monocr´ıtica 4.1. Computaci´on de la homolog´ıa multipersistente Una v´ıa para el estudio de la homolog´ıa multipersistente abandona la b´usqueda de invariantes y trata de calcular MHnde una forma m´as directa, suponiendo que la multifiltraci´on Kes monocr´ıtica y por tanto los operadores borde ∂nse pueden expresar en forma de matriz. Toda la teor´ıa de esta secci´on se desarrolla en el art´ıculo [18]. El primer paso es calcular una base para el subm´odulo de los ciclos ker ∂n. La matriz Mntiene mncolumnas, a las que llamaremos F={f1, f2,··· , fmn}. El conjunto de vectores polin´omicos (q1, q2,··· , qmn)∈Amn ptales que mn X i=1 qifi= 0 forman el m´odulo de sicigias de hFi, denotado Syz(f1, f2,··· , fmn). Nuestro objetivo ser´a encontrar una base de generadores para Syz(f1, f2,··· , fmn). Para ello debemos calcular una base Gde Gr¨obner para hFiy a partir de ella usar el algoritmo de Schreyer, basado principalmente sobre el algoritmo de divisi´on de polinomios, para obtener una base de Gr¨obner para el subm´odulo Syz(f1, f2,··· , fmn) [28]. En segundo lugar hemos de calcular una base para el subm´odulo de los bordes im ∂n+1. La matriz Mn+1 tiene mn+1 columnas, a las que llamaremos F={f1, f2,··· , fmn+1 }.Fforma una base del subm´odulo im ∂n+1. El tercer paso ser´ıa comprobar si los vectores de Gpertenecen al subm´odulo hFi. Para comprobar si f∈ hFi, aplicamos el algoritmo de divisi´on ordenadamente entre los elementos de F. Al final del proceso tenemos que f= (Pmn+1 i=1 qifi)+r. Si el resto res igual a 0, entonces f∈ hFi. En cambio, el rec´ıproco no siempre es cierto en m´odulos sobre anillos de polinomios en varias variables y podr´ıa darse que r6= 0 aunque f∈ hFi. Para solventar este problema, debemos cambiar la base Fpor una base de Gr¨obner. Si Fes una base de Gr¨obner y f∈ hFi, siempre obtenemos r= 0 al aplicar el algoritmo de divisi´on, y el problema de pertenencia a un subm´odulo queda resuelto. Dada una base Fde un subm´odulo M, siempre podemos calcular una base de Gr¨obner de Ma partir de Fmediante el algoritmo de Buchberger [28]. Una vez que tenemos una base de Gr¨obner para im ∂n+1, este tercer paso es sencillo. Basta tomar cada generador de ker ∂ny dividirlo por dicha base. Si el resto res nulo, ese generador es un borde y no genera ninguna clase de MHn. Si en cambio r6= 0, podemos a˜nadir rtanto a im ∂n+1 como a MHn. 4.2. INVARIANTE DEL RANGO 37       ac ac bc bd cd a−x2−x2000 b x20−x2−x2 10 c0x2x20−x1 d000x11       (a) M1         abc ab x1 ac −x1 bc x1 bd 0 cd 0         (b) M2 Figura 4.2: Matrices de ∂1y∂2para la bifiltraci´on de la figura 4.1. Los algoritmos que resuelven estos problemas en general requieren un tiempo de computaci´on y una memoria exponencial respecto al n´umero de generadores, pero los m´odulos procedentes de multifiltraciones monocr´ıticas gozan de unas propiedades de homogeneidad que nos permiten reducir el tiempo y la memoria hasta un orden polin´omico. La soluci´on computada por este m´etodo para la homolog´ıa multipersistente es compacta y completa, pero no es un invariante dado que las bases de Gr¨obner calculadas son diferentes seg´un la base con la que se inicie el algoritmo de Buchberger y la relaci´on de orden que se defina entre monomios. No obstante, una base de Gr¨oebner siempre se puede convertir en una base reducida de Gr¨obner, que s´ı es ´unica para un orden fijado entre los monomios [28]. 4.2. Invariante del rango Dada una multifiltraci´on de tipo finito cualquiera, podemos dar un invariante que compare la homolog´ıa de todas las parejas de complejos posibles, al que llamaremos el invariante del rango. La descripci´on de este invariante se detalla en el art´ıculo [19]. En lo que resta de cap´ıtulo supondremos que la multifiltraci´on tiene par´ametros naturales, pero tambi´en es posible considerar par´ametros reales. Definici´on 4.1. Sea ¯ N=N∪{∞} y sea Dp={(u, v)|u∈Np, v ∈¯ Np, u ≤v}. Dados dos pares (u, v),(u0, v0)∈Dp, definimos (u, v)(u0, v0)si y solo si u≤u0yv0≤v. En el caso uniparam´etrico, D1es el espacio donde podemos situar los puntos de un diagrama de persistencia y (b, d)(b0, d0) si y solo si Tb0,d0⊂ Tb,d. Definici´on 4.2. Sea MHnun m´odulo de homolog´ıa multipersistente. Definimos el invariante del rango de MHncomo la aplicaci´on ρn:Dp→Ndada por ρn(u, v) = rank(ι(u,v):Hu n→Hv n). El invariante del rango es la versi´on multiparam´etrica de los n´umeros persistentes de Betti. Una propiedad importante del invariante del rango es que conserva el orden [19]. Lema 4.1. Si (u, v)(u0, v0), entonces ρn(u, v)≤ρn(u0, v0). Demostraci´on. Si (u, v)(u0, v0), entonces se tiene la cadena de desigualdades u≤u0≤v0≤v. Por tanto se tiene que ι(u,v)=ι(v0,v)◦ι(u0,v0)◦ι(u,u0). La prueba concluye trivialmente sabiendo que dadas dos transformaciones lineales fyg, se tiene que rank(g◦f)≤rank(f),rank(g). Mientras que en el caso uniparam´etrico el invariante del rango es completo, para p≥2 no lo es y podr´ıan haber dos m´odulos de persistencia no isomorfos que compartan la misma ρ. El invariante del rango cuenta con su propio teorema de estabilidad, que nos permite saber que es robusto frente a peque˜nas perturbaciones en S[20]. 38 CAP´ ITULO 4. HOMOLOG´ IA MULTIPERSISTENTE Teorema 4.2 (de estabilidad del invariante del rango).Sea Xun espacio triangulable y sean f, g :X→Rpdos funciones continuas. Supongamos que se construyen sendas filtraciones por conjuntos de subnivel de fygy sean ρn(X, f)yρn(X, g)sus invariantes del rango de dimensi´on nasociados. Entonces existe una funci´on de distancia dmatch dentro del conjunto {ρn(X, φ)|φ:X→Rpcontinua}tal que: dmatch(ρn(X, f), ρn(X, g)) ≤m´ax x∈X||f(x)−g(x)||∞(4.2) Para calcular el invariante del rango, debemos considerar todos los pares (u, v)∈Ddy calcular ρn(u, v) para cada uno. Para ello podemos considerar la filtraci´on que contiene ´unicamente a los complejos KuyKvy aplicar el algoritmo de persistencia 3.1. El rango ρn(u, v) coincidir´ıa con el n´umero de intervalos en su c´odigo de barras que pasen a la vez por 0 y 1 [18]. Otra v´ıa para calcular el invariante del rango consiste en calcular la homolog´ıa multipersistente con los algoritmos de la secci´on 4.1 y calcular cada rango ρn(u, v) con un algoritmo espec´ıfico que toma como valores de entrada una base de Gr¨obner Gde im ∂n+1 y una base Z de sicigias de ker ∂n[18]. En ambos casos, el invariante del rango completo requiere de un gran espacio de almacenaje y un gran tiempo de c´alculo. En un caso extremo donde cada coordenada cr´ıtica de nuestra multifiltraci´on solo aporta un s´ımplice nuevo y hay ms´ımplices en total, pueden existir hasta un total de mpcomplejos simpiciales distintos dentro de la multifiltraci´on. En total hay O(m2p) pares de coordenadas comparables y el invariante del rango requiere un espacio de almacenaje exponencial. Por tanto es preferible solo calcularlo para alg´un subconjunto de pares de Dp, como el conjunto de los pares triviales (u, u), que extiende las curvas de Betti al caso multiparam´etrico. 4.3. Diagrama de persistencia multiparam´etrico Los diagramas de persistencia multiparam´etricos se construir´an como una especie de “derivada” del invariante del rango ρn, de forma tal que el diagrama indica en qu´e pares de Dpla funci´on ρncambia. Esta secci´on se basa en el art´ıculo [47], y para su desarrollo cambiaremos la definici´on 4.1 para Dp, eliminando los puntos infinitos. Definici´on 4.3. Definimos el orden parcial ≤en Dpde forma tal que (u, v)≤(u0, v0)si existen ε, δ ≥0tales que u0=u+ε·~ 1yv0=v+δ·~ 1 Definici´on 4.4. Sea A(Dp)el conjunto de todas las matrices α:Dp×Dp→Ztales que α(I, J) = 0 si I6≤ J. La multiplicaci´on de dos matrices αyβtiene por f´ormula: αβ(I, K) = X I≤J≤K α(I, J)·β(J, K) (4.3) El conjunto A(Dp) es cerrado para la suma, el producto por un escalar y la multiplicaci´on. Definici´on 4.5. Hay dos matrices de A(Dp)particularmente interesantes. La funci´on zeta, dada por ζ(I, K) = 1si I ≤K 0cc. (4.4) La matriz de M¨obius µ, dada de forma inductiva por µ(I, K) =    1si I =K −PI≤J≤Kµ(I, J)si I ≤K 0cc. (4.5) 4.3. DIAGRAMA DE PERSISTENCIA MULTIPARAM ´ ETRICO 39 Ambas funciones son invertibles, y de hecho ζµ =µζ = 1. Definici´on 4.6. Sea MHnun m´odulo de homolog´ıa multipersistente y sea ρnsu invariante del rango. Definimos el diagrama de persistencia de MHncomo la inversi´on de M¨obius de ρn, a la que denotaremos ∂ρn, y que tiene por f´ormula: ∂ρn(K) = X J≤K µ(J, K)ρn(J) (4.6) El diagrama de persistencia es nulo salvo para un n´umero finito de elementos de Dp, y estaremos interesados en encontrar el soporte de ∂ρn. En el caso uniparam´etrico, la funci´on ∂ρn vale 1 en los puntos simples del diagrama D, 2 en los puntos dobles y as´ı sucesivamente. En el caso general, ∂ρnes no nula en los pares (u, v)∈Dpdonde ues una coordenada de nacimiento de una clase y ves una coordenada de muerte, aunque no necesariamente est´en ambas referidas a la misma clase. El diagrama de persistencia ∂ρnes una derivada de ρnen tanto que ρnse puede reconstruir a partir de ∂ρnmediante sumas, y de hecho se tiene que ρn(K) = PJ≤K∂ρn(J). Esta funci´on tambi´en es un invariante incompleto para MHn, y al igual que el invariante del rango, tambi´en requiere en el peor de los casos un tiempo de c´alculo y una memoria exponenciales. No obstante, como solo nos interesan sus valores no nulos que en general son una minor´ıa, el espacio de almacenamiento puede ser mucho menor. 40 CAP´ ITULO 4. HOMOLOG´ IA MULTIPERSISTENTE Cap´ıtulo 5 Inferencia estad´ıstica aplicada a la homolog´ıa persistente En el cap´ıtulo 3 vimos la persistencia como una herramienta para analizar la evoluci´on de la homolog´ıa de un complejo simplicial a lo largo de una filtraci´on. Su estructura algebraica se correspond´ıa con un m´odulo graduado, y al calcular sus coeficientes sobre el dominio de ideales principales k[t], encontramos dos parametrizaciones discretas equivalentes: el c´odigo de barras y el diagrama de persistencia. Nuestro objetivo ahora ser´a aplicar t´ecnicas de inferencia estad´ıstica sobre los m´odulos de persistencia. La principal desventaja de los c´odigos de barras o los diagramas de persistencia a la hora de hacer estad´ıstica con ellos es que no est´an sujetos a nociones de aditividad o producto por un escalar, haciendo dif´ıcil incluso el c´alculo de una media aritm´etica entre diagramas. Algunos autores propusieron usar la media de Fr´echet para resolver este problema. Definici´on 5.1. Dados los diagramas de persistencia D1, D2,··· , Dp, se dice que el diagrama ¯ Des una media de Fr´echet de este conjunto si verifica: ¯ D= arg m´ın Y 1 p p X i=1 dB(Di, Y ) (5.1) La b´usqueda de la media de Fr´echet requiere aplicar un algoritmo de descenso del gradiente, que no siempre alcanza su ´optimo. Adem´as, esta definici´on de media es inestable frente a peque˜nas perturbaciones y en general no es ´unica [65]. Por ello definiremos en esta secci´on nuevos objetos que nos permitan analizar la homolog´ıa persistente desde un punto de vista estad´ıstico. 5.1. Paisajes de persistencia y siluetas Un paisaje de persistencia es un conjunto de funciones lineales a trozos que presenta una parametrizaci´on alternativa a los diagramas de persistencia. Su definici´on [24] se basa en unas ciertas funciones triangulares extra´ıdas de un diagrama de persistencia. Definici´on 5.2. Sea Dun diagrama de persistencia y sea p= (b, d)∈D. Si pes un punto finito no diagonal, definimos su funci´on triangular Λpcomo: Λp(t) =            t−b si t ∈[b, (b+d)/2] d−t si t ∈[(b+d)/2, d] 0c.c. (5.2) 41 48CAP´ ITULO 5. INFERENCIA ESTAD´ ISTICA APLICADA A LA HOMOLOG´ IA PERSISTENTE (a) Banda de confianza para un paisaje (b) Banda de confianza para una silueta Figura 5.4: Bandas de confianza basadas en bootstrap multiplicativo. La muestra de paisajes y siluetas se ha obtenido a˜nadiendo peque˜nas perturbaciones al diagrama de la figura 5.1. V´ease que en ambos casos la banda de confianza contiene completamente a la funci´on esperada. Generamos ξ1,··· , ξn∼N(0,1) independientes. Calculamos θ∗=1 √n||Pn i=1 ξiλi(t)−¯ λn(t) ˆσ(t)||∞. Repetimos los dos pasos anteriores Bveces, obteniendo θ∗ 1,··· , θ∗ B. Calculamos el cuantil qα= ´ınf{q:1 BPB j=1 I(θ∗ j≥q)≤α}. La banda de confianza es C(¯ λn)=[¯ λn−qαˆσ(t) √n,¯ λn+qαˆσ(t) √n] Esta banda de confianza adaptativa tambi´en es asint´oticamente correcta [24], y se ajusta m´as a la variabilidad de los paisajes en cada tdel dominio. 5.5. Conjuntos de confianza En este apartado trataremos el problema de la detecci´on de ruido en un diagrama de persistencia. Para ello debemos encontrar el equivalente a un intervalo de confianza para diagramas de persistencia, al que llamaremos un conjunto de confianza al nivel (1 −α). Supongamos que los diagramas provienen de filtraciones basadas en funciones como en las secciones anteriores. Sea Dfel diagrama de persistencia no observable asociado a la filtraci´on por la funci´on f. Nuestro conjunto Snos da unos estimadores ˆ fyˆ D. Aplicando el algoritmo bootstrap 5.1 para funciones, encontramos un valor qαtal que asint´oticamente P(√N· ||f− ˆ f||∞> qα)→α. Gracias al teorema de estabilidad 3.3 para diagramas, tenemos que: P√N·dB(Df,ˆ D)> qα≤P√N·||f−ˆ f||∞> qα→α, (5.12) y hemos encontrado un conjunto de confianza asint´oticamente correcto para ˆ D[37]. Si llamamos cα=qα/√N, el conjunto de confianza C(cα) es la bola cerrada C(cα) = {D:dB(ˆ D, D)≤cα}(5.13) El conjunto C(cα) se puede visualizar centrando en cada punto de ˆ Duna caja de lado 2cα, como se hace en la figura 5.5a. Un diagrama Dpertenece a C(cα) si tiene un punto en cada caja. Si la caja de un punto concreto interseca a la diagonal, entonces existe un diagrama en C(cα) donde este punto se integra en la diagonal y desaparece. Entonces, dado nuestro diagrama ˆ D, podemos eliminar todos los puntos cuyas cajas corten con la diagonal, y por la construcci´on 5.6. AN ´ ALISIS DE M ´ AXIMA PERSISTENCIA 49 (a) Conjunto de confianza al 95 % (b) Bandas de confianza al 95 % Figura 5.5: Conjunto y bandas de confianza para el diagrama de la figura 5.3. que hemos hecho sabemos que el nuevo diagrama no es significativamente distinto al original porque est´a dentro de C(cα). Una alternativa a las cajas es dibujar sobre la diagonal una banda con altura 2cα, como se hace en la figura 5.5b, y eliminar los puntos que caigan sobre esta banda, es decir, que tengan una persistencia menor que 2cα. Podemos hacer un conjunto de confianza m´as preciso si en vez de usar el algoritmo bootstrap 5.1 sobre las funciones lo usamos directamente sobre la distancia bottleneck . Este algoritmo se llama bottleneck bootstrap. Algoritmo 5.4 (bottleneck bootstrap para diagramas).Sea S={x1,··· , xN}y sean ˆ fyˆ Dlos estimadores de fyDfcalculados sobre S. Tomamos una submuestra x∗ 1,··· , x∗ Nde Scon reemplazamiento y sobre ella calculamos el estimador ˆ f∗. A partir de ˆ f∗calculamos el diagrama de persistencia D∗. Teniendo D∗, calculamos θ∗=√N·dB(D∗,ˆ D). Repetimos los dos primeros pasos Bveces, obteniendo D∗ 1,··· , D∗ Byθ∗ 1,··· , θ∗ B. Calculamos el cuantil qα= ´ınf{z:1 BPB j=1 I(θ∗ j≥q)≤α}, y cα=qα/√N. El conjunto de confianza es C(cα) = {D:dB(ˆ D, D)≤cα} Este algoritmo es correcto independientemente de que fsea una funci´on DTM o un estimador de densidad, y da un conjunto de confianza m´as preciso porque el algoritmo bootstrap 5.1 para funciones aprovecha que gracias al teorema 3.3 de estabilidad dB(Df,ˆ D)≤ ||f−ˆ f||∞y encuentra una cota superior para ||f−ˆ f||∞, lo cual da un conjunto de confianza muy conservativo [22]. 5.6. An´alisis de m´axima persistencia A lo largo de todo este cap´ıtulo hemos supuesto que las filtraciones se constru´ıan a trav´es de funciones. En particular nos hemos centrado en la funci´on DTM, que depende de un par´ametro m0, y del estimador de densidad con n´ucleo gaussiano, que depende de un par´ametro hque define su varianza. La elecci´on de un buen par´ametro para estas funciones es a´un un problema sin resolver. 50CAP´ ITULO 5. INFERENCIA ESTAD´ ISTICA APLICADA A LA HOMOLOG´ IA PERSISTENTE Definici´on 5.7. Sea una filtraci´on dada por una funci´on fque depende de un par´ametro m. Para cada m, sean l1(m), l2(m),··· las persistencias de los diferentes puntos de su diagrama, y sea c(m)la altura de su banda de confianza. Definimos: N(m) = # {i:l(i)> c(m)} S(m) = Pi[l(i)−c(m)]+ N(m) cuenta el n´umero de clases de homolog´ıa significativas, mientras que S(m) nos indica cu´anto persisten estas clases. En el art´ıculo [22] se propone escoger el par´ametro mque maximice N(m), S(m) o bien la vida significativa media S(m)/N(m). Cap´ıtulo 6 Aplicaciones del an´alisis topol´ogico de datos A pesar de lo reciente de esta rama de estudio, el an´alisis topol´ogico de datos ya ha sido aplicado en todo tipo de ´areas, tales como el reconocimiento de formas [31], an´alisis de redes [44], an´alisis del lenguaje natural [72], an´alisis de im´agenes [17], biolog´ıa computacional [42], dise˜no de f´armacos [69], etc. En este cap´ıtulo haremos una breve introducci´on a algunas de sus aplicaciones. 6.1. Machine learning basado en homolog´ıa persistente En los cap´ıtulos anteriores hemos descrito diferentes objetos que representan la homolog´ıa persistente, tales como los c´odigos de barras, las curvas de Betti, los diagramas de persistencia, los paisajes y las siluetas. Estos objetos pueden ser utilizados como variable de entrada en diversos algoritmos de aprendizaje autom´atico, bien en su forma original o bien tras aplicarles alguna transformaci´on. La forma m´as simple de extraer caracter´ısticas a partir de la homolog´ıa persistente es hacer un resumen estad´ıstico del c´odigo de barras. Por ejemplo extraer el m´aximo, el m´ınimo, la media y la varianza de los tiempos de nacimiento, los tiempos de muerte o las persistencias de cada intervalo [14]. Otros autores proponen construir vectores num´ericos a partir de nuestros objetos con el m´etodo del empaquetado o binning. La idea es discretizar el rango de la filtraci´on en varios intervalos del mismo tama˜no y asignar a cada uno un valor que resuma una propiedad suya. Para un c´odigo de barras, podemos dividir su dominio en nintervalos y crear tres vectores con ncomponentes que cuenten respectivamente cu´antos intervalos nacen, cu´antos mueren y cu´antos hay en total en cada intervalo. Para la curva de Betti, el binning es tan simple como almacenar sus valores para una serie de puntos equidistantes. Para un paisaje se puede repetir este proceso una vez por cada componente, obteniendo un vector distinto para cada una y para siluetas se puede aplicar tambi´en. Para un diagrama de persistencia rotado, podemos dar una malla con n×ncuadrados y crear un vector de n2componentes que cuente cu´antos puntos hay en cada cuadrado. Estos vectores num´ericos se pueden usar para algoritmos de clustering, m´aquinas de vector soporte, ´arboles de decisi´on, redes neuronales, etc [56]. En el caso de una bifiltraci´on o una trifiltraci´on, podemos crear matrices en dos y tres dimensiones si almacenamos el n´umero de Betti para cada pareja o terna de par´ametros. Al contrario que las representaciones vectoriales antes explicadas, estas matrices pueden ser tratadas como im´agenes y se pueden combinar con redes neuronales convolucionales [15]. 51 52 CAP´ ITULO 6. APLICACIONES DEL AN ´ ALISIS TOPOL ´ OGICO DE DATOS Otro enfoque para poder aplicar m´aquinas de vector soporte a diagramas o paisajes de persistencia pasa por definir sobre ellos diferentes medidas de semejanza o funciones n´ucleo cuyos valores se puedan almacenar en una matriz [56]. Las medidas de semejanza m´as conocidas son la distancia pde Wasserstein y la distancia bottleneck entre diagramas, cuyas f´ormulas se describen en la definici´on 3.12. Tambi´en podemos definir n´ucleos entre diagramas de persistencia. Definici´on 6.1. Sean dos diagramas de persistencia DyD0. El n´ucleo PSSK (persistent scale space kernel [58]) tiene por f´ormula: κPSSK(D, D0, σ) = 1 8πσ X x∈D,x0∈D0 e−||x−x0||2 ∞ σ−e−||¯x−x0||2 ∞ σ,(6.1) donde ¯xes el punto sim´etrico de xrespecto a la diagonal. El n´ucleo GTK (geodesic topological kernel [52]) tiene por f´ormula: κGTK(D, D0, σ) = e{1 hW2(D,D0)2}(6.2) El n´ucleo GLK (geodesic laplacian kernel [52]) tiene por f´ormula: κGLK(D, D0, σ) = e{1 hW2(D,D0)}(6.3) Tambi´en es posible definir n´ucleos para paisajes o siluetas. Definici´on 6.2. Sean dos paisajes de persistencia m-´esimos λmyλ0 m. El n´ucleo GPHK (global persistent homology kernel [56]) tiene por f´ormula: κGPHK(λm, λ0 m) = hλm, λ0 mi=ZR λm(x)·λ0 m(x)dx (6.4) Utilizando herramientas como ´estas, podemos combinar nuestros resultados topol´ogicos con cualquier algoritmo de aprendizaje autom´atico que necesitemos. 6.2. An´alisis de formas En los ´ultimos a˜nos se est´an creando bases de datos de formas tridimensionales cada vez m´as grandes, siendo necesario desarrollar nuevos algoritmos para su procesamiento. En el art´ıculo [11], se emplea el m´etodo de las bolsas de palabras para la clasificaci´on de estas formas. Supongamos que tenemos una base con Mformas 3Dalmacenadas como complejos simpliciales de dimensi´on 2 con v´ertices Vi. El primer paso consiste en definir unas ciertas funciones fi:Vi→Rnque describan la geometr´ıa local de la forma en torno a cada v´ertice. Despu´es estas funciones se proyectan sobre un conjunto de puntos W={wk|k= 1,··· , K}al que llamaremos diccionario o codebook. La composici´on de las ficon la proyecci´on sobre Wson las llamadas funciones de palabras Ci:Vi→RK, y a los puntos de RKque forman su imagen se les llama simplemente palabras. El ´ultimo paso de este proceso consiste en encontrar un buen resumen de todas estas funciones Cien forma de un ´unico vector compacto P= (Pi)i=1,···,M que pueda ser usado por ejemplo por una m´aquina de vector soporte. Este paso se conoce como pooling y se puede realizar de varias formas. Dos muy comunes son el pooling de la suma y el del m´aximo [70], cuyas componentes Pivienen dadas por: 6.3. AN ´ ALISIS DE SERIES TEMPORALES 53 Figura 6.1: Ejemplo de una funci´on de palabras sobre tres formas distintas de la base SHREC 2014. Imagen extra´ıda de [11] Pi= (SumPool(Ci,1),··· , SumPool(Ci,K) = X x∈Vi Ci(x)1,··· ,X x∈Vi Ci(x)K! Pi= (MaxPool(Ci,1),··· , MaxPool(Ci,K) = m´ax x∈Vi Ci(x)1,··· ,m´ax x∈Vi Ci(x)K El pooling del m´aximo es m´as robusto que el de la suma, y hace un resumen de los m´aximos locales de las fi, sin importar la “prominencia” del pico que representan. Los autores del art´ıculo [11] proponen un nuevo m´etodo m´as robusto que los anteriores basado en homolog´ıa persistente de dimensi´on 0, llamado pooling topol´ogico, que s´ı tiene en cuenta la prominencia de estos m´aximos. Definici´on 6.3. Dada una funci´on fsobre un grafo G, definimos: TopoPool(f)i=pi(Df), donde Dfes el diagrama de persistencia de dimensi´on 0 asociado a los conjuntos de supernivel de fypies la persistencia de su i-´esimo punto m´as alejado de la diagonal. Si Dfno tuviera i puntos, se asigna el valor 0. Gracias al teorema de estabilidad 3.3 de diagramas, el pooling topol´ogico es estable y ofrece mejores prestaciones que los m´etodos convencionales. La comparaci´on entre los tres m´etodos de pooling se hizo usando la base de datos SHREC 2014 [55], que contiene 400 formas de 40 personas en 10 poses distintas. 6.3. An´alisis de series temporales Existen diferentes v´ıas para usar la homolog´ıa persistente y sus representaciones en el estudio de series temporales tanto num´ericas como categ´oricas [57]. Los datos de una serie temporal no pueden ser usados directamente para construir una filtraci´on, sino que deben ser transformados previamente en una nube de puntos. Para ello contamos con dos m´etodos distintos: la inmersi´on de retrasos de Takens y el metodo de ventanas deslizantes. Sea una serie temporal {xt|t= 1,2,··· , T}. La inmersi´on de retrasos de Takens [63] convierte la serie en una nube de puntos en Rddados por vi= (xi, xi+τ, xi+2τ,··· , xi+(d−1)τ), donde dindica la dimensi´on de los puntos creados y τes un par´ametro de retraso. Existen multitud de propuestas para escoger unos valores d, τ apropiados. 54 CAP´ ITULO 6. APLICACIONES DEL AN ´ ALISIS TOPOL ´ OGICO DE DATOS (a) Serie temporal (b) Nube de puntos (c) Diagrama de persistencia Figura 6.2: Detecci´on de se˜nales peri´odicas con homolog´ıa persistente. La nube de puntos se ha construido con el m´etodo de retrasos de Takens. El m´etodo de las ventanas deslizantes es una alternativa al m´etodo de Takens para detectar periodicidades en series con ruido, con demasiada separaci´on temporal o con separaciones desiguales [54]. El primer paso, opcional, es eliminar ruido aplicando por ejemplo una media m´ovil. El segundo es interpolar esta serie suavizada con un spline c´ubico g. El tercero es seleccionar Npuntos equidistantes 1 = i1< i2<··· < iN=T, dando lugar a una nueva serie {g(i1),··· , g(iN)}. El cuarto paso es aplicar el m´etodo de Takens a esta nueva serie con unos d, τ adecuados. El ´ultimo paso, tambi´en opcional, es estandarizar la nube de puntos obtenida. Si nuestra serie temporal es categ´orica o contiene fuertes discontinuidades, podemos aplicarle una transformada de Walsh-Fourier antes de aplicarle alguno de los m´etodos anteriores [60]. Una vez construida una nube de puntos, podemos calcular una filtraci´on sobre ella con cualquiera de los m´etodos explicados en el cap´ıtulo 1 y a partir de ella obtener un diagrama de persistencia, un paisaje, una silueta o cualquier otra representaci´on que nos pueda ser ´util. La homolog´ıa persistente ha resultado ser una herramienta muy interesante en el estudio de series temporales, dando buenos resultados en la detecci´on de se˜nales peri´odicas ocultas, clustering, clasificaci´on de series temporales y detecci´on temprana de cambios estructurales. Por ejemplo, un ciclo de dimensi´on 1 de larga persistencia indica la existencia de una se˜nal peri´odica en la serie, como se puede ver en la figura 6.2. Se puede encontrar un resumen de las otras aplicaciones en el art´ıculo [57], que sirve de introducci´on a toda esta teor´ıa. 6.4. Algoritmo Mapper en biomedicina El an´alisis topol´ogico de datos fue usado en [51] para el estudio de informaci´on gen´omica en c´elulas tumorales de pacientes con c´ancer de mama, obteniendo resultados nuevos que por ejemplo el clustering no fue capaz de detectar. Este estudio se realiz´o con una base de datos de 295 pacientes con c´ancer de mama [67]. A cada una de ellas se les analizaron 262 genes de sus c´elulas tumorales, obteniendo los vectores T1,··· , T295. Cada vector se compar´o con un modelo lineal que representa una situaci´on de salud, obteniendo la descomposici´on Ti=Ni+Di, donde Nies la componente normal de TiyDies su componente de enfermedad. Unos Dialejados del vector nulo se corresponden con fuertes desviaciones con respecto al estado normal de salud. Una vez separados los datos en componente normal y componente enferma, se descartaron los genes que no presentaban desviaciones significativas del estado de salud. Por ´ultimo se aplic´o el algoritmo Mapper 1.2 a los datos usando como funciones de filtro potencias k-´esimas de la norma de Lp. Para calcular los clusters del algoritmo Mapper, usamos la distancia de correlaci´on entre los datos. Podemos ver el grafo resultante al usar p= 2 y k= 4 en la figura 6.3. Este grafo tiene tres ramas principales. La rama de la izquierda, con los puntos en azul, se corresponde con los tumores que tienen un tipo gen´omico parecido al estado normal de salud. La rama superior 6.5. SOFTWARE PARA EL AN ´ ALISIS TOPOL ´ OGICO DE DATOS 55 Figura 6.3: Realizaci´on del algoritmo Mapper sobre la base de datos del c´ancer de mama. Imagen extra´ıda de [51] corresponde con los tumores de la secuencia ER−, con una gran desviaci´on con respecto al estado normal y de los que la mayor´ıa son tumores basales. Existe una tercera secuencia larga de tumores con una gran desviaci´on con respecto al estado normal, la secuencia ER+. Se puede observar dentro de la secuencia ER+un subgrupo de puntos rodeado por dos ´areas m´as dispersas, al que los autores llamaron c−MY B+. El subgrupo c−MY B+, una vez detectado, ha sido estudiado y ha demostrado tener bastante uniformidad en su estructura molecular, un 100 % de supervivencia y un riesgo nulo de reaparici´on. Tambi´en se puede detectar usando otros par´ametros para el algoritmo Mapper y se ha validado aplicando el mismo procedimiento a otros conjuntos de datos de c´ancer de mama. Es decir, estos tumores realmente forman un subgrupo totalmente diferenciado de los ya conocidos con anterioridad. En cambio, el subgrupo c−MY B+no puede ser distinguido con t´ecnicas m´as antiguas como los algoritmos de clustering. En la figura 6.4, se puede ver c´omo el clustering dispersa los tumores del subgrupo c−MY B+en diversos clusters sin encontrar una relaci´on significativa entre ellos. Esto nos da a entender que el an´alisis topol´ogico de datos puede aportar informaci´on realmente novedosa y vale la pena profundizar en este tipo de t´ecnicas para otros estudios similares. 6.5. Software para el an´alisis topol´ogico de datos Terminamos este ´ultimo cap´ıtulo haciendo un repaso a los distintos softwares que se pueden usar para el c´alculo de las herramientas que hemos definido a lo largo de todo este trabajo. Entre ellos se encuentran Javaplex [3], Perseus [50], Dipha [7], Dionysus [1], jHoles [9], GUDHI [64], Ripser [6], PHAT [8], DIPHA [7] y R. Centraremos nuestra atenci´on en el software R, enfocado a usuarios con conocimientos de estad´ıstica y extensamente estudiado en este doble 56 CAP´ ITULO 6. APLICACIONES DEL AN ´ ALISIS TOPOL ´ OGICO DE DATOS Figura 6.4: Dendrograma de la base de datos tras aplicar un algoritmo de clustering. Las barras rojas representan tumores del subgrupo c−MY B+. Imagen extra´ıda de [51] grado. La librer´ıa TDA de R fue desarrollada con vistas al tratamiento estad´ıstico de la homolog´ıa persistente [36]. Con ella se pueden calcular filtraciones basadas en el complejo alpha y Vietoris-Rips, filtraciones basadas en funciones, homolog´ıa persistente con coeficientes en Z/Z2, c´odigos de barras y diagramas de persistencia, distancias bottleneck y de Wasserstein, paisajes, siluetas, bandas de confianza para funciones, bandas para paisajes y siluetas con el bootstrap multiplicativo, conjuntos de confianza con el bottleneck bootstrap y se puede hacer an´alisis de m´axima persistencia. Casi todas las im´agenes que figuran en este trabajo han sido creadas con las funciones que ofrece esta librer´ıa. La librer´ıa TDAmapper [53] sirve para crear el grafo Mapper de un conjunto de datos para una o dos funciones de filtro. Los resultados de estas funciones se pueden combinar con las funciones de la librer´ıa igraph para obtener una representaci´on visual del grafo. Existen algunas funcionalidades que a´un no han sido implementadas en este lenguaje de programaci´on, como el c´alculo de un complejo de ˇ Cech o un complejo witness, el c´alculo de la homolog´ıa simplicial con coeficientes arbitrarios o el c´alculo de multifiltraciones. Para el c´alculo de los complejos de ˇ Cech y witness podemos usar la librer´ıa GUDHI de C++, o bien la librer´ıa JavaPlex de Java. Para el c´alculo de la homolog´ıa simplicial con coeficientes en un anillo cualquiera podemos usar la librer´ıa CHomP de C++ o el software SageMath [73]. Para el c´alculo de la homolog´ıa multipersistente tal y como se hace en la secci´on 4.1 podemos usar la librer´ıa CoCoA de C++ [2] y el software Macaulay2 [39]. Concluimos presentando en el anexo una propuesta de elaboraci´on propia para la implementaci´on en R del c´alculo de complejos de ˇ Cech, complejos witness y una funci´on para representar complejos simpiciales cuyos puntos est´an en el plano. Bibliograf´ıa [1] Dionysus, the persistent homology software. Available at https://www.mrzv.org/ software/dionysus/. [2] J. Abbott, A. M. Bigatti, and L. Robbiano. CoCoA: a system for doing Computations in Commutative Algebra. Available at http://cocoa.dima.unige.it. [3] H. Adams, A. Tausz, and M. Vejdemo-Johansson. Javaplex: A research software package for persistent (co) homology. In International Congress on Mathematical Software, pages 129–136. Springer, 2014. [4] American Mathematical Society and A. J. Zomorodian, editors. Advances in applied and computational topology: AMS Short Course on Computational Topology, January 4- 5, 2011, New Orleans, Louisiana. Number v. 70 in Proceedings of symposia in applied mathematics. American Mathematical Society, Providence, R.I, 2012. [5] M. Artin. Algebra. Pearson, 2 edition, 2018. [6] U. Bauer. Ripser: a lean c++ code for the computation of vietoris–rips persistence barcodes. Software available at https://github. com/Ripser/ripser, 436, 2017. [7] U. Bauer, M. Kerber, and J. Reininghaus. Distributed computation of persistent homology. In 2014 proceedings of the sixteenth workshop on algorithm engineering and experiments (ALENEX), pages 31–38. SIAM, 2014. [8] U. Bauer, M. Kerber, J. Reininghaus, and H. Wagner. Phat–persistent homology algorithms toolbox. Journal of symbolic computation, 78:76–90, 2017. [9] J. Binchi, E. Merelli, M. Rucco, G. Petri, and F. Vaccarino. jholes: A tool for understanding biological complex networks via clique weight rank persistent homology. Electron. Notes Theor. Comput. Sci., 306:5–18, 2014. [10] J.-D. Boissonnat, O. Devillers, and S. Hornus. Incremental construction of the delaunay triangulation and the delaunay graph in medium dimension. In Proceedings of the twentyfifth annual symposium on Computational geometry, pages 208–216, 2009. [11] T. Bonis, M. Ovsjanikov, S. Oudot, and F. Chazal. Persistence-based pooling for shape pose recognition. In International workshop on computational topology in image context, pages 19–29. Springer, 2016. [12] N. Bourbaki. Elements of mathematics. algebra, part i: Chapters 1-3. Translated from the French, Hermann, Paris, 1974. [13] P. Bubenik. Statistical topological data analysis using persistence landscapes. The Journal of Machine Learning Research, 16(1):77–102, 2015. 57 64 AP ´ ENDICE A. FUNCIONES DE R meb Dado un conjunto Ade puntos en Rny un nivel de precisi´on, esta funci´on te da el centro y el radio aproximados de la menor bola eucl´ıdea que envuelve al conjunto A.Adebe ir en forma de matriz donde cada fila representa un punto. El algoritmo est´a basado en [71]. meb = function(A , epsilon =0.0001){ ref=rep(0,dim(A)[2]) sqdist=function (u) sum ((uref )^2) alpha = which .max(apply (A,1,sqdist)) ref=A[alpha ,] beta=which .max (apply (A,1,sqdist)) center =(A[alpha ,]+A[beta ,])/2 gamma = sqdist (A[ beta ,]) ref = center cappa = which .max(apply (A,1,sqdist)) delta =( sqdist (A[ cappa ,])/ gamma ) -1 while ( delta > (1+epsilon)^2-1){ lambda = delta /(2+2* delta ) gamma =(1-lambda)*gamma +lambda*(1-lambda )*( 1+ delta )* gamma center=(1- lambda )* center + lambda *A[cappa ,] ref = center cappa = which .max(apply (A,1,sqdist)) delta = sqdist (A[ cappa ,])/ gamma -1 } result=list( center = center , radius = sqrt(gamma *(1+ delta ))) result } 65 cechFiltration Dado un conjunto de puntos X, un εm´aximo y una dimensi´on m´axima para los s´ımplices, esta funci´on te calcula la filtraci´on de ˇ Cech en un formato compatible con las funciones de la librer´ıa TDA. Est´a basado en los algoritmos de [74] cechFiltration = function (X, maxdimension , maxscale ){ #En primer lugar definimos una serie de funciones auxiliares #1. Aristas del complejo cech . edges = function (X, max .scale ){ edges =c() l=dim(X)[1] for (i in 1:(l-1)){ for (j in (i+1):l){ d=sqrt(sum((X[i,]-X[j ,])^ 2)) if (d <=2*max .scale ){ edges = rbind (edges ,c(j,i,d/2)) } } } if (! is.null( edges )){ edges = matrix(edges,ncol=3) edges = edges [ order ( edges [,3]),]} edges } #2. Vecinos menores get_lower _ nbrs = function( edgemat , num ){ edgemat [ edgemat [,1]== num ,2]} #3. Engrandecer ( encontrar s´ımplices de una dimensi´on mayor ) cech . enlarge = function ( simplices , edges ){ d=dim( simplices ) result=c() for (i in 1:d[1]){ near=get_lower _ nbrs (edges , simplices [i,1]) for (j in 2:(d[2]-1)){ near =intersect (near ,get _lower _ nbrs (edges , simplices [i,j])) if ( length ( near )== 0)break } for (n in near ) { r=meb(X[c( simplices [i,-d[2]],n) ,])$ radius if(r<= maxscale ){ result=rbind ( result ,c( simplices [i,-d[2]],n,r)) } } } if (! is.null(result))result=matrix(result , ncol=d[2]+1) result } #4. Generar todos los s´ımplices a partir de las aristas fill _ cliques = function (edges,max .dim ){ if(is.null( edges )) simplices = NULL 66 AP ´ ENDICE A. FUNCIONES DE R else{ simplices = list( edges ) if (max.dim ==0)return( simplices ) else{ matrix_ to _ enlarge = edges for (i in 1:max .dim){ new_ simplices = cech . enlarge ( matrix_ to_ enlarge , edges ) if (length(new_ simplices )== 0)break simplices =c(simplices ,list(new_ simplices )) matrix_to_ enlarge = new _ simplices } }} simplices } # Ahora podemos crear una lista de matrices , que almacenan los # v´ertices de cada s´ımplice y la escala a la que entran matrices.list= fill _ cliques ( cech . edges (X , max.scale = maxscale ), max.dim=maxdimension) #Por ´ultimo almacenamos esta informaci´on en un formato compatible #con las funciones de la librer´ıa TDA cmplx =c(1:dim(X)[1], list()) values=rep(0,dim(X)[1]) if (! is.null( matrices . list)){ for (l in 1: length ( matrices . list)){ for (i in 1:dim( matrices . list[[l ]])[ 1]) { s= matrices . list[[l ]][i ,] values=c(values ,s[l+2]) cmplx =c(cmplx,list(s[-(l+2)])) } }} cmplx = cmplx [ order (values)] values=sort(values) result=list( cmplx = cmplx , values = values , increasing =T , coordinates =X) } 67 choose.witness Dado un conjunto de puntos y una proporci´on entre 0 y 1, esta funci´on te separa los puntos en marcadores y testigos con el m´etodo maxmin y te da la matriz de distancias entre un subconjunto y otro. Basado en [29]. choose.witness = function (set,prop){ # Definimos una funci´on auxiliar que nos d´e una matriz de distancias distmat= function (set 1,set 2){ a=dim(set1) b=dim(set2) if(a[2]!= b[2]) { warning(" Points dimension differ .")} m= matrix(0,nrow=a[1], ncol=b[1]) for (i in 1:a[1]){ for (j in 1:b[1]){ m[i,j]= sqrt(sum(( set1[i ,]- set 2[j ,])^ 2)) } } m } #Ya podemos usar el m´etodo maxmin d=dim(set) n= floor (d[1]* prop) r= sample(d[1],1) L= matrix(set[r ,] , nrow=1) W=set[-r ,] dm= distmat (L ,W) nW= dim(W)[1] for (i in 2:n){ mini=c() for (j in 1:nW ){ mini [j ]= min (dm[,j]) } t=which .max( mini ) newrow = distFct ( matrix(W[t,] , nrow=1),W[-t,]) L= rbind (L,W[t,]) W=W[-t,] nW=nW -1 dm=dm[,-t] dm= rbind (dm , newrow ) } rownames( dm )= c() result = list(L=L, W=W, distmatrix = dm) } 68 AP ´ ENDICE A. FUNCIONES DE R witnessFiltration Dada una matriz como las obtenidas con choose.witness, una dimensi´on m´axima para los s´ımplices y una escala m´axima, esta funci´on nos calcula la filtraci´on witness con ν= 0 o ν= 1. Basado en las ideas de [29] y los algoritmos de [74]. witnessFiltration= function ( wlist , maxdimension , maxscale ,v){ if(v!=0& v!=1)warning("v should be either 0 or 1") #En primer lugar definimos una serie de funciones auxiliares #1. Aristas del complejo witness . edges = function(wlist, max.scale, v){ distmat = wlist $ distmatrix d=dim(distmat) if(v==1){ for (j in 1:d[2]){ m=min( distmat [,j ]) distmat [,j]= distmat [,j]-m}} edges =c() for (j in 1:d[2]){ verticesj = which ( distmat [,j] <= max.scale ) if ( length ( verticesj )<2)next for (k in 1:( length ( verticesj )-1)){ for (l in (k+1):( length ( verticesj ))){ edges = rbind (edges ,c( verticesj [l], verticesj [k ]))}}} edges = unique( edges ) if (! is.null( edges )){ values=c() for (n in 1:dim( edges )[ 1]){ values=c(values ,min(apply ( distmat [ edges [n ,] ,] ,2,max )))} edges = cbind ( edges , values ) edges = edges [ order ( values ),] colnames( edges )= c() edges = matrix(edges, ncol=3)} edges } #2. Vecinos menores get_lower _ nbrs = function( edgemat , num ){ edgemat [ edgemat [,1]== num ,2]} #3. Engrandecer ( encontrar s´ımplices de una dimensi´on mayor ) enlarge=function ( simplices , edges ){ d=dim( simplices ) result=c() values=c() for (i in 1:d[1]){ near=get_lower _ nbrs (edges , simplices [i,1]) for (j in 2:(d[2]-1)){ near =intersect (near ,get _lower _ nbrs (edges , simplices [i,j])) if ( length ( near )== 0)break } for (n in near ) {r= simplices [i,d[2]] m2=c() 69 for (j in 1:(d[2]-1)) {m2=c(m2, edges [ edges [,1]== simplices [i ,j ]& edges [,2]==n,3])} r=max(r,max(m2)) result=rbind ( result ,c( simplices [i,-d[2]],n,r ))}} if (! is.null(result)) result=matrix( result , ncol=d[2]+1) result } #4. Generar todos los s´ımplices a partir de las aristas fill _ cliques = function (edges,max .dim ){ if(is.null( edges )) simplices = NULL else{ simplices = list( edges ) if (max.dim ==0)return( simplices ) else{ matrix_ to _ enlarge = edges for (i in 1:max .dim){ new_ simplices = enlarge ( matrix_ to_ enlarge , edges ) if (length(new_ simplices )== 0)break simplices =c(simplices ,list(new_ simplices )) matrix_to_ enlarge = new _ simplices } }} simplices } # Ahora podemos crear una lista de matrices , que almacenan los # v´ertices de cada s´ımplice y la escala a la que entran matrices.list= fill _ cliques ( witness . edges ( wlist , max.scale =maxscale,v), max.dim=maxdimension) #Por ´ultimo almacenamos esta informaci´on en un formato compatible #con las funciones de la librer´ıa TDA cmplx =c(1:dim( wlist $L )[ 1], list()) values=rep(0,dim( wlist $L )[ 1]) if (! is.null( matrices . list)){ for (l in 1: length ( matrices . list)){ for (i in 1:dim( matrices . list[[l ]])[ 1]) { s= matrices . list[[l ]][i ,] values=c(values ,s[l+2]) cmplx =c(cmplx,list(s[-(l+2)])) } }} cmplx = cmplx [ order (values)] values=sort(values) result=list( cmplx = cmplx , values = values , increasing =T , coordinates = wlist $L )}