Faites pivoter les composants PCA pour égaliser la variance de chaque composant

9

J'essaie de réduire la dimensionnalité et le bruit d'un ensemble de données en effectuant l'ACP sur l'ensemble de données et en jetant les derniers PC. Après cela, je veux utiliser certains algorithmes d'apprentissage automatique sur les PC restants, et donc je veux normaliser les données en égalisant la variance des PC pour améliorer le fonctionnement des algorithmes.

Une façon simple consiste à normaliser simplement la variance par rapport aux valeurs unitaires. Cependant, le premier PC contient plus de variance par rapport à l'ensemble de données d'origine que les suivants, et je veux toujours lui donner plus de "poids". Je me demandais donc: existe-t-il un moyen simple de simplement diviser sa variance et de la partager avec les PC avec moins de variances?

Une autre façon consiste à mapper les PC à l'espace de fonctionnalité d'origine, mais dans ce cas, la dimensionnalité augmenterait également à la valeur d'origine.

Je suppose qu'il vaut mieux garder les colonnes résultantes orthogonales, mais ce n'est pas nécessaire pour le moment.

feilong
la source
1
Non ... varimax maximise la somme des variances au carré des chargements, donc il essaie de les rendre aussi inégaux que possible. Aussi, pourquoi voudriez-vous égaliser les composants? Le but est de capturer autant de variations que possible dans le moins de composants possible.
2
La standardisation des scores des composants en variances unitaires ne vous convient-elle pas? Pourquoi alors? Quel type de résultat souhaitez-vous - les colonnes résultantes devraient-elles être non corrélées en plus des variances égales?
ttnphns
2
D'après votre description, il semble que vous souhaitiez simplement "sphérer" les données (de dimensionnalité réduite). Elle est souvent effectuée comme une étape de prétraitement dans l'apprentissage automatique. Pour y parvenir, il vous suffit d'effectuer l'ACP, de choisir certains composants et de les normaliser. Je suppose qu'il est possible de trouver une rotation orthogonale (telle que varimax) qui fait tourner les composants normalisés de telle sorte qu'ils restent non corrélés mais expliquent exactement la même quantité de variance; c'est une question intéressante, je dois y réfléchir. Mais je n'ai jamais vu cela se faire, certainement pas dans l'apprentissage automatique.
amoeba
2
Soit dit en passant, quels sont "certains algorithmes d'apprentissage automatique" que vous souhaitez appliquer après PCA? Cela pourrait être pertinent.
amoeba
1
Notez que si vous faites pivoter vos PC standardisés, les distances ne changeront pas du tout! Donc, cela ne devrait vraiment pas avoir d'importance pour tout algorithme basé sur la distance ultérieur.
amibe du

Réponses:

10

Il n'est pas tout à fait clair pour moi que ce que vous demandez est ce dont vous avez vraiment besoin: une étape de prétraitement courante dans l'apprentissage automatique est la réduction de dimensionnalité + le blanchiment, ce qui signifie faire de l'ACP et standardiser les composants, rien d'autre. Mais je vais néanmoins me concentrer sur votre question telle qu'elle est formulée, car elle est plus intéressante.


Soit la matrice de données n × d centrée avec des points de données en lignes et des variables en colonnes. PCA équivaut à une décomposition en valeurs singulières X = U S VU k S k V k , où pour effectuer la réduction de dimensionnalité, nous ne gardons que k composantes. Une «rotation factorielle» orthogonale de ces composantes implique de choisir une matrice orthogonale k × k R et de la brancher dans la décomposition: XU k S k VXn×

X=USVUkSkVk,
kk×kRIci
XUkSkVk=UkRRSkVk=n-1UkRTournéscores normalisésRSkVk/n-1Charges tournées.
sont des composants normalisés tournés et le deuxième terme représente des chargements tournés transposés. La variance de chaque composante après rotation est donnée par la somme des carrés du vecteur de chargement correspondant; avant la rotation, il s'agit simplement des 2 i /(n-1). Après la rotation, c'est autre chose.n1UkRsi2/(n1)

