Jak przeprowadzić walidację krzyżową dla PCA w celu ustalenia liczby głównych składników?

Próbuję napisać własną funkcję do analizy głównych składników, PCA (oczywiście jest już dużo napisanych, ale jestem zainteresowany tylko implementacją różnych rzeczy). Głównym problemem, jaki napotkałem, jest krok weryfikacji krzyżowej i obliczanie przewidywanej sumy kwadratów (PRASA). Nie ma znaczenia, z której walidacji krzyżowej korzystam, chodzi głównie o teorię, ale zastanów się nad weryfikacją krzyżową z pominięciem jednej zasady (LOOCV). Z teorii dowiedziałem się, że aby wykonać LOOCV, musisz:

1. usuń obiekt
2. skaluj resztę
3. wykonać PCA z pewną liczbą składników
4. skaluj usunięty obiekt zgodnie z parametrami uzyskanymi w (2)
5. przewidzieć obiekt zgodnie z modelem PCA
6. obliczyć PRASA dla tego obiektu
7. ponownie wykonaj ten sam algorytm dla innych obiektów
8. zsumuj wszystkie wartości PRASY
9. zysk

Ponieważ jestem bardzo nowy w tej dziedzinie, aby mieć pewność, że mam rację, porównuję wyniki z danymi wyjściowymi z posiadanego oprogramowania (również w celu napisania kodu postępuję zgodnie z instrukcjami w oprogramowaniu). Mam zupełnie takie same wyniki obliczania pozostałego sumę kwadratów i $R^2$ , ale obliczenia PRESS jest problem.

Czy możesz mi powiedzieć, czy to, co wdrażam na etapie weryfikacji krzyżowej, jest prawidłowe, czy nie:

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

W oprogramowaniu ( PLS_Toolbox ) działa to tak:

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

Tak więc dokonują dodatkowej normalizacji za pomocą tej tempRepmatzmiennej: jedynym powodem, dla którego znalazłem, było zastosowanie LOOCV dla solidnego PCA. Niestety zespół wsparcia nie chciał odpowiedzieć na moje pytanie, ponieważ mam tylko wersję demonstracyjną ich oprogramowania.

To, co robisz, jest złe: nie ma sensu obliczać PRESS dla PCA w ten sposób! W szczególności problem leży w kroku 5.

---

Naiwne podejście do PRESS dla PCA

$n$$d$$\mathbf x^{(i)} \in \mathbb R^d, \, i=1 \dots n$$\mathbf x^{(i)}$$\mathbf X^{(-i)}$$k$$\mathbf U^{(-i)}$$\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.$$

Dla uproszczenia ignoruję tutaj kwestie centrowania i skalowania.

Naiwne podejście jest złe

Problem powyżej polega na tym, że używamy do obliczenia prognozy , i to jest bardzo zła rzecz.$\mathbf x^{(i)}$$\hat{\mathbf x}^{(i)}$

Zwróć uwagę na zasadniczą różnicę w przypadku regresji, w której wzór na błąd rekonstrukcji jest zasadniczo taki sam , ale predykcja jest obliczana przy użyciu zmiennych predykcyjnych, a nie przy użyciu . Nie jest to możliwe w PCA, ponieważ w PCA nie ma zmiennych zależnych i niezależnych: wszystkie zmienne są traktowane razem.$\left\|\mathbf y^{(i)} - \hat{\mathbf y}^{(i)}\right\|^2$$\hat{\mathbf y}^{(i)}$$\mathbf y^{(i)}$

W praktyce oznacza to, że obliczone powyżej PRESS może się zmniejszać wraz ze wzrostem liczby składników i nigdy nie osiągać minimum. Co doprowadziłoby do wniosku, że wszystkie składniki są znaczące. A może w niektórych przypadkach osiąga minimum, ale nadal ma tendencję do przeładowywania i przeceniania optymalnej wymiarowości.$k$$d$

Prawidłowe podejście

Istnieje kilka możliwych podejść, patrz Bro i in. (2008) Walidacja krzyżowa modeli komponentów: krytyczne spojrzenie na obecne metody przeglądu i porównania. Jednym podejściem jest pominięcie jednego wymiaru jednego punktu danych na raz (tj. zamiast ), aby dane treningowe stały się macierzą z jedną brakującą wartością , a następnie przewidzieć („przypisać”) tę brakującą wartość za pomocą PCA. (Oczywiście można losowo wypuścić nieco większą część elementów macierzy, np. 10%). Problem polega na tym, że obliczanie PCA z brakującymi wartościami może być obliczeniowo dość powolne (opiera się na algorytmie EM), ale musi być tutaj wielokrotnie powtarzane. Aktualizacja: patrz http://alexhwilliams.info/itsneuronalblog/2018/02/26/crossval/$x^{(i)}_j$$\mathbf x^{(i)}$ dla miłej dyskusji i implementacji w języku Python (PCA z brakującymi wartościami jest implementowany przez naprzemienne najmniejsze kwadraty).

Podejście, które uważam za bardziej praktyczne, polega na pomijaniu jednego punktu danych na raz, obliczaniu PCA na danych treningowych (dokładnie tak jak powyżej), ale następnie zapętlaniu wymiarów , pomiń je pojedynczo i oblicz resztę przy użyciu reszty. Na początku może to być mylące, a formuły stają się dość nieuporządkowane, ale implementacja jest raczej prosta. Pozwól mi najpierw podać (nieco przerażającą) formułę, a następnie krótko ją wyjaśnić:$\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.$$

