Entradas

Incrustación estocástica de vecinos (SNE) y su corrección en t-SNE

Autor: Matteo Alberti

 

 

Incrustación estocástica de vecinos (SNE) y su corrección en t-SNE

 

En este tutorial nos disponemos a enfrentar el problema de reducción de la dimensionalidad, con una herramienta significativa:

 

Stochastic Neighbor Embedding o “SNE”, como es normalmente conocido.

 

La definición del problema es:

 

“Tenemos un gran conjunto de datos y queremos encontrar una forma de reducir su dimensionalidad, tanto para el pre-procesamiento como para la clasificación/visualización”.

 

Podrían existir elementos ambiguos. Por lo tanto, la primera pregunta que nos planteamos es: ¿Qué queremos decir con ambigüedad?

 

Ejemplo: Estamos tratando con un problema de clasificación de texto y abordamos la siguiente palabra: bank.

 

Nos concentramos en asociar esta palabra a otras como ‘finance’ o ‘money’ pero podría estar relacionada con ‘river’ (riverbank/orilla del río) al mismo tiempo.

 

Dentro de las herramientas más populares para la reducción de la dimensionalidad, se encuentran:

 

* Multidimensional scaling (MDS).

* PCA (Principal Component Analysis).

* Local Geometry preserved tools (Local Linear Embedding).

 

 

Dado que, cada objeto de un espacio de alta dimensionalidad puede ser asociado con un solo objeto dentro de un espacio de baja dimensionalidad. El SNE trata de mapear objetos de alta dimesionalidad en objetos de poca dimensionalidad, preservando la estructura de parentesco de la vecindad, sin importar el intercambio que resulta de una clasificación incorrecta alrededor de objetos lejanos. El SNE no tiene nada que ver con las medidas de disimilaridad pero está basado en dos conceptos: Entropy y Probability Divergences.

 

La idea principal es centrar una distribución gaussiana para cada valor de entrada (a lo largo del espacio de alta dimensionalidad) a fin de usar su densidad para definir una distribución de probabilidad de todos los vecinos. El objetivo es aproximar esta distribución de probabilidad tanto como sea posible replicando la estructura de parentesco en un espacio de baja dimensionalidad.

 

Presentando el formalismo matemático apropiado:

 

 

 

Dado X1….Xn ∈ RK ,  intentamos construir una función (𝛗) ∈ R2.

Vamos a centrar una distribución gaussiana para cada valor y asociarle una probabilidad.

De esta manera, se espera tener un vistazo en lo que ocurre durante la transición de Rk a R2

 R^k\{x_1, ... , x_k\} \frac{\varphi} {\to} \{Y_1, ..., Y_n\} \in R^2 \\ ... \\ R^k\{x_1, ... , x_k\} \leftrightarrow \{Q_1, ..., Q_n\} \in R^2

 

 

Esta relación ocurre cada que {P1,…, Pk} y {Q1,…, Qn} son lo más similares posible. Con esto, nuestra función objetivo (K-L) es minimizar:

 

C= \sum_{i=1}^k D(P_i||Q_i) = \sum_{i,j} P_{i|j}ln(\frac{P_{j|i}}{Q_{j|i}})

 

Nota explicativa:

Kullback-Leibler Divergence o Cross-Entropy:

Dado dos funciones de probabilidad:

Pi,…, Pk y Qi,…,Qk definimos la entropía cruzada, como:

Con las siguientes propiedades:

1.      D(P||Q) ≥ 0

2.      D(P||Q) = 0                 ↔ Q=P

3.      D(P||Q) ≠ D(Q||P)      No simétrica


 

¿Qué significa esto?

 

  •  Si Pi incrementa, el objeto está más próximo.
  •  Si Qi aumenta mientras Pi disminuye, estamos aproximando dos objetos lejanos en R^k (lo cual no es bueno). Pero es un error despreciable porque nuestro Qi-ésimo elemento está ponderado por Pj|i. De esta forma el SNE preserva la estructura del vecino local sin importar la estructura global.

 

¿Cómo cambia C, mientras nos movemos a lo largo del espacio? ¿Cómo nos podemos mover?

Esto se puede realizar, al estudiar las derivadas parciales:

Primero, recordemos que la densidad de los puntos está dada por ‘sigma’ (Gaussian dilation y Smoothing parameter). Esto define si nuestros valores buscan más lejos o más cerca dentro del espacio.

 

¿Cómo podemos elegir ‘sigma‘?

Tenemos que definir un nuevo objeto llamado perplejidad:

Perplejidad: 2H(Pi). Donde:

Si la perplejidad es constante entonces nos permite cambiar los vecinos al mutar la única sigma.

 

Problema: En caso de alta dimensionalidad, los valores lejanos son acumulados en un ‘bloque periférico’. Este deslizamiento se debe a nuestra disposición de representar de manera adecuada la estructura del vecino.

 

Resolución: Su correción en t-SNE

Por lo tanto, nos disponemos a corregir esta concentración de clasificaciones erróneas. Para hacerlo, operaremos en la estructura general de la distribución gaussiana. La idea principal es muy simple, intenamos elevar las colas de la distribución. Al hacer esto, los puntos originalmente cerrados (cercano a uno) se deslizarán en un factor insignificante. Por otro lado, los valores lejanos experimentarán un deslizamiento considerable, distanciándolos unos de otros.

 

 

 

 

Hagamos uso de él con Python.

Únicamente tenemos que llamar a la biblioteca Scikit, en donde está implementado:

from sklearn.manifold import TSNE

Respecto a t-SNE, vamos a discutir los parámetros más relevantes:

TSNE(n_components=2, perplexity=30.0, early_exaggeration=12.0, learning_rate=200.0, n_iter=1000, n_iter_without_progress=300, min_grad_norm=1e-07, metric=’euclidean’, init=’random’, verbose=0, random_state=None, method=’barnes_hut’, angle=0.5)

 

  • n_components: por defecto (igual a 2), es el número de dimesiones finales.
  • perplejidad: es el parámetro principal. Cuanto mayor sea la perplejidad, mayor será la cantidad de vecinos. Una buena regla es ajustarlo entre 5 y 50.

¡Atención!

 

Si la perplejidad es demasiado pequeña, no capturaremos suficientes vecinos: al contrario, un número excesivo acercará valores demasiado lejanos. Con grandes conjuntos de datos se recomienda mejorar este parámetro.

 

* learning_rate : sugerimos un rango entre 10 y 1000.

 

Pero, ¿cómo deberíamos de inicializarlo de manera correcta en un rango tan grande?

Consideremos dos casos:

 

  • Los datos se concentran en un espacio casi equidistante el uno del otro, esto significa, que emerge un valor demasiado grande.
  • La conenctración fuerte de datos con solo algunos valores atípicos. El parámetro está ajustado muy bajo.

 

 

Ahora, queremos aplicarlo a los datos reales.

El conjunto de datos MNIST, consiste de más de 60 mil imágenes de digitos escritos a mano.

Importamos los paquetes y datos requeridos:

%pylab inline
from sklearn.datasets import fetch_mldata
from sklearn.decomposition import PCA

# Load MNIST dataset
mnist = fetch_mldata("MNIST original")
X, y = mnist.data / 255.0, mnist.target

Hacemos una primera reducción de dimensionalidad con un análisis de componentes principales (PCA): útil cuando existen muchos datos comprimidos, a diferencia de cuando existen datos de sobra, en donde SVD puede ser la mejor alternativa.

indices = arange(X.shape[0])
random.shuffle(indices)
n_train_samples = 5000
X_pca = PCA(n_components=50).fit_transform(X)
X_train = X_pca[indices[:n_train_samples]]
y_train = y[indices[:n_train_samples]]

Empecemos a usar el t-SNE. La perplejidad ha sido definida como 40, después de algunos intentos.

X_train_embedded = TSNE(n_components=2, perplexity=40, verbose=2).fit_transform(X_train)

Representando nuestros resultados:

matplotlib.rc('font', **{'family' : 'sans-serif',
                         'weight' : 'bold',
                         'size'   : 18})
matplotlib.rc('text', **{'usetex' : True})

def plot_mnist(X, y, X_embedded, name, min_dist=10.0):
    fig = figure(figsize=(10, 10))
    ax = axes(frameon=False)
    title("\\textbf{MNIST dataset} -- Two-dimensional "
          "embedding of 70,000 handwritten digits with %s" % name)
    setp(ax, xticks=(), yticks=())
    subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=0.9,
                    wspace=0.0, hspace=0.0)
    scatter(X_embedded[:, 0], X_embedded[:, 1],
            c=y, marker="x")

    if min_dist is not None:
        from matplotlib import offsetbox
        shown_images = np.array([[15., 15.]])
        indices = arange(X_embedded.shape[0])
        random.shuffle(indices)
        for i in indices[:5000]:
            dist = np.sum((X_embedded[i] - shown_images) ** 2, 1)
            if np.min(dist) < min_dist:
                continue
            shown_images = np.r_[shown_images, [X_embedded[i]]]
            imagebox = offsetbox.AnnotationBbox(
                offsetbox.OffsetImage(X[i].reshape(28, 28),
                                      cmap=cm.gray_r), X_embedded[i])
            ax.add_artist(imagebox)
