Wie kann eine Kreuzvalidierung für PCA durchgeführt werden, um die Anzahl der Hauptkomponenten zu bestimmen?

Ich versuche, meine eigene Funktion für die Hauptkomponentenanalyse, PCA, zu schreiben (natürlich ist bereits viel geschrieben, aber ich bin nur daran interessiert, Dinge selbst zu implementieren). Das Hauptproblem, auf das ich gestoßen bin, ist der Kreuzvalidierungsschritt und die Berechnung der vorhergesagten Quadratsumme (PRESS). Es spielt keine Rolle, welche Kreuzvalidierung ich verwende, es geht hauptsächlich um die dahinter stehende Theorie, aber erwägen Sie eine einmalige Kreuzvalidierung (LOOCV). Aus der Theorie habe ich herausgefunden, dass Sie zur Durchführung von LOOCV Folgendes tun müssen:

1. ein Objekt löschen
2. skalieren Sie den Rest
3. Führen Sie eine PCA mit einer bestimmten Anzahl von Komponenten durch
4. Skalieren Sie das gelöschte Objekt gemäß den in (2) erhaltenen Parametern.
5. Vorhersage des Objekts gemäß dem PCA-Modell
6. Berechnen Sie PRESS für dieses Objekt
7. Führen Sie denselben Algorithmus für andere Objekte erneut aus
8. Fassen Sie alle PRESS-Werte zusammen
9. profitieren

Da ich auf diesem Gebiet sehr neu bin, vergleiche ich die Ergebnisse mit der Ausgabe einer Software, die ich habe, um sicherzugehen, dass ich Recht habe (auch um Code zu schreiben, folge ich den Anweisungen in der Software). Ich erhalte die gleichen Ergebnisse bei der Berechnung der Restsumme von Quadraten und $R^2$ , aber die Berechnung von PRESS ist ein Problem.

Könnten Sie mir bitte sagen, ob das, was ich im Kreuzvalidierungsschritt implementiere, richtig ist oder nicht:

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);

In der Software ( PLS_Toolbox ) funktioniert das folgendermaßen:

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

Daher führen sie mit dieser tempRepmatVariablen eine zusätzliche Normalisierung durch : Der einzige Grund, den ich gefunden habe, war, dass sie LOOCV für robuste PCA anwenden. Leider wollte das Support-Team meine Frage nicht beantworten, da ich nur eine Demoversion seiner Software habe.

Was Sie tun, ist falsch: Es macht keinen Sinn, PRESS für PCA so zu berechnen! Insbesondere liegt das Problem in Schritt 5.

---

Naiver PRESS-Ansatz für PCA

Der Datensatz bestehe aus Punkten im d- dimensionalen Raum: x ( i ) ∈ R d$n$$d$ . Um den Rekonstruktionsfehler für einen einzelnen Testdatenpunkt x ( i ) zu berechnen, führen Sie eine PCA für den Trainingssatz X ( - i ) durch, wobei dieser Punkt ausgeschlossen ist. Nehmen Sie eine bestimmte Anzahl k von Hauptachsen als Spalten von U ( - i ) und findet die Rekonstruktionsfehler als ‖ x ( i ) - $\mathbf x^{(i)} \in \mathbb R^d, \, i=1 \dots n$$\mathbf x^{(i)}$$\mathbf X^{(-i)}$$k$$\mathbf U^{(-i)}$ - i ) [ U ( - i ) ] ≤ x. PRESS ist dann gleich Summe über alle Testprobeni, daher scheint die vernünftige Gleichung zu sein:PRESS ? = n ≤ i = 1 ≤ x ( i ) - U $\left \|\mathbf x^{(i)} - \hat{\mathbf x}^{(i)}\right\|^2 = \left\|\mathbf x^{(i)} - \mathbf U^{(-i)} [\mathbf U^{(-i)}]^\top \mathbf x^{(i)}\right\|^2$$i$

$$\mathrm{PRESS} \overset{?}{=} \sum_{i=1}^n \left\|\mathbf x^{(i)} - \mathbf U^{(-i)} [\mathbf U^{(-i)}]^\top \mathbf x^{(i)}\right\|^2.$$

Der Einfachheit halber ignoriere ich hier die Probleme der Zentrierung und Skalierung.

Der naive Ansatz ist falsch

Das Problem oben ist, dass wir $\mathbf x^{(i)}$ die Vorhersage zu berechnen x ( i ) , und das ist eine sehr schlechte Sache.$\hat{\mathbf x}^{(i)}$