Rozważ tutaj wewnętrzną pętlę. Pominęliśmy jeden punkt i obliczyliśmy głównych składników danych treningowych, . Teraz trzymamy każdą wartość jako test i używamy pozostałych wymiarów do wykonania prognozy . Prognozowanie jest współrzędną „rzutu” (w sensie najmniejszych kwadratów) na podprzestrzeń rozpiętą autor: . Aby to obliczyć, znajdź punkt w przestrzeni komputera który jest najbliżej$\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}$ obliczając gdzie to z -tym rzędem wyrzucony, a oznacza pseudoinwersję. Teraz mapuj powrotem do pierwotnego miejsca: i weź -tą współrzędną . $\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$

Przybliżenie prawidłowego podejścia

Nie do końca rozumiem dodatkową normalizację stosowaną w PLS_Toolbox, ale oto jedno podejście, które idzie w tym samym kierunku.

Istnieje inny sposób odwzorowania na przestrzeń głównych komponentów: , tzn. po prostu weź transpozycję zamiast pseudo-odwrotnej. Innymi słowy, wymiar, który jest pominięty do testowania, w ogóle nie jest liczony, a odpowiadające mu ciężary są po prostu wykreślane. Myślę, że powinno to być mniej dokładne, ale często może być do zaakceptowania. Dobrą rzeczą jest to, że wynikową formułę można teraz wektoryzować w następujący sposób (pomijam obliczenia):$\mathbf x^{(i)}_{-j}$$\hat{\mathbf z}_\mathrm{approx} = \left [\mathbf U^{(-i)}_{-j}\right]^\top\mathbf x^{(i)}_{-j}$

$$\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}$$

gdzie napisałem jako dla zwięzłości, a oznacza ustawienie wszystkich elementów niediagonalnych na zero. Zauważ, że ta formuła wygląda dokładnie tak samo jak pierwsza (naiwne PRASA) z małą korektą! Zauważ również, że ta korekta zależy tylko od przekątnej , jak w kodzie PLS_Toolbox. Jednak formuła wciąż różni się od tego, co wydaje się być zaimplementowane w PLS_Toolbox, a tej różnicy nie mogę wyjaśnić. U d i a g {⋅} U U$\mathbf U^{(-i)}$$\mathbf U$$\mathrm{diag}\{\cdot\}$$\mathbf U \mathbf U^\top$

Aktualizacja (luty 2018 r.): Powyżej nazwałem jedną procedurę „poprawną”, a drugą „przybliżoną”, ale nie jestem już tak pewien, czy to ma sens. Obie procedury mają sens i myślę, że żadna nie jest bardziej poprawna. Naprawdę podoba mi się, że procedura „przybliżona” ma prostszą formułę. Pamiętam też, że miałem zbiór danych, w którym procedura „przybliżona” dała wyniki, które wyglądały na bardziej znaczące. Niestety nie pamiętam już szczegółów.

---

Przykłady

Oto porównanie tych metod dla dwóch dobrze znanych zestawów danych: zestawu danych Iris i zestawu danych wina. Zauważ, że metoda naiwna tworzy krzywą monotonicznie malejącą, podczas gdy dwie pozostałe metody dają krzywą z minimum. Zauważ ponadto, że w przypadku Iris metoda przybliżona sugeruje 1 PC jako liczbę optymalną, ale metoda pseudoinwersyjna sugeruje 2 PC. (Patrząc na dowolny wykres rozrzutu PCA dla zestawu danych Iris, wydaje się, że oba pierwsze komputery PC niosą jakiś sygnał.) W przypadku wina metoda pseudoinwersyjna wyraźnie wskazuje na 3 komputery, podczas gdy metoda przybliżona nie może zdecydować między 3 a 5.

---

Kod Matlab do przeprowadzania weryfikacji krzyżowej i wykreślania wyników

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

Aby dodać jeszcze bardziej ogólny punkt do ładnej odpowiedzi @ amoeba:

Jedną praktyczną i kluczową różnicą między modelami nadzorowanymi i bez nadzoru jest to, że w przypadku modeli bez nadzoru należy znacznie bardziej przemyśleć, co uznamy za równoważne, a co nie.

Modele nadzorowane zawsze mają swoje ostateczne wyjście w sposób, w którym nie trzeba się przejmować o tym: z definicji i konstrukcji, twierdzi, że ma takie samo znaczenie jak , dzięki czemu można bezpośrednio porównać. y$\hat y$$\hat y$$y$

Aby skonstruować znaczące miary wydajności, musisz pomyśleć, jakie rodzaje swobody modelu są pozbawione znaczenia dla twojej aplikacji, a które nie. Doprowadziłoby to do NACIŚNIĘCIA wyników, być może (zwykle?) Po pewnego rodzaju rotacji / przerzucaniu przypominającym Procrustes.

NACISNĄĆ na x Domyślam się (nie mam teraz czasu, aby dowiedzieć się, co robią ich 2 wiersze kodu - ale może mógłbyś przejść przez te wiersze i spojrzeć?):

Aby uzyskać miarę, która jest przydatna do określenia dobrej złożoności modelu na podstawie miary, która daje dobroć dopasowania, która zwykle rośnie, aż do osiągnięcia modelu pełnej rangi, należy karać za zbyt złożone modele. Co z kolei oznacza, że ​​ta kara jest a) kluczowa i b) dostosowanie kary dostosuje wybraną złożoność.

---

Uwaga dodatkowa: Chciałbym dodać, że byłbym bardzo ostrożny z tego rodzaju automatyczną optymalizacją złożoności modelu. Z mojego doświadczenia wynika, że ​​wiele z tych algorytmów zapewnia jedynie pseudoobiektywność i często wiąże się to z kosztem dobrej pracy tylko w przypadku niektórych rodzajów danych.