plot_mnist(X[indices[:n_train_samples]], y_train, X_train_embedded,"t-SNE", min_dist=20.0)

 

El t-SNE obtuvo el resultado en 17.51s.

Si queremos comprararlo con otras ténicas de reducción de la dimensionalidad (como se mencionó antes), hay dos cuestiones fundamentales a explicar:

  1. Habilidad discriminatoria efectiva.

Establecer los parámetros de manera adecuada no siempre es fácil. Por lo tanto, aquí hay un enlace que muestra cómo modificar los parámetros principales y mostrar el diseño de datos en tiempo real. Hay que recordar que los datos que se encuentran en el enlace son simulados y provienen de un ambiente controlado.

 

https://distill.pub/2016/misread-tsne/

 

  1. Tiempo de ejecución:

A diferencia de otras herramientas, el t-SNE requiere un número elevado de tiempo y recursos. Como se puede ver, mientras el PCA toma aproximadamente 0.01, el t-SNE es 17s más lento. Tenemos que considerar esto, si  consideramos el lapso de tiempo como un factor crucial.

 

Sin embargo, es posible proporcionar algunas ideas sobre cómo abordar este problema. Aquí hay dos artículos del Departamento de Ciencias de la Computación de la Universidad de California:

 

Fast Optimization for t-SNE, Laurens van der Maaten (University of California)

http://cseweb.ucsd.edu/~lvdmaaten/workshops/nips2010/papers/vandermaaten.pdf

Accelerating t-SNE using Tree-Based Algorithms, Laurens van der Maaten (University of California)

http://lvdmaaten.github.io/publications/papers/JMLR_2014.pdf

 

 

 

References:

http://scikit-learn.org/stable/index.html (Official Website of Scikit libraries)

http://www.cs.toronto.edu/~fritz/absps/sne.pdf (documentation about SNE)

http://www.cs.toronto.edu/~hinton/absps/tsne.pdf (documentation about t-SNE)

 

Métodos para la reducción de la dimensionalidad manifold-based: el caso ISOMAP

Autor: Matteo Alberti


 

En este tutorial veremos una miniserie dedicada a métodos de reducción de dimensionalidad basados en una estructura matemática llamada Manifold.

Entonces, primero que todo comprendamos qué es un Manifold, cuándo es útil y cómo aplicarlo, pero sin entrar en detalles sobre la estructura y las propiedades matemáticas.

“Manifold no es más que un espacio matemático donde viene recreado localmente un espacio euclidiano (de una dimensión específica) pero no a nivel global”

Entonces, ¿cuándo y cómo vamos a aplicar métodos basados en Manifold?

La respuesta es simple: cuando fallan los métodos lineales

Investigemos con más profundidad:

 

Ejemplo:

En los métodos lineales (como PCA, MDS, SVD), nuestros datos pueden redimensionarse solo mediante combinaciones lineales de nuestras características, esto significa que no se encuentran estructuras más complejas y no lineales, y por lo tanto se discriminan.

Veamos un ejemplo gráfico:

Metodos Lineares reduciríamos los datos al perder la estructura geométrica. Como se muestra en la figura, se abordarán los puntos X1 y X2  estos se ven mas cerca.(y, por lo tanto, toda la estructura sería “aplanada”)

Metodos Manifold,

los principales son:

  • ISOMAP
  • Local Linear Embedding (LLE)

 

Esta idea básica nos permite resolver muchos problemas no solo en la fase de preproceso (en la computadora Vision, por ejemplo, nos permite resolver problemas de: estimación de posición, compresión de datos, reducción de ruido en la imagen e interpolación de datos faltantes)

 

 

ISOMAP

Comencemos con la primera herramienta:

El algoritmo consta de tres fases:

(1)

Determinar el número de vecinos (neighbors) en Manifold ( basadonde según las distancias euclidianas). Para cada punto, asignaremos los vecinos (determinados según el número de vecinos fijados como en el algoritmo K-NN o en función de un radio fijo)

 

El radio/ el número de vecinos es constante

(2)

En este paso, calcularemos la geodésica. Por lo tanto, calcularemos la distancia de cada par de puntos pero no trabajará en un espacio euclidiano, corresponderá con el camino más corto entre todos los caminos posibles para conectar dos puntos en el colector. Un ejemplo de uso común puede ser el algoritmo Djkstra, muy común en el cálculo de la navegación por carretera.

