Comment effectuer une validation croisée pour PCA afin de déterminer le nombre de composants principaux?

13

J'essaie d'écrire ma propre fonction pour l'analyse des composants principaux, PCA (bien sûr, il y a déjà beaucoup écrit mais je suis juste intéressé à implémenter des choses par moi-même). Le principal problème que j'ai rencontré est l'étape de validation croisée et le calcul de la somme prédite des carrés (PRESS). Peu importe la validation croisée que j'utilise, c'est une question principalement sur la théorie derrière, mais pensez à la validation croisée avec absence de contact (LOOCV). D'après la théorie, j'ai découvert que pour effectuer LOOCV, vous devez:

  1. supprimer un objet
  2. redimensionner le reste
  3. effectuer PCA avec un certain nombre de composants
  4. mettre à l'échelle l'objet supprimé en fonction des paramètres obtenus en (2)
  5. prédire l'objet selon le modèle PCA
  6. calculer la PRESSE pour cet objet
  7. ré-exécuter le même algorithme sur d'autres objets
  8. résumer toutes les valeurs PRESS
  9. profit

Parce que je suis très nouveau dans le domaine, pour être sûr d'avoir raison, je compare les résultats avec la sortie de certains logiciels que j'ai (également pour écrire du code, je suis les instructions du logiciel). J'obtiens complètement les mêmes résultats en calculant la somme résiduelle des carrés et , mais le calcul de la PRESSE est un problème.R2

Pourriez-vous me dire si ce que j'implémente dans l'étape de validation croisée est correct ou non:

case 'loocv'

% # n - number of objects
% # p - number of variables
% # vComponents - the number of components used in CV
dataSets = divideData(n,n); 
         % # it is just a variable responsible for creating datasets for CV 
         % #  (for LOOCV datasets will be equal to [1, 2, 3, ... , n]);'
tempPRESS = zeros(n,vComponents);

for j = 1:n
  Xmodel1 = X; % # X - n x p original matrix
  Xmodel1(dataSets{j},:) = []; % # delete the object to be predicted
  [Xmodel1,Xmodel1shift,Xmodel1div] = skScale(Xmodel1, 'Center', vCenter, 
                                              'Scaling', vScaling); 
          % # scale the data and extract the shift and scaling factor
  Xmodel2 = X(dataSets{j},:); % # the object to be predicted
  Xmodel2 = bsxfun(@minus,Xmodel2,Xmodel1shift); % # shift and scale the object
  Xmodel2 = bsxfun(@rdivide,Xmodel2,Xmodel1div);
  [Xscores2,Xloadings2] = myNipals(Xmodel1,0.00000001,vComponents); 
          % # the way to calculate the scores and loadings
                % # Xscores2 - n x vComponents matrix
                % # Xloadings2 - vComponents x p matrix
  for i = 1:vComponents
    tempPRESS(j,i) = sum(sum((Xmodel2* ...
       (eye(p) - transpose(Xloadings2(1:i,:))*Xloadings2(1:i,:))).^2));
  end
end
PRESS = sum(tempPRESS,1);

Dans le logiciel ( PLS_Toolbox ) cela fonctionne comme ceci:

for i = 1:vComponents
    tempPCA = eye(p) - transpose(Xloadings2(1:i,:))*Xloadings2(1:i,:);
    for kk = 1:p
        tempRepmat(:,kk) = -(1/tempPCA(kk,kk))*tempPCA(:,kk);
          % # this I do not understand
        tempRepmat(kk,kk) = -1; 
          % # here is some normalization that I do not get
    end 
    tempPRESS(j,i) = sum(sum((Xmodel2*tempRepmat).^2)); 
end

Donc, ils font une normalisation supplémentaire en utilisant cette tempRepmatvariable: la seule raison pour laquelle j'ai trouvé était qu'ils appliquent LOOCV pour une PCA robuste. Malheureusement, l'équipe de support n'a pas voulu répondre à ma question car je n'ai qu'une version de démonstration de leur logiciel.

