¿Por qué el Lik scikit-learn de Python no funciona correctamente y cómo calcula LDA a través de SVD?

Estaba utilizando el Análisis discriminante lineal (LDA) de la scikit-learnbiblioteca de aprendizaje automático (Python) para la reducción de la dimensionalidad y tenía un poco de curiosidad por los resultados. Ahora me pregunto qué scikit-learnestá haciendo la LDA para que los resultados se vean diferentes de, por ejemplo, un enfoque manual o una LDA realizada en R. Sería genial si alguien pudiera darme algunas ideas aquí.

Lo que básicamente es más preocupante es que scikit-plotmuestra una correlación entre las dos variables donde debería haber una correlación 0.

Para una prueba, utilicé el conjunto de datos de Iris y los primeros 2 discriminantes lineales se veían así:

IMG-1. LDA a través de scikit-learn

Esto es básicamente consistente con los resultados que encontré en la documentación de scikit-learn aquí.

Ahora, revisé el LDA paso a paso y obtuve una proyección diferente. Intenté diferentes enfoques para averiguar qué estaba pasando:

IMG-2. LDA en datos sin procesar (sin centrado, sin estandarización)

Y aquí estaría el enfoque paso a paso si estandarizara (normalización del puntaje z; varianza de la unidad) los datos primero. Hice lo mismo solo con centrado en la media, lo que debería conducir a la misma imagen de proyección relativa (y que de hecho lo hizo).

IMG-3. LDA paso a paso después del centrado medio o estandarización

IMG-4. LDA en R (configuración predeterminada)

LDA en IMG-3, donde centré los datos (que sería el enfoque preferido) también se ve exactamente igual al que encontré en una publicación de alguien que hizo el LDA en R

Código de referencia

No quería pegar todo el código aquí, pero lo he subido como un cuaderno de IPython desglosado en varios pasos que utilicé (ver más abajo) para la proyección LDA.

1. Paso 1: Calcular los vectores medios d-dimensionales $$\mathbf m_i = \frac{1}{n_i} \sum\limits_{\mathbf x \in D_i}^n \; \mathbf x_k$$

2. Paso 2: Calcular las matrices de dispersión

   2.1 La matriz de dispersión dentro de clase se calcula mediante la siguiente ecuación:$S_W$
   

   $$S_W = \sum\limits_{i=1}^{c} S_i = \sum\limits_{i=1}^{c} \sum\limits_{\mathbf x \in D_i}^n (\mathbf x - \mathbf m_i)\;(\mathbf x - \mathbf m_i)^T$$

   2.2 La matriz de dispersión entre clases se calcula mediante la siguiente ecuación: donde es la media general.$S_B$
   

   $$S_B = \sum\limits_{i=1}^{c} n_i (\mathbf m_i - \mathbf m) (\mathbf m_i - \mathbf m)^T$$ $\mathbf m$

3. Paso 3. Resolver el problema de valor propio generalizado para la matriz$S_{W}^{-1}S_B$

   3.1. Ordenar los vectores propios disminuyendo los valores propios

   3.2. Elección de k vectores propios con los valores propios más grandes. Combinando los dos vectores propios con los valores propios más altos para construir nuestra matriz de vectores propios -dimensional$d \times k$$\mathbf W$

4. Paso 5: Transformando las muestras en el nuevo subespacio

   $$\mathbf y = \mathbf W^T \times \mathbf x.$$

Actualización: Gracias a esta discusión, scikit-learnse actualizó y ahora funciona correctamente. Su código fuente LDA se puede encontrar aquí . El problema original se debió a un error menor (vea esta discusión de github ) y mi respuesta en realidad no apuntaba correctamente (disculpas por cualquier confusión causada). Como todo eso ya no importa (el error está solucionado), edité mi respuesta para centrarme en cómo se puede resolver LDA a través de SVD, que es el algoritmo predeterminado en scikit-learn.