(3)

En este último paso aplicaremos un MDS (escalado multidimensional) a la matriz de distancia de las geodésicas para reconstruirlo en un espacio euclidiano que mantendrá la misma estructura geométrica. Los vectores de coordenadas en el espacio euclidiano se eligen para minimizar la siguiente función de costo:

  • Dy corresponde a las matriz de las distancias Euclidee
  • L2 corresponde a la raíz cuadrada de la suma de los elementos
  • τ no es más que una función de conversión (las distancias se reemplazan en productos internos) para utilizar algoritmos de optimización eficientes

 

 

Implementémoslo en Python:

En primer lugar, importamos los paquetes necesarios de  scikit-learn


import numpy as np

import matplotlib.pyplot as plt

from sklearn import manifold

 

entonces:

OUR_DATA= manifold.Isomap(n_neighbors=N, n_components=C, eigen_solver=’auto’, tol=0, max_iter=None, path_method=’auto’, neighbors_algorithm=’auto’, n_jobs=1)

 

Analicemos los parámetros principales:
• n_ neighbors = Número entero, que corresponde al número de vecinos
• n_ components = Integer, generalmente establecido de manera predeterminada en 2. Corresponde al número de coordenadas en el espacio del Colector
• path_method = default ‘auto’, FW para Floyd-Warshall, D para el algoritmo Djkstra establece cómo calcular las distancias en el gráfico

Si queremos ir a imprimir un gráfico con los tiempos de ejecución relativos:

 

from time import time

t0=time()

ISOMAP_DATA= manifold.Isomap(n_neighbors, n_components=2).fit_transform(OUR_DATA)

t1=time()

print("%s: %.2g sec" % ('ISO' t1-t0))

plot.embedding(ISOMAP_DATA, "ISOMAP ottenuto in " % (t1-t0)

plot.show()

Hacemos dos últimas consideraciones de uso práctico:

ISOMAP resulta ser una herramienta de reducción muy efectiva que nos permite eludir muchos problemas relacionados con la linealidad, pero es capaz de preservar una estructura geométrica local naturalmente es a expensas de algunos factores, en particular:
• Sensible a outliers
• Pocos parámetros libres

Esto se debe a que la funcionalidad del algoritmo depende principalmente de la elección de algunos parámetros (la elección del número de vecinos), de hecho, incluso unos pocos valores atípicos pueden, de alguna manera, reunir porciones de datos que de otra manera deberían ser discriminados.
Por este motivo, ISOMAP es muy recomendable cuando ya tenemos algunas ideas sobre la estructura de datos porque, comparado con otros métodos lineales, es computacionalmente más intensivo.

En el próximo tutorial abordaremos LLE, el segundo de nuestros métodos basados en Manifold.

 

Referencias
• http://scikit-learn.org/stable/index.html (Official Website of Scikit libraries)
• http://web.mit.edu/cocosci/Papers/sci_reprint.pdf (documentazione su ISOMAP)

 

Métodos para la reducción de la dimensionalidad manifold-based: el caso ISOMAP

Autor: Matteo Alberti

 

 

 

 

Entre metodologías principales  que hay en la  reducción lineal de dimensionalidad, tenemos las PCA o los componentes principales,estas son ciertamente las principales herramientas del Statistical Machine Learning.

Nos enfocamos muy a menudo en instrumentos que capturan el componente no lineal, el análisis de los componentes principales es el punto de partida para muchos análisis (y como una herramienta de preprocesamiento) y su conocimiento se vuelve imperativo en caso  que las condiciones en la linealidad sean satisfactorios.

En este tutorial vamos a presentar a nivel matemático la formación de los principales componentes de forma clara y concisa, como se implementa en python pero sobre todo su interpretación

“La idea básica es encontrar un sistema de referencia que maximice la varianza de las variables representadas”

Esto se hace dividiendo la varianza total en un número igual de variables de inicio, pero será posible reducir el número en función de la contribución que cada PC proporciona en la construcción de nuestra varianza total.

Nos gustaría recordarle que la aplicación de Análisis de Componentes Principales es útil cuando las variables de inicio no son independientes

 

Dado un conjunto de p variables cuantitativas X1, X2 ,. . . , Xp (variables centradas o estandarizadas) queremos determinar un nuevo conjunto de k variables t.c k≤p indicadas con Ys (s = 1,2, … k) que tienen las siguientes propiedades:

 

sin correlación, reproducir la mayor parte posible de la varianza restante después de la construcción de los primeros s-1 componentes, construido en orden creciente con respecto al porcentaje de varianza total que se reproducen, En promedio, no son nada

Como resultado, la combinación lineal será:

Por lo tanto, debemos encontrar los coeficientes v t.a.c que satisfacen nuestras propiedades.Este es un problema de restricción máxima donde la primera restricción se llama normalización:

Nuestro sistema es así:

Donde la varianza se puede expresar de la siguiente manera:

Y se resuelve mediante el método multiplicador de Lagrange:

Cálculo del gradiente de L1 y su cancelación obteniendo:

Donde   Este sistema admite infinitas soluciones (que respetan la restricción) al disminuir el rango de la matriz de coeficientes del sistema:

que corresponden a λs llamados valores propios de Σ.

Analogamente per la costruzione della seconda PCA (e cosi per tutte le altre) subentra al nostro sistema il Vincolo di Ortogonalità, dato dalla nostra richiesta che le PC siano incorrelate, esprimibile nel seguente modo:

De manera similar, para la construcción del segundo PCA (y para todos los demás), se agrega a nuestro sistema la Restricción de Ortogonalidad, le perdimos nuestra solicitud a los PC que estan sin correlación,podemos dar un ejemplo  de la siguiente manera:

Entonces al establecer la lagrangiana en p + 2 variables obtenemos:

De donde obtenemos el segundo valor propio de Y2 donde recordamos que se aplica la siguiente propiedad:

Propiedades de computador:

a) Cada valor propio de Σ tiene un papel en varianza respectiva del computador.