Kirill
la source
Vérification supplémentaire si je comprends l'extrait de normalisation supplémentaire: quel est le rôle de la tempRepmat(kk,kk) = -1ligne? La ligne précédente ne garantit-elle pas déjà que tempRepmat(kk,kk)-1 est égal à? Aussi, pourquoi les inconvénients? L'erreur va être quadrillée de toute façon, alors ai-je bien compris que si les inconvénients sont supprimés, rien ne changera?
Amoeba dit Reinstate Monica
Je le vérifiais dans le passé et rien ne changera. C'est correct. Je n'ai trouvé que quelques parallèles avec des implémentations PCA robustes car chaque valeur PRESS calculée dans une telle implémentation (avant de tout résumer) a son propre poids.
Kirill
Je cherche le code R équivalent au code MATLAB fourni dans la réponse et j'ai mis en place une prime.
AIM_BLB
3
@CSA, demander du code est hors sujet ici (même, probablement, via des commentaires et des primes). Vous devriez pouvoir demander cela sur Stack Overflow : vous pouvez copier le code, citer la source avec un lien ici, et demander une traduction. Je crois que tout cela serait là-dessus.
gung - Rétablir Monica

Réponses:

21

Ce que vous faites est faux: cela n'a aucun sens de calculer PRESS for PCA comme ça! Plus précisément, le problème réside dans votre étape # 5.


Approche naïve de PRESS pour PCA

Soit l'ensemble de données composé de points dans l'espace dimensionnel: . Pour calculer l'erreur de reconstruction pour un seul point de données de test , vous effectuez l'ACP sur l'ensemble d'apprentissage sans ce point, prenez un certain nombre d'axes principaux sous forme de colonnes de , et recherchez l'erreur de reconstruction comme . PRESS est alors égal à la somme de tous les échantillons de testd x ( i )R d ,ndx ( i ) X ( - i ) k U ( - i ) x ( i ) - x ( i ) 2 = x ( i ) - U ( - i ) [ U ( - i ) ] x ( i ) R E S Sx(i)Rd,i=1nx(i)X(i)kU(i) i Px(i)x^(i)2=x(i)U(i)[U(i)]x(i)2i, donc l'équation raisonnable semble être:

PRESS=?i=1nx(i)U(i)[U(i)]x(i)2.

Par souci de simplicité, j'ignore les problèmes de centrage et de mise à l'échelle ici.

L'approche naïve est fausse

Le problème ci-dessus est que nous utilisons pour calculer la prédiction , et c'est une très mauvaise chose.xx(i)x^(i)

Notez la différence cruciale dans un cas de régression, où la formule de l'erreur de reconstruction est fondamentalement la même , mais la prédiction est calculée en utilisant les variables prédictives et non en utilisant . Cela n'est pas possible en PCA, car en PCA il n'y a pas de variables dépendantes et indépendantes: toutes les variables sont traitées ensemble.y ( i )y(i)y^(i)2y^(i)y(i)

En pratique, cela signifie que la PRESSE telle que calculée ci-dessus peut diminuer avec un nombre croissant de composants et n'atteindre jamais un minimum. Ce qui laisserait penser que tous les composants sont significatifs. Ou peut-être que dans certains cas, il atteint un minimum, mais tend toujours à sur-adapter et à surestimer la dimensionnalité optimale.kd

Une approche correcte

Il existe plusieurs approches possibles, voir Bro et al. (2008) Validation croisée des modèles de composants: un regard critique sur les méthodes actuelles pour un aperçu et une comparaison. Une approche consiste à laisser de côté une dimension d'un point de données à la fois (c'est-à-dire au lieu de ), de sorte que les données d'apprentissage deviennent une matrice avec une valeur manquante , puis de prédire ("imputer") cette valeur manquante avec PCA. (On peut bien sûr tenir au hasard une fraction plus importante d'éléments de matrice, par exemple 10%). Le problème est que le calcul de l'ACP avec des valeurs manquantes peut être assez lent sur le plan informatique (il repose sur l'algorithme EM), mais doit être répété plusieurs fois ici. Mise à jour: voir http://alexhwilliams.info/itsneuronalblog/2018/02/26/crossval/ xxj(i)x(i) pour une discussion agréable et l'implémentation de Python (PCA avec des valeurs manquantes est implémenté via l'alternance des moindres carrés).

Une approche que j'ai trouvée beaucoup plus pratique consiste à omettre un point de données à la fois, à calculer l'ACP sur les données d'entraînement (exactement comme ci-dessus), puis à parcourir les dimensions de , les laisser un par un et calculer une erreur de reconstruction en utilisant le reste. Cela peut être assez déroutant au début et les formules ont tendance à devenir assez désordonnées, mais la mise en œuvre est plutôt simple. Permettez-moi d'abord de donner la formule (quelque peu effrayante), puis de l'expliquer brièvement:x(i)x(i)

PRESSPCA=i=1nj=1d|xj(i)[U(i)[Uj(i)]+xj(i)]j|2.