Beachten Sie den entscheidenden Unterschied zu einem Regressions Fall, in dem die Formel für den Rekonstruktionsfehler ist im Grunde die gleiche , aber Prädiktion y ( i ) berechnet , die unter Verwendung von Vorhersagevariablen und nicht unter Verwendung von y ( i $\left\|\mathbf y^{(i)} - \hat{\mathbf y}^{(i)}\right\|^2$$\hat{\mathbf y}^{(i)}$$\mathbf y^{(i)}$ . Dies ist in PCA nicht möglich, da es in PCA keine abhängigen und unabhängigen Variablen gibt: Alle Variablen werden zusammen behandelt.

In der Praxis bedeutet dies, dass die oben berechnete PRESSE mit zunehmender Anzahl von Komponenten abnehmen und niemals ein Minimum erreichen kann. Was zu der Annahme führen würde, dass alle d- Komponenten von Bedeutung sind. Oder vielleicht erreicht es in einigen Fällen ein Minimum, neigt aber immer noch dazu, die optimale Dimensionalität zu überanpassen und zu überschätzen.$k$$d$

Ein korrekter Ansatz

Es gibt mehrere mögliche Ansätze, siehe Bro et al. (2008) Kreuzvalidierung von Komponentenmodellen: Ein kritischer Blick auf aktuelle Methoden für einen Überblick und Vergleich. Ein Ansatz besteht darin, jeweils eine Dimension eines Datenpunkts wegzulassen (dh anstelle von x ( i ) http://alexhwilliams.info/itsneuronalblog/2018/02/26/crossval/ für eine nette Diskussion und Python-Implementierung (PCA mit fehlenden Werten wird über alternierende kleinste Quadrate implementiert).$x^{(i)}_j$$\mathbf x^{(i)}$ ), so dass die Trainingsdaten zu einer Matrix mit einem fehlenden Wert werden, und dann vorherzusagen ("unterstellen") ") dieser fehlende Wert mit PCA. (Man kann natürlich zufällig einen größeren Teil der Matrixelemente heraushalten, z. B. 10%). Das Problem ist, dass das Berechnen von PCA mit fehlenden Werten sehr langsam sein kann (es basiert auf dem EM-Algorithmus), aber hier viele Male wiederholt werden muss. Update: siehe

Ein Ansatz, den ich als viel praktischer empfand, besteht darin, einen Datenpunkt $\mathbf x^{(i)}$$\mathbf x^{(i)}$

$$\mathrm{PRESS_{PCA}} = \sum_{i=1}^n \sum_{j=1}^d \left|x^{(i)}_j - \left[\mathbf U^{(-i)} \left [\mathbf U^{(-i)}_{-j}\right]^+\mathbf x^{(i)}_{-j}\right]_j\right|^2.$$

$\mathbf x^{(i)}$$k$$\mathbf U^{(-i)}$$x^{(i)}_j$$\mathbf x^{(i)}_{-j} \in \mathbb R^{d-1}$$\hat{x}^{(i)}_j$$j$$\mathbf x^{(i)}_{-j}$$\mathbf U^{(-i)}$$\hat{\mathbf z}$$\mathbb R^k$$\mathbf x^{(i)}_{-j}$$\hat{\mathbf z} = \left [\mathbf U^{(-i)}_{-j}\right]^+\mathbf x^{(i)}_{-j} \in \mathbb R^k$$\mathbf U^{(-i)}_{-j}$$\mathbf U^{(-i)}$$j$$[\cdot]^+$$\hat{\mathbf z}$$\mathbf U^{(-i)} \left [\mathbf U^{(-i)}_{-j}\right]^+\mathbf x^{(i)}_{-j}$$j$$[\cdot]_j$

Eine Annäherung an den richtigen Ansatz

Ich verstehe die zusätzliche Normalisierung, die in der PLS_Toolbox verwendet wird, nicht ganz, aber hier ist ein Ansatz, der in die gleiche Richtung geht.

$\mathbf x^{(i)}_{-j}$$\hat{\mathbf z}_\mathrm{approx} = \left [\mathbf U^{(-i)}_{-j}\right]^\top\mathbf x^{(i)}_{-j}$Nehmen Sie einfach die Transponierte anstelle der Pseudo-Inversen. Mit anderen Worten, die zum Testen ausgelassene Abmessung wird überhaupt nicht gezählt, und die entsprechenden Gewichte werden ebenfalls einfach herausgeschmissen. Ich denke, dies sollte weniger genau sein, könnte aber oft akzeptabel sein. Das Gute ist, dass die resultierende Formel nun wie folgt vektorisiert werden kann (ich lasse die Berechnung weg):

$$\begin{align} \mathrm{PRESS_{PCA,\,approx}} &= \sum_{i=1}^n \sum_{j=1}^d \left|x^{(i)}_j - \left[\mathbf U^{(-i)} \left [\mathbf U^{(-i)}_{-j}\right]^\top\mathbf x^{(i)}_{-j}\right]_j\right|^2 \\ &= \sum_{i=1}^n \left\|\left(\mathbf I - \mathbf U \mathbf U^\top + \mathrm{diag}\{\mathbf U \mathbf U^\top\}\right) \mathbf x^{(i)}\right\|^2, \end{align}$$

$\mathbf U^{(-i)}$$\mathbf U$$\mathrm{diag}\{\cdot\}$$\mathbf U \mathbf U^\top$

Update (Februar 2018): Oben habe ich eine Prozedur als "korrekt" und eine andere als "ungefähr" bezeichnet, bin mir aber nicht mehr so ​​sicher, ob dies sinnvoll ist. Beide Verfahren sind sinnvoll und ich denke, keines ist korrekter. Ich mag es wirklich, dass das "ungefähre" Verfahren eine einfachere Formel hat. Ich erinnere mich auch, dass ich einen Datensatz hatte, in dem das "ungefähre" Verfahren zu Ergebnissen führte, die aussagekräftiger aussahen. Leider erinnere ich mich nicht mehr an die Details.

---

Beispiele

Diese Methoden werden für zwei bekannte Datensätze verglichen: Iris-Datensatz und Wein-Datensatz. Beachten Sie, dass die naive Methode eine monoton abnehmende Kurve erzeugt, während die beiden anderen Methoden eine Kurve mit einem Minimum ergeben. Beachten Sie ferner, dass im Fall Iris die ungefähre Methode 1 PC als optimale Anzahl vorschlägt, die pseudoinverse Methode jedoch 2 PCs vorschlägt. (Wenn man sich ein PCA-Streudiagramm für den Iris-Datensatz ansieht, scheint es, dass beide ersten PCs ein Signal übertragen.) Im Weinfall zeigt die pseudoinverse Methode eindeutig auf 3 PCs, während die ungefähre Methode nicht zwischen 3 und 5 entscheiden kann.

---

Matlab-Code zur Durchführung einer Kreuzvalidierung und zur Darstellung der Ergebnisse

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')

Um @ noebas nette Antwort noch allgemeiner zu gestalten:

Ein praktischer und entscheidender Unterschied zwischen überwachten und unbeaufsichtigten Modellen besteht darin, dass Sie bei unbeaufsichtigten Modellen viel genauer überlegen müssen, was Sie als gleichwertig betrachten und was nicht.

Überwachte Modelle haben immer ihre endgültige Ausgabe $\hat y$ in einer Weise, in der Sie sich nicht viel darum kümmern müssen: per Definition und Konstruktion, $\hat y$ behauptet, die gleiche Bedeutung zu haben wie $y$, so können Sie es direkt vergleichen.

Um aussagekräftige Leistungskennzahlen zu erstellen, müssen Sie sich überlegen, welche Arten von Freiheit des Modells für Ihre Anwendung bedeutungslos sind und welche nicht. Das würde zu einer PRESSE der Partituren führen, möglicherweise (normalerweise?) Nach einer Art Procrustes-ähnlichen Rotation / Flip.

PRESSE auf x Meine Vermutung ist (ich habe jetzt keine Zeit herauszufinden, was ihre 2 Codezeilen bewirken - aber vielleicht könnten Sie durch die Zeilen gehen und einen Blick darauf werfen?):

Um eine Kennzahl zu erhalten, die nützlich ist, um eine gute Modellkomplexität aus einer Kennzahl zu ermitteln, die eine Anpassungsgüte ergibt, die sich normalerweise erhöht, bis das Modell mit vollem Rang erreicht ist, müssen Sie zu komplexe Modelle bestrafen. Was wiederum bedeutet, dass diese Bestrafung a) entscheidend ist und b) die Anpassung der Strafe die gewählte Komplexität anpasst.

---

Randnotiz: Ich möchte nur hinzufügen, dass ich bei dieser Art der automatisierten Optimierung der Modellkomplexität sehr vorsichtig sein würde. Nach meiner Erfahrung ergeben viele dieser Algorithmen nur Pseudoobjektivität und gehen oft zu Lasten einer guten Arbeitsweise nur für bestimmte Datentypen.