su valor es positivo

Varianza total

Varianza General

 

CRITERIOS DE SELECCIÓN:

Para la elección del número k (con k <p) de PC que se mantendrá en el análisis no existe un criterio universalmente aceptado y válido. Por lo tanto, es una buena práctica usarlos de forma conjunta y siempre tener en cuenta las necesidades del análisis. Expodremos los principales:

  1. a) Porcentaje acumulado de la varianza total reproducida:

 

                                                                                   

medida absoluta                           Normalización                                     Porcentaje acumulado

 

Establezca un umbral T en el porcentaje acumulado y mantenga el primer kPC en análisis que garantice que excede el umbral

 

b) Screen-Plot

Representa los valores propios con respecto al número de orden del componente

Donde se selecciona el primer k PC basado en la reducción de pendiente (también llamado criterio de codo). En este caso específico, los PC que se mantendrán en análisis serían los dos primeros.

 

c) Criterio Kaiser

También conocido como criterio de valor proprio mayor que 1 (válido solo para variables estandarizadas)

 

 

Vamos a implementar:

Importamos los paquetes necesarios de scikit-learn

import numpy as np

from sklearn.decomposition import PCA

La clase tiene los siguientes atributos:

Sklearn.decomposition.PCA(n_components=None, copy=True, whiten=False, svd_solver=’auto’, tol=0.0, iterated_power=’auto’, random_state=None)

Comentemos los principales parámetros:

  • n_components = número de componentes a analizar. Aconsejamos en una fase de elección no establecer este valor y evaluar en función de los criterios anteriores.
  • svd_solver = nos da algunas de las alternativas principales (arpack, randomized ..) queremos recordar que PCA no admite datos dispersos (para lo cual deberá cargar TruncatedSVD)

Vamos a probarlo en nuevos datos reales, por ejemplo, en los datos de Wine que se pueden importar a través del comando

from sklearn import datasets

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

from sklearn import decomposition

En este punto vamos a implementar PCA y hacer un plot:

np.random.seed(123)
wine = datasets.load_wine()

X = wine.data

y = wine.target



fig = plt.figure(1, figsize=(5, 4))

plt.clf()

ax = Axes3D(fig, rect=[1, 0, 1, 0.9], elev=30, azim=222)



plt.cla()

pca = decomposition.PCA(n_components=None)

pca.fit(X)

X = pca.transform(X)





for name, label in [('Setosa', 0), ('Versicolour', 1), ('Virginica', 2)]:

ax.text3D(X[y == label, 0].mean(),

X[y == label, 1].mean() + 1.5,

X[y == label, 2].mean(), name,

bbox=dict(alpha=.5, edgecolor='b', facecolor='w'))

# Reorder the labels to have colors matching the cluster results

y = np.choose(y, [1, 2, 0]).astype(np.float)

ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=y, cmap=plt.cm.spectral,

edgecolor='r')



ax.w_xaxis.set_ticklabels([])

ax.w_yaxis.set_ticklabels([])

ax.w_zaxis.set_ticklabels([])



plt.show()

Este será el gráfico que obtendremos:

Por lo tanto, vamos a seleccionar solo las primeras 3 variables de referencia