如何對 PCA 進行交叉驗證以確定主成分的數量？

我正在嘗試編寫我自己的主成分分析函數，PCA（當然已經寫了很多，但我只是對自己實現一些東西感興趣）。我遇到的主要問題是交叉驗證步驟和計算預測平方和（PRESS）。我使用哪種交叉驗證並不重要，這主要是關於背後理論的問題，但請考慮留一法交叉驗證（LOOCV）。從理論上我發現，為了執行 LOOCV，您需要：

1. 刪除一個對象
2. 縮放其餘部分
3. 使用一些組件執行 PCA
4. 根據（2）中獲得的參數縮放刪除的對象
5. 根據 PCA 模型預測對象
6. 計算此對象的 PRESS
7. 對其他對象重新執行相同的算法
8. 總結所有 PRESS 值
9. 利潤

因為我是這個領域的新手，為了確保我是對的，我將結果與我擁有的一些軟件的輸出進行了比較（也為了編寫一些代碼，我按照軟件中的說明進行操作）。我得到完全相同的結果計算殘差平方和，但計算 PRESS 是個問題。

您能否告訴我我在交叉驗證步驟中實施的內容是否正確：

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

在軟件（PLS_Toolbox）中，它的工作方式如下：

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

因此，他們使用此tempRepmat變量進行了一些額外的標準化：我發現的唯一原因是他們將 LOOCV 應用於穩健的 PCA。不幸的是，支持團隊不想回答我的問題，因為我只有他們軟件的演示版本。

您所做的是錯誤的：像這樣為 PCA 計算 PRESS 是沒有意義的！具體來說，問題出在您的第 5 步。

---

PRESS for PCA 的樸素方法

讓數據集由點在維空間：. 計算單個測試數據點的重建誤差，你在訓練集上執行 PCA排除這一點，取一定數量主軸作為列，並找到重構誤差為. 然後 PRESS 等於所有測試樣本的總和，所以合理的等式似乎是：

為簡單起見，我在這裡忽略了居中和縮放的問題。

天真的方法是錯誤的

上面的問題是我們使用計算預測，這是一件非常糟糕的事情。

請注意與回歸案例的關鍵區別，其中重建誤差的公式基本相同, 但預測 使用預測變量而不是使用. 這在 PCA 中是不可能的，因為在 PCA 中沒有因變量和自變量：所有變量都被一起處理。

在實踐中，這意味著上面計算的 PRESS 會隨著組件數量的增加而減少並且永遠不會達到最低限度。這會導致人們認為所有組件很重要。或者在某些情況下它確實達到了最小值，但仍然傾向於過度擬合和高估最佳維度。

正確的做法

有幾種可能的方法，參見Bro 等人。(2008) 組件模型的交叉驗證：對當前方法的批判性研究，以進行概述和比較。一種方法是一次省略一個數據點的一維（即代替)，使訓練數據成為一個有一個缺失值的矩陣，然後用 PCA 預測（“估算”）這個缺失值。（當然可以隨機保留一些較大比例的矩陣元素，例如 10%）。問題是用缺失值計算 PCA 的計算速度可能非常慢（它依賴於 EM 算法），但這裡需要多次迭代。更新：請參閱http://alexhwilliams.info/itsneuronalblog/2018/02/26/crossval/以獲得很好的討論和 Python 實現（具有缺失值的 PCA 是通過交替最小二乘法實現的）。

我發現一種​​更實用的方法是省略一個數據點一次，在訓練數據上計算 PCA（完全如上），然後循環遍歷，一次只留下一個，然後使用其餘部分計算重建誤差。這在開始時可能會很混亂，並且公式往往會變得很混亂，但實現起來相當簡單。先給出（有點嚇人的）公式，然後簡單解釋一下：

考慮這裡的內循環。我們遺漏了一點併計算訓練數據的主要成分，. 現在我們保留每個值作為測試並使用剩餘的尺寸來執行預測。預測是個-“投影”的坐標（在最小二乘意義上）到跨越的子空間. 要計算它，找到一個點在 PC 空間最接近通過計算在哪裡是和-第 行被踢出，並且代表偽逆。現在地圖回到原來的空間：並採取它-th 坐標.

正確方法的近似值

我不太了解 PLS_Toolbox 中使用的額外規範化，但這是一種朝著相同方向發展的方法。

還有另一種映射方式在主成分空間上：，即簡單地採用轉置而不是偽逆。換句話說，測試時遺漏的維度根本不計算在內，相應的權重也被簡單地踢掉了。我認為這應該不太準確，但通常可以接受。好消息是生成的公式現在可以按如下方式向量化（我省略了計算）：

我寫的地方作為為了緊湊，和意味著將所有非對角元素設置為零。請注意，這個公式看起來與第一個公式（樸素的 PRESS）完全一樣，只是稍加修正！另請注意，此校正僅取決於對角線，就像在 PLS_Toolbox 代碼中一樣。但是，該公式與 PLS_Toolbox 中似乎實現的公式仍然不同，我無法解釋這種差異。

**更新（2018 年 2 月）：**上面我稱一個程序“正確”和另一個“近似”，但我不太確定這是否有意義。這兩個程序都有意義，我認為兩者都不正確。我真的很喜歡“近似”程序有一個更簡單的公式。另外，我記得我有一些數據集，其中“近似”過程產生的結果看起來更有意義。不幸的是，我已經不記得細節了。

---

例子

以下是這些方法對兩個著名數據集的比較：Iris 數據集和 wine 數據集。請注意，樸素方法產生單調遞減曲線，而其他兩種方法產生具有最小值的曲線。進一步注意，在 Iris 案例中，近似方法建議 1 個 PC 作為最佳數量，但偽逆方法建議 2 個 PC。（查看 Iris 數據集的任何 PCA 散點圖，似乎兩個第一台 PC 都帶有一些信號。）在葡萄酒的情況下，偽逆方法清楚地指向 3 台 PC，而近似方法無法確定 3 和 5 之間。

---

用於執行交叉驗證並繪製結果的 Matlab 代碼

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

引用自：https://stats.stackexchange.com/questions/93845