Rotación de componentes PCA para igualar la varianza en cada componente

Question

Estoy tratando de reducir la dimensionalidad y el ruido de un conjunto de datos mediante la realización de PCA en el conjunto de datos y se deshace de los últimos Equipos. Después de eso, quiero usar algunos algoritmos de aprendizaje automático, en el resto de PCs, y por lo tanto quiero normalizar los datos mediante la nivelación de la varianza de la Pc para hacer el trabajo de los algoritmos mejor.

Una forma sencilla es simplemente normalizar la varianza de los valores unitarios. Sin embargo, el primer PC contiene más de la varianza del conjunto de datos original de la siguiente, y todavía quiero darle más "peso". Por lo tanto, me pregunto: ¿hay una manera sencilla, recientemente desvinculado de su varianza y compartirlo con los equipos con menos variaciones?

Otra forma es el mapa de la Pc de nuevo a la original espacio de características, pero en caso de que la dimensionalidad también aumentaría el valor original.

Supongo que es mejor de mantener la resultante de las columnas ortogonales, pero no es necesario en este momento.

Answer 1

No es completamente claro para mí que lo que están pidiendo es lo que realmente se necesita: una común preprocesamiento paso en el aprendizaje de máquina es la reducción de dimensionalidad + blanqueamiento, lo que significa hacer de la PCA y la estandarización de los componentes, nada más. Pero yo, sin embargo, se centran en su pregunta, ya que está formulado, porque es más interesante.

---

Deje $\mathbf X$ ser el centrado $n\times d$ matriz de datos con puntos de datos en filas y las variables en las columnas. PCA cantidades a la descomposición de valor singular $$\mathbf X = \mathbf{USV}^\top \approx \mathbf U_k \mathbf S_k \mathbf V_k^\top,$$ where to perform the dimensionality reduction we keep only $k$ components. An orthogonal "factor rotation" of these components implies choosing an orthogonal $k \times k$ matrix $\mathbf R$ and plugging it into the decomposition: $$\mathbf X \approx \mathbf U_k \mathbf S_k \mathbf V_k^\top = \mathbf U_k \mathbf {RR}^\top \mathbf S_k \mathbf V_k^\top = \underbrace{\sqrt{n-1}\mathbf U_k^\phantom\top \mathbf {R}}_{\substack{\text{Rotated}\\\text{standardized scores}}} \cdot \underbrace{\mathbf R^\top \mathbf S_k \mathbf V_k^\top/\sqrt{n-1}}_{\text{Rotated loadings}^\top}.$$ Here $\sqrt{n-1}\mathbf U_k \mathbf R$ are rotated standardized components and the second term represents rotated loadings transposed. The variance of each component after rotation is given by the sum of squares of the corresponding loading vector; before rotation it is simply $s_i^2/(n-1)$. After rotation it is something else.

Now we are ready to formulate the problem in mathematical terms: given unrotated loadings $\mathbf L = \mathbf V_k \mathbf S_k / \sqrt{n-1}$, find rotation matrix $\mathbf R$ such that the rotated loadings, $\mathbf L \mathbf R$, has equal sum of squares in each column.

Let's solve it. Column sums of squares after rotation are equal to the diagonal elements of $$(\mathbf {LR})^\top \mathbf{LR} = \mathbf R^\top \frac{\mathbf S^2}{n-1} \mathbf R.$$ This makes sense: rotation simply redistributes the variances of components, which are originally given by $s_i^2/(n-1)$, between them, according to this formula. We need to redistribute them such they all become equal to their average value $\mu$.

I don't think there is a closed form solution to this, and in fact there are many different solutions. But a solution can be easily built in a sequential fashion:

1. Take the first component and the $k$-th component. The first one has variance $\sigma_\text{max}>\mu$ and the last one has the variance $\sigma_\text{min}<\mu$.
2. Girar sólo estos dos tales que la varianza de la primera es igual a $\mu$. Matriz de rotación en 2D sólo depende de un parámetro $\theta$ y es fácil escribir la ecuación y calcular la necesaria $\theta$. De hecho, $$\mathbf R_\text{2D} = \left(\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin\theta & \cos \theta\end{array}\right)$$ and after transformation the first PC will get variance $$\cos^2\theta \cdot \sigma_\text{max} + \sin^2\theta \cdot \sigma_\text{min} = \cos^2\theta \cdot \sigma_\text{max} + (1-\cos^2\theta)\cdot \sigma_\text{min} =\mu,$$ from which we immediately obtain $$\cos^2\theta = \frac{\mu-\sigma_\text{min}}{\sigma_\text{max}-\sigma_\text{min}}.$$
3. El primer componente se hace ahora, se ha varianza $\mu$.
4. Proceder a la siguiente pareja, tomando el componente con la mayor varianza y la uno con la menor varianza. Goto #2.