Considérez la boucle intérieure ici. Nous avons omis un point et calculé composantes principales sur les données d'apprentissage, . Maintenant, nous gardons chaque valeur comme test et utilisons les dimensions restantes pour effectuer la prédiction . La prédiction est la ème coordonnée de "la projection" (dans le sens des moindres carrés) de sur le sous-espace étendu par . Pour le calculer, trouvez un point dans l'espace PC plus proche dex(i)kU(i)xj(i)xj(i)Rd1x^j(i)jxj(i)U(i)z^Rkxj(i) en calculant où est avec -ième ligne expulsé et signifie pseudoinverse. Mappez maintenant vers l'espace d'origine: et prendre sa -ème coordonner . z^=[Uj(i)]+xj(i)RkUj(i) j[ ] + z U ( - i ) [U(i)j[]+z^U(i)[Uj(i)]+xj(i)j[]j

Une approximation de la bonne approche

Je ne comprends pas très bien la normalisation supplémentaire utilisée dans la PLS_Toolbox, mais voici une approche qui va dans le même sens.

Il existe une autre façon de mapper sur l'espace des composants principaux: , c'est-à-dire simplement prendre la transposition au lieu du pseudo-inverse. En d'autres termes, la dimension qui est laissée de côté pour les tests n'est pas comptée du tout et les poids correspondants sont également simplement supprimés. Je pense que cela devrait être moins précis, mais pourrait souvent être acceptable. La bonne chose est que la formule résultante peut maintenant être vectorisée comme suit (j'omet le calcul):jxj(i)z^approx=[Uj(i)]xj(i)

PRESSPCA,approx=i=1nj=1d|xj(i)[U(i)[Uj(i)]xj(i)]j|2=i=1n(IUU+diag{UU})x(i)2,

où j'ai écrit tant que pour la compacité, et signifie mettre tous les éléments non diagonaux à zéro. Notez que cette formule ressemble exactement à la première (PRESSE naïve) avec une petite correction! Notez également que cette correction ne dépend que de la diagonale de , comme dans le code PLS_Toolbox. Cependant, la formule est toujours différente de ce qui semble être implémenté dans PLS_Toolbox, et je ne peux pas expliquer cette différence. U d i a g {} U UU(i)Udiag{}UU

Mise à jour (février 2018): ci- dessus, j'ai appelé une procédure "correcte" et une autre "approximative" mais je ne suis plus si sûr que cela ait un sens. Les deux procédures ont du sens et je pense que ni l'une ni l'autre n'est plus correcte. J'aime vraiment que la procédure "approximative" ait une formule plus simple. De plus, je me souviens que j'avais un ensemble de données où la procédure "approximative" a donné des résultats qui semblaient plus significatifs. Malheureusement, je ne me souviens plus des détails.


Exemples

Voici comment ces méthodes se comparent pour deux jeux de données bien connus: le jeu de données Iris et le jeu de données Wine. Notez que la méthode naïve produit une courbe décroissante monotone, tandis que les deux autres méthodes produisent une courbe avec un minimum. Notez en outre que dans le cas d'Iris, la méthode approximative suggère 1 PC comme nombre optimal, mais la méthode pseudoinverse suggère 2 PC. (Et en regardant n'importe quel nuage de points PCA pour le jeu de données Iris, il semble que les deux premiers PC portent un signal.) Et dans le cas du vin, la méthode pseudoinverse pointe clairement vers 3 PC, tandis que la méthode approximative ne peut pas décider entre 3 et 5.

entrez la description de l'image ici


Code Matlab pour effectuer une validation croisée et tracer les résultats

function pca_loocv(X)

%// loop over data points 
for n=1:size(X,1)
    Xtrain = X([1:n-1 n+1:end],:);
    mu = mean(Xtrain);
    Xtrain = bsxfun(@minus, Xtrain, mu);
    [~,~,V] = svd(Xtrain, 'econ');
    Xtest = X(n,:);
    Xtest = bsxfun(@minus, Xtest, mu);

    %// loop over the number of PCs
    for j=1:min(size(V,2),25)
        P = V(:,1:j)*V(:,1:j)';        %//'
        err1 = Xtest * (eye(size(P)) - P);
        err2 = Xtest * (eye(size(P)) - P + diag(diag(P)));
        for k=1:size(Xtest,2)
            proj = Xtest(:,[1:k-1 k+1:end])*pinv(V([1:k-1 k+1:end],1:j))'*V(:,1:j)'; 
            err3(k) = Xtest(k) - proj(k);
        end

        error1(n,j) = sum(err1(:).^2);
        error2(n,j) = sum(err2(:).^2);
        error3(n,j) = sum(err3(:).^2);
    end    
end

error1 = sum(error1);
error2 = sum(error2);
error3 = sum(error3);
%// plotting code
figure
hold on
plot(error1, 'k.--')
plot(error2, 'r.-')
plot(error3, 'b.-')
legend({'Naive method', 'Approximate method', 'Pseudoinverse method'}, ...
    'Location', 'NorthWest')
legend boxoff
set(gca, 'XTick', 1:length(error1))
set(gca, 'YTick', [])
xlabel('Number of PCs')
ylabel('Cross-validation error')
amibe dit réintégrer Monica
la source
Merci pour votre réponse! Je connais ce papier. Et j'ai appliqué la validation croisée par ligne décrite ici (cela semble correspondre exactement au code que j'ai fourni). Je compare au logiciel PLS_Toolbox mais ils ont une ligne de code en LOOCV que je ne comprends vraiment pas (j'ai écrit dans la question d'origine).
Kirill
Oui, ils l'appellent "validation croisée par ligne" et votre implémentation semble correcte, mais veuillez noter que c'est une mauvaise façon d'effectuer la validation croisée (comme indiqué et également démontré empiriquement dans Bro et al.) Et, fondamentalement, vous devriez ne l'utilisez jamais! Concernant la ligne de code que vous ne comprenez pas - allez-vous mettre à jour votre question? Je ne sais pas de quoi vous parlez.
amibe dit Reinstate Monica
Le fait est que cette implémentation semble avoir la capacité d'atteindre le minimum en CV.
Kirill
La PRESSION de est très proche de la variance expliquée par LOO% - je ne dirais pas que c'est bon ou mauvais, mais c'est définitivement quelque chose dont il faut être conscient. Et comme la variance expliquée en% approchera 1 lorsque le modèle PCA se rapproche du rang de l'ensemble de données, ce X PRESS approchera 0.x^x
cbeleites mécontent de SX
@Kirill: Merci, l'extrait de code a du sens maintenant (vous pouvez peut-être supprimer les commentaires ci-dessus qui sont désormais obsolètes). Je n'ai aucune idée de ce qu'il est censé faire, mais si vous dites que l'erreur de prédiction calculée atteint un minimum, cela introduit probablement une pénalité pour les plus grands (nombre de composants; c'est-à-dire variable dans votre code). Honnêtement, je serais sceptique à l'égard d'une telle méthode (à moins qu'il n'y ait une justification théorique), d'autant plus qu'il existe de meilleures approches que celle que j'ai décrite dans ma réponse. ki
amibe dit Réintégrer Monica le
1

Pour ajouter un point encore plus général à la belle réponse de @ amoeba:

Une différence pratique et cruciale entre les modèles supervisés et non supervisés est que pour les modèles non supervisés, vous devez réfléchir beaucoup plus sérieusement à ce que vous considérez comme équivalent ou non.

Les modèles supervisés ont toujours leur sortie finale d'une manière où vous n'avez pas besoin de vous en préoccuper: par définition et par construction, prétend avoir la même signification que , vous pouvez donc le comparer directement. y yy^y^y

Afin de construire des mesures de performances significatives, vous devez réfléchir aux types de liberté du modèle qui n'ont pas de sens pour votre application et lesquels ne le sont pas. Cela conduirait à une PRESSE sur les scores, peut-être (généralement?) Après une sorte de rotation / retournement de type Procrustes.

APPUYEZ sur x Ma supposition est (je n'ai pas le temps maintenant de découvrir ce que font leurs 2 lignes de code - mais peut-être pourriez-vous parcourir les lignes et y jeter un œil?):

Afin d'obtenir une mesure utile pour déterminer une bonne complexité de modèle à partir d'une mesure qui donne une qualité d'ajustement qui augmentera généralement jusqu'à ce que le modèle de classement complet soit atteint, vous devez pénaliser les modèles trop complexes. Ce qui à son tour signifie que cette pénalisation est a) cruciale et b) l'ajustement de la pénalité ajustera la complexité choisie.


Note latérale: Je voudrais juste ajouter que je serais très prudent avec ce type d'optimisation automatisée de la complexité du modèle. D'après mon expérience, bon nombre de ces algorithmes ne produisent que de la pseudo-objectivité et ont souvent pour effet de ne bien fonctionner que pour certains types de données.

cbeleites mécontents de SX
la source