Nous sommes maintenant prêts à formuler le problème en termes mathématiques: compte tenu des chargements non rotatifs , trouver la matrice de rotationRtelle que les chargements tournés,LR, ont une somme égale de carrés dans chaque colonne.L=VkSk/n1RLR

Résolvons-le. Les sommes des colonnes des carrés après rotation sont égales aux éléments diagonaux de Cela a du sens: la rotation redistribue simplement les variances des composantes, qui sont à l'origine données pars 2 i /(n-1), entre elles, selon cette formule. Nous devons les redistribuer de sorte qu'ils deviennent tous égaux à leur valeur moyenneμ.

(LR)LR=RS2n1R.
si2/(n1)μ

Je ne pense pas qu'il existe une solution de forme fermée à cela, et en fait, il existe de nombreuses solutions différentes. Mais une solution peut être facilement construite de manière séquentielle:

  1. Prenez le premier composant et le composant -ème. Le premier a la variance σ max > μ et le dernier a la variance σ min < μ .kσmax>μσmin<μ
  2. Tournez seulement ces deux pour que la variance du premier devienne égale à . La matrice de rotation en 2D ne dépend que d'un seul paramètre θ et il est facile d'écrire l'équation et de calculer le θ nécessaire . En effet, R 2D = ( cos θ sin θ - sin θ cos θ ) et après transformation le premier PC obtiendra la variance cos 2 θ σ max + sin 2 θ σ min = cos 2 θ σμθθ
    R2D=(cosθsinθsinθcosθ)
    partir duquel on obtient immédiatement cos 2 θ = μ - σ min
    cos2θσmax+péché2θσmin=cos2θσmax+(1-cos2θ)σmin=μ,
    cos2θ=μ-σminσmax-σmin.
  3. Le premier composant est maintenant terminé, il a la variance .μ
  4. Passez à la paire suivante, en prenant le composant avec la plus grande variance et celui avec la plus petite variance. Allez à # 2.

(k-1)R


Exemple

S2/(n-1)

(dix000060000300001).
5
  1. 51+(dix-5)=6

  2. 53+(6-5)=4

  3. 54+(6-1)=5

  4. Terminé.

J'ai écrit le script Matlab qui implémente cet algorithme (voir ci-dessous). Pour cette matrice d'entrée, la séquence des angles de rotation est:

48.1897   35.2644   45.0000

Variations des composants après chaque étape (en lignes):

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