Esto va a redistribuir todos los varianzas iguales por una secuencia de $(k-1)$ 2D rotaciones. La multiplicación de todas estas matrices de rotación juntos producirán un total $\mathbf R$.

---

Ejemplo

Considere el siguiente $\mathbf S^2/(n-1)$ matriz: $$\left(\begin{array}{cccc}10&0&0&0\\0&6&0&0\\0&0&3&0\\0&0&0&1\end{array}\right).$$ The mean variance is $5$. My algorithm will proceed as follows:

1. Step 1: rotate PC1 and PC4 so that PC1 gets variance $5$. As a result, PC4 gets variance $1+(10-5)=6$.

2. Paso 2: gire el PC2 (nueva máxima varianza) y PC3, de modo que la PC2 obtiene la varianza $5$. Como resultado, PC3 obtiene la varianza $3+(6-5)=4$.

3. Paso 3: gire la PC4 (nueva máxima varianza) y PC3 para que PC4 obtiene la varianza $5$. Como resultado, PC3 obtiene la varianza $4+(6-1)=5$.

4. Hecho.

Escribí la secuencia de comandos de Matlab que implementa este algoritmo (ver más abajo). Para esta matriz de entrada, la secuencia de ángulos de rotación es:

48.1897   35.2644   45.0000

Componente de varianzas después de cada paso (en filas):

10     6     3     1
 5     6     3     6
 5     5     4     6
 5     5     5     5

El final de la matriz de rotación (producto de tres matrices de rotación 2D):

 0.6667         0    0.5270    0.5270
      0    0.8165    0.4082   -0.4082
      0   -0.5774    0.5774   -0.5774
-0.7454         0    0.4714    0.4714

Y el final de la $(\mathbf{LR})^\top \mathbf{LR}$ de la matriz es:

5.0000         0    3.1623    3.1623
     0    5.0000    1.0000   -1.0000
3.1623    1.0000    5.0000    1.0000
3.1623   -1.0000    1.0000    5.0000

Aquí está el código:

S = diag([10 6 3 1]);
mu = mean(diag(S));
R = eye(size(S));

vars(1,:) = diag(S);
Supdated = S;

for i = 1:size(S,1)-1
    [~, maxV] = max(diag(Supdated));
    [~, minV] = min(diag(Supdated));

    w = (mu-Supdated(minV,minV))/(Supdated(maxV,maxV)-Supdated(minV,minV));
    cosTheta = sqrt(w);
    sinTheta = sqrt(1-w);

    R2d = eye(size(S));
    R2d([maxV minV], [maxV minV]) = [cosTheta sinTheta; -sinTheta cosTheta];
    R = R * R2d;

    Supdated = transpose(R2d) * Supdated * R2d;    

    vars(i+1,:) = diag(Supdated);
    angles(i) = acosd(cosTheta);
end

angles                %// sequence of 2d rotation angles
round(vars)           %// component variances on each step
R                     %// final rotation matrix
transpose(R)*S*R      %// final S matrix

Aquí está el código en Python proporcionada por @feilong:

def amoeba_rotation(s2):
    """
    Parameters
    ----------
    s2 : array
        The diagonal of the matrix S^2.

    Returns
    -------
    R : array
        The rotation matrix R.

    Examples
    --------
    >>> amoeba_rotation(np.array([10, 6, 3, 1]))
    [[ 0.66666667  0.          0.52704628  0.52704628]
     [ 0.          0.81649658  0.40824829 -0.40824829]
     [ 0.         -0.57735027  0.57735027 -0.57735027]
     [-0.74535599  0.          0.47140452  0.47140452]]

    http://stats.stackexchange.com/a/177555/87414
    """
    n = len(s2)
    mu = s2.mean()
    R = np.eye(n)
    for i in range(n-1):
        max_v, min_v = np.argmax(s2), np.argmin(s2)
        w = (mu - s2[min_v]) / (s2[max_v] - s2[min_v])
        cos_theta, sin_theta = np.sqrt(w), np.sqrt(1-w)
        R[:, [max_v, min_v]] = np.dot(
            R[:, [max_v, min_v]],
            np.array([[cos_theta, sin_theta], [-sin_theta, cos_theta]]))
        s2[[max_v, min_v]] = [mu, s2[max_v] + s2[min_v] - mu]
    return R

---

Tenga en cuenta que este problema es totalmente equivalente a la siguiente: dado $k$ variables correlacionadas con las variaciones $\sigma_i^2$, encontramos una rotación (es decir, una nueva base ortogonal) que producirá $k$ variables con igualdad de varianzas (pero, por supuesto, no correlacionadas).

Answer 2

En su perspicaz y una respuesta exhaustiva @ameba ha demostrado - como parte de la respuesta -, ¿cómo se puede girar dos variables correlacionadas (como el de componentes principales, por ejemplo) para lograr la quería desviaciones para ellos (aunque a costa de perder uncorrelatedness, por supuesto). Vamos ortogonal variables $X$ $Y$ tienen variaciones $\sigma^2_{max}$ (más grande) y $\sigma^2_{min}$ (un menor), respectivamente. Hacerlos girar de manera que $X$ obtendrá arbitraria, la disminución de la varianza $\mu^2$ (mientras $Y$, en consecuencia, será de la varianza $\sigma^2_{max}+\sigma^2_{min}-\mu^2$).

@ameba muestra la fórmula a partir de la cual podemos calcular el ángulo de rotación, $\cos\theta$:

$$\mu^2 = \cos^2\theta (\sigma^2_{max}) + \sin^2\theta (\sigma^2_{min})$$

pero no ha demostrado donde esta ecuación proviene, probablemente pensando que es obvio, sin explicación. Obvio o no, creo que vale la pena la aclaración - de alguna manera. Mi respuesta se presenta de una manera.

Y así, tenemos un elipsoidal, centrada en la nube de datos en el espacio de la no correlación de las variables de $X$$Y$. Tenemos que rotar los ejes ángulo de $\theta$. Un punto de datos en la nube (como el que se muestra como el punto verde en la imagen) con $X$ coordinar $x$ tendrá este coordinar como $x^*$ después de la rotación.

Observar que la proyección de la coordenada $x$ muesca en el eje rotado $X^*$ está dado por $x'=x\cos\theta$ (cathetus como la hipotenusa y el ángulo entre ellos). Observar también que $x^*$ es de menos de $x'$ por el corte de la longitud de la $x'-x^*$ computables a partir de coordinar $y$: $y\sin\theta$ (otro cathetus y la hipotenusa). Y así,

$$x^* = x' - (x'-x^*) = x\cos\theta-y\sin\theta$$

Sabemos (véase el principio) las desviaciones (o sumas de cuadrados) de las dos variables y la varianza (suma de cuadrados) $\mu^2$$X^*$. A continuación, de la siguiente manera:

$$\mu^2=\sum x^{*2} = \sum(x\cos\theta-y\sin\theta)^2 = \sum(x^2\cos^2\theta+y^2\sin^2\theta-2xy\cos\theta\sin\theta) = \cos^2\theta\sum x^2 + \sin^2\theta\sum y^2 - \underbrace{ 2\cos\theta\sin\theta\sum xy}_{\text{=0 (X and Y are uncorrelated)}} = \cos^2\theta (\sigma^2_{max}) + \sin^2\theta (\sigma^2_{min})$$

A partir de la cual se puede estimar $\cos\theta$ @ameba ha demostrado, y realizar el giro.

Answer 3

Si he de interpretar las cosas correctamente, quiere decir que el primer componente principal (autovalor) explica la mayor parte de la varianza de los datos. Esto puede suceder cuando el método de compresión es lineal. Sin embargo, no puede ser no-lineal de las dependencias en su función de espacio.

TL/DR: PCA es un método lineal. Uso Autoencoders (no lineal pca) para la reducción de dimensionalidad. Si la máquina de aprendizaje es parte de aprendizaje supervisado, a continuación, simplemente monitor de su función de pérdida, mientras que el ajuste de la (hiper)parámetros para la autoencoder. De esta manera usted va a terminar con una mejor versión comprimida de los datos originales.

He aquí un scikit ejemplo donde se hace el grid de búsqueda para encontrar el número óptimo de componentes principales para mantener (hyper-parámetro) uso de PCA. Finalmente se aplica la Regresión Logística en el menor espacio de dimensiones: http://scikit-learn.org/stable/auto_examples/plot_digits_pipe.html#example-plot-digits-pipe-py

Protip: Autoencoders no tiene una solución de forma cerrada (afaik) así que si el contexto es de transmisión de datos, esto significa que usted puede actualizar continuamente su autoencoder (representación comprimida) y por lo tanto puede compensar para cosas tales como el concepto de la deriva. Con pca tendrá que re-entrenar a modo de proceso por lotes de tiempo en tiempo, a medida que nuevos datos.

Como para dar algunas de las características más "peso", ver regularización ( me gustaría empezar de normas https://en.wikipedia.org/wiki/Norm_(matemáticas) ). Usted también puede estar sorprendido de lo similar que la regresión logística es el tipo perceptrón.