Después de definir las matrices de dispersión dentro y entre las clases y , el cálculo estándar de LDA, como se señala en su pregunta, es tomar vectores propios de como ejes discriminantes ( ver, por ejemplo, aquí ). Sin embargo, los mismos ejes se pueden calcular de una manera ligeramente diferente, explotando una matriz de blanqueamiento:$\boldsymbol \Sigma_W$$\boldsymbol \Sigma_B$$\boldsymbol \Sigma_W^{-1} \boldsymbol \Sigma_B$

1. Calcule . Esta es una transformación de blanqueamiento con respecto a la covarianza agrupada dentro de la clase (vea mi respuesta vinculada para más detalles).$\boldsymbol \Sigma_W^{-1/2}$

   Tenga en cuenta que si tiene descomposición propia , entonces . Tenga en cuenta también que uno calcula lo mismo haciendo SVD de datos agrupados dentro de la clase: .$\boldsymbol \Sigma_W = \mathbf{U}\mathbf{S}\mathbf{U}^\top$$\boldsymbol \Sigma_W^{-1/2}=\mathbf{U}\mathbf{S}^{-1/2}\mathbf{U}^\top$$\mathbf{X}_W = \mathbf{U} \mathbf{L} \mathbf{V}^\top \Rightarrow \boldsymbol\Sigma_W^{-1/2}=\mathbf{U}\mathbf{L}^{-1}\mathbf{U}^\top$

2. Encuentre vectores propios de , llamémoslos .$\boldsymbol \Sigma_W^{-1/2} \boldsymbol \Sigma_B \boldsymbol \Sigma_W^{-1/2}$$\mathbf{A}^*$

   Una vez más, tenga en cuenta que se puede calcular haciendo SVD de datos entre clases , transformado con , es decir, datos entre clases blanqueados con respecto a dentro de la clase covarianza$\mathbf{X}_B$$\boldsymbol \Sigma_W^{-1/2}$

3. Los ejes discriminantes estarán dados por , es decir, por los ejes principales de los datos transformados , transformados nuevamente .$\mathbf A$$\boldsymbol \Sigma_W^{-1/2} \mathbf{A}^*$

   De hecho, si es un vector propio de la matriz anterior, entonces y multiplicando desde la izquierda por y definiendo , obtenemos inmediatamente :$\mathbf a^*$

   $$\boldsymbol \Sigma_W^{-1/2} \boldsymbol \Sigma_B \boldsymbol \Sigma_W^{-1/2}\mathbf a^* = \lambda \mathbf a^*,$$ $\boldsymbol \Sigma_W^{-1/2}$$\mathbf a = \boldsymbol \Sigma_W^{-1/2}\mathbf a^*$ $$\boldsymbol \Sigma_W^{-1} \boldsymbol \Sigma_B \mathbf a = \lambda \mathbf a.$$

En resumen, LDA es equivalente a blanquear la matriz de los medios de clase con respecto a la covarianza dentro de la clase, hacer PCA en los medios de clase y volver a transformar los ejes principales resultantes en el espacio original (sin blanquear).

Esto se señala, por ejemplo, en Los elementos del aprendizaje estadístico , sección 4.3.3. En scikit-learnesta es la forma predeterminada para calcular LDA porque SVD de una matriz de datos es numéricamente más estable que eigen-descomposición de su matriz de covarianza.

Tenga en cuenta que uno puede usar cualquier transformación de blanqueamiento en lugar de y todo seguirá funcionando exactamente igual. En se usa (en lugar de ), y funciona bien (al contrario de lo que se escribió originalmente en mi respuesta). L - 1 U ⊤ U L - 1 U$\boldsymbol \Sigma_W^{-1/2}$scikit-learn $\mathbf{L}^{-1}\mathbf{U}^\top$$\mathbf{U}\mathbf{L}^{-1}\mathbf{U}^\top$

Solo para cerrar esta pregunta, el problema discutido con la LDA se ha solucionado en scikit-learn 0.15.2 .