La matrice de rotation finale (produit de trois matrices de rotation 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

(LR)LR

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

Voici le code:

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

Voici le code en Python fourni par @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

kσje2k

amibe
la source
Je suppose que pour deux paires de composants (leurs scores), l'angle de rotation serait de 45 degrés pour égaliser leurs variances. Cependant, je ne peux pas imaginer comment effectuer l'ensemble de la tâche avec 3+ composants en binôme.
ttnphns
1
@feilong, je pense que l'égalisation de la variance d'une paire de composants à la fois est un algorithme très sous-optimal. Ce que j'ai suggéré, c'est de choisir les rotations de telle sorte que la variance d'une composante devienne exactement égale à la variance moyenne globale. Ensuite, ce composant est "terminé", et on peut s'occuper du reste. Ceci est garanti pour égaliser toutes les variances dans un nombre fini d'étapes. Voir mon commentaire précédent pour un exemple.
amoeba
1
@amoeba Vous avez raison, c'est une meilleure solution et vous devriez terminer par n-1 étapes.
feilong
1
@amoeba J'ai ajouté mon implémentation minimale en utilisant Python. J'ai modifié la partie en multipliant la matrice entière, car cela peut prendre du temps pour les grandes matrices.
feilong
1
@amoeba Spécifiquement pour les composants principaux, il est possible de gagner plus de temps en supprimant la pièce recherchant le maximum et le minimum. Nous pouvons simplement faire pivoter les 1er et 2e composants (pour que le 1er composant ait une variance moyenne), puis les 2e et 3e, et ainsi de suite. Nous devons juste nous assurer que la variance totale de chaque paire est supérieure à mu.
feilong
2

XOuiσmuneX2σmjen2Xμ2OuiσmuneX2+σmjen2-μ2

cosθ

μ2=cos2θ(σmuneX2)+péché2θ(σmjen2)

mais n'a pas démontré d'où vient cette équation; pensant probablement que c'est évident sans explication. Évident ou non, je crois que cela vaut la peine d'être élucidé - d'une certaine manière. Ma réponse présente une façon.

XOuiθXXX

illustration de la rotation

X XX=XcosθXXX-Xyypéchéθ

X=X-(X-X)=Xcosθ-ypéchéθ

μ2X

μ2=X2=(Xcosθ-ypéchéθ)2=(X2cos2θ+y2péché2θ-2Xycosθpéchéθ)=cos2θX2+péché2θy2-2cosθpéchéθXy= 0 (X et Y ne sont pas corrélés)=cos2θ(σmuneX2)+péché2θ(σmjen2)

cosθ

ttnphns
la source
2
(cosθpéchéθ-péchéθcosθ)(σmax200σmin2)(cosθpéchéθ-péchéθcosθ),
amoeba
Et je pense que votre explication géométrique et votre calcul "direct" (sans matrices) sont plus faciles à comprendre et très utiles pour développer les bonnes intuitions.
amoeba
0

Si j'interprète les choses correctement, vous voulez dire que la première composante principale (valeur propre) explique la majeure partie de la variance dans les données. Cela peut se produire lorsque votre méthode de compression est linéaire. Cependant, il peut y avoir des dépendances non linéaires dans votre espace d'entités.

TL / DR: PCA est une méthode linéaire. Utilisez les encodeurs automatiques (pca non linéaire) pour réduire la dimensionnalité. Si la partie d'apprentissage automatique est un apprentissage supervisé, surveillez simplement votre fonction de perte tout en ajustant les paramètres (hyper) pour l'encodeur automatique. De cette façon, vous vous retrouverez avec une bien meilleure version compressée de vos données d'origine.

Voici un exemple de scikit où ils effectuent une recherche dans la grille pour trouver le nombre optimal de composants principaux à conserver (hyper-paramètre) à l'aide de PCA. Enfin, ils appliquent une régression logistique sur l'espace dimensionnel inférieur: http://scikit-learn.org/stable/auto_examples/plot_digits_pipe.html#example-plot-digits-pipe-py

Protip: Les encodeurs automatiques n'ont pas de solution de formulaire fermé (afaik), donc si votre contexte est en train de diffuser des données, cela signifie que vous pouvez mettre à jour en continu votre encodeur automatique (représentation compressée) et ainsi compenser des choses telles que la dérive du concept. Avec pca, vous devez réentraîner le mode batch de temps en temps à mesure que de nouvelles données entrent.

Quant à donner à certaines fonctionnalités plus de "poids", voir la régularisation (je commencerais par les normes https://en.wikipedia.org/wiki/Norm_(mathematics) ). Vous pourriez également être surpris de voir à quel point la régression logistique est similaire à celle du perceptron.

shuriken x bleu
la source
Je ne vois pas comment cela répond à la question du PO; votre réponse ne semble pas du tout liée à la question.
amibe du
Je me demandais donc: existe-t-il un moyen simple de simplement diviser sa variance et de la partager avec les PC avec moins de variances? OP veut réduire la dimensionnalité. J'ai proposé une alternative pour résoudre son problème, car en fin de compte, ce que OP souhaite ne garantit pas de meilleures performances à moins que les performances ne soient mesurées. Travailler dans des espaces hilbert / espaces normés ne garantit pas de meilleurs résultats. La mesure des performances conduit à de meilleurs résultats.
shuriken x blue