傳統的OLS(普通最小二乘)方法無法解決樣本資料的共線性(multicollinearity)問題，如果你的資料樣本中每個特徵變數具有共線性，那麼使用基於PCA的PCR和PLSR方法對資料樣本進行迴歸建立模型將會是一個不錯的選擇。PCA是一種資料降維方式，但同時保持了原始資料降維後的特性；PCR是在降維後的資料空間(英文裡常稱為score)上進行OLSR(普通最小二乘迴歸)，然後將回歸係數矩陣轉化為原始空間；PLSR則可以看成改進版的PCR，該方法通過X和Y資料集的交叉投影方法使得迴歸模型兼顧到了X和Y資料集的本質關聯，同時相比於PCR，在使用少數主成分的情況下具有更好的預測結果。
本文所有測試用資料集均來自Matlab，並使用Matlab封裝的迴歸方法，對自己實現的程式碼做了驗證，本文參考文獻及資料如下：
Reference:
[1] GELADI P, KOWALSKI B R. Partial least-squares regression: a tutorial [J]. Analytica chimica acta, 1986, 185(1-17).
[2] WU F Y, ASADA H H. Implicit and intuitive grasp posture control for wearable robotic fingers: A data-driven method using partial least squares [J]. IEEE Transactions on Robotics, 2016, 32(1): 176-86.
[3] https://en.wikipedia.org/wiki/Ordinary_least_squares
[4] https://en.wikipedia.org/wiki/Principal_component_analysis
[5] https://en.wikipedia.org/wiki/Principal_component_regression
[6] https://en.wikipedia.org/wiki/Partial_least_squares_regression
完整Matlab程式碼實現: https://github.com/ShieldQiQi/PCA-PCR-PLSR-Matlab-code
一、OLSR
即為普通最小二乘迴歸，對此我們應該十分熟悉，各種大物材料力學實驗都會用到這種方法，只不過我們當時使用的單變數的資料，當資料集涉及到矩陣，多維變數的形式時，就需要使用更加普遍適用的模型，我們設原始資料自變數(independent value)矩陣為$ X∈R_{n{\times}m} $，即X資料集含有n個樣本，每個樣本有m個特徵變數；設原始資料因變數(dependent value)矩陣為$ Y∈R_{n{\times}p} $，即Y資料集含有n個樣本，每個樣本有p個特徵變數。構建的最小二乘迴歸模型為：
$$ Y=X{\cdot}B+E \tag{1} $$
上式中$ B∈R_{m{\times}p} $為迴歸模型的係數矩陣，$ E∈R_{n{\times}p} $為模型預測的殘差。B的通用解法參考維基百科，為：
$$ B=(X^{T}X)^{-1}X^{T}Y \tag{2} $$
二、PCA
PCA本質上是一種建立一種維度小於原始資料維度(特徵變數數)正交基底空間，將原始資料投影到新的低維空間，以達到資料降維而保持原有特性的方法。PCA的步驟為：
1.對原始資料進行列居中處理: X(:,j) = X(:,j) - mean(X(:,j))
2.計算協方差矩陣$ X^{T}X $的前num大個特徵值和對應的特徵向量(此處num即為我們需要使用的主成分個數)
3.取前num個特徵向量(作為列向量)組成係數矩陣P
4.通過公式 $ T=XP $ 即可求得在新空間下的降維後的(原來維度為m，降維後為num)資料矩陣T，英文裡稱為score，P稱為loading
至於為什麼這樣做，PCA的原理可以參考維基百科，或者我的這篇博文: https://www.cnblogs.com/QiQi-Robotics/p/14303718.html
在實際應用中，計算協方差矩陣的特徵向量常採用迭代計算的方式，常用的方法為NIPALS，Matlab精簡程式碼(Matlab使用的為散佈矩陣，而我的程式碼為協方差矩陣，所以特徵值會相差(n-1)倍)實現如下：
1 % 迭代得到num個成份 2 for h = 1:num 3 % step(1) 4 % --------------------------------------------------------------------- 5 % 取T(:,h)為任意一個X_centered中的列向量，此處直接取第一列 6 T(:,h) = X_iteration(:,1); 7 8 % step(2) to step(5) 9 % 迭代直到收斂到容忍度內的主成分 10 while(1) 11 P(:,h) = X_iteration'*T(:,h)/(T(:,h)'*T(:,h)); 12 % 歸一化P(:,h) 13 P(:,h) = P(:,h)/sqrt(P(:,h)'*P(:,h)); 14 t_temp = T(:,h); 15 T(:,h) = X_iteration*P(:,h)/(P(:,h)'*P(:,h)); 16 17 % 檢查當前T(:,h)與上一步T(:,h)是否相等以決定是否繼續迭代 18 if max(abs(T(:,h)-t_temp)) <= tolerance 19 % 儲存按順序排列的特徵值 20 % 注意此處的特徵值為協方差矩陣的特徵值，而matlab PCA方法使用的為散佈矩陣(離散度矩陣)，故後者的特徵值為前者的(n-1)倍 21 eigenValues(h) = P(:,h)'*(X_centered'*X_centered)*P(:,h); 22 break; 23 else 24 end 25 end 26 27 % 計算殘差，更新資料矩陣 28 % --------------------------------------------------------------------- 29 X_iteration = X_iteration - T(:,h)*P(:,h)'; 30 end
三、PCR
PCR使用的迴歸方法是OLSR，只不過迴歸的模型是建立在主成分空間，以防止原始資料的共線性問題導致模型建立不準確，步驟如下：
1.執行PCA對原始資料進行降維處理
2.對新資料矩陣T(score)(選多少列，就是利用多少個主成分)和居中(mean-centered)後的Y建立OLSR迴歸模型，得到主成分空間中的迴歸係數矩陣$ B^{'} $
3.最終原始空間的係數矩陣$ B=P{\cdot}B^{'} $，該步可以將 $ T=XP $ 代入到式(1)中推導而得(利用$ PP^{T}=E $)
4.當我們需要回歸新的到的資料X*時，將該資料對減去原始模型資料X的均值，代入到迴歸模型，得到預測的$Y^{'}$，然後該矩陣加上原始模型資料Y的均值即為最終的結果
Matlab精簡程式碼如下：
1 % 定義測試集樣本的數量 2 r = n; 3 % 將原始資料降維到主成分空間(T)後，使用OLS最小二乘迴歸獲取係數矩陣 4 B_inPca = inv(T'*T)*T'*Y_centered; 5 %B_inPca = regress(Y-mean(Y), T(:,1:num)); 6 % 將係數矩陣從主成分空間轉化到原始空間 7 B_estimated = P*B_inPca; 8 9 % 定義測試集，此處直接使用原始資料的前r行 10 X_validate = zeros(r,m); 11 % 對原始資料集居中列平均化 12 for j = 1:m 13 % 注意，此處減去的平均值應該為模型資料集的平均值，而非新資料的平均值 14 X_validate(:,j) = X(1:r,j) - mean(X(:,j)); 15 end 16 17 Y_estimated = X_validate*B_estimated; 18 for i = 1:p 19 % 注意此處最終的輸出需要加上資料集Y的均值 20 Y_estimated(:,i) = Y_estimated(:,i) + mean(Y(:,i)); 21 end
四、PLSR
PLSR相對於PCR的一個優點在於在使用更少的主成分可以獲得更具有魯棒性的預測結果(具體可以檢視Matlab中關於PLSR的幫助文件)，具體步驟查閱論文 [1]。精簡版Matlab程式碼如下：
1.建立模型部分
1 % 迭代得到num個成份 2 for h = 1:num 3 % step(1) 4 % --------------------------------------------------------------------- 5 % 取u_h為任意一個Y_centered中的列向量，此處直接取第一列 6 U(:,h) = Y_centered(:,1); 7 8 % step(2) to step(8) 9 % --------------------------------------------------------------------- 10 while 1 11 % 在資料矩陣X_centered中 12 W(:,h) = X_centered'*U(:,h)/(U(:,h)'*U(:,h)); 13 % 對資料進行歸一化 14 W(:,h) = W(:,h)/sqrt(W(:,h)'*W(:,h)); 15 t_temp = T(:,h); 16 T(:,h) = X_centered*W(:,h)/(W(:,h)'*W(:,h)); 17 18 % 在資料矩陣Y_centered中 19 Q(:,h) = Y_centered'*T(:,h)/(T(:,h)'*T(:,h)); 20 % 對資料進行歸一化 21 Q(:,h) = Q(:,h)/sqrt(Q(:,h)'*Q(:,h)); 22 U(:,h) = Y_centered*Q(:,h)/(Q(:,h)'*Q(:,h)); 23 24 % 檢查T(:,h)與T(:,h)的前一步是否相等，若小於某個數值則該PLS成份迭代完成，否則返回繼續迭代 25 if max(abs(T(:,h)-t_temp)) <= tolerance 26 break; 27 else 28 end 29 end 30 31 % step(9) to step(13) 32 % --------------------------------------------------------------------- 33 P(:,h) = X_centered'*T(:,h)/(T(:,h)'*T(:,h)); 34 % 對資料進行歸一化 35 p_norm = sqrt(P(:,h)'*P(:,h)); 36 P(:,h) = P(:,h)/p_norm; 37 T(:,h) = T(:,h)*p_norm; 38 W(:,h) = W(:,h)*p_norm; 39 B(h) = U(:,h)'*T(:,h)/(T(:,h)'*T(:,h)); 40 41 % 計算殘差，更新資料矩陣 42 % --------------------------------------------------------------------- 43 X_centered = X_centered - T(:,h)*P(:,h)'; 44 Y_centered = Y_centered - B(h)*T(:,h)*Q(:,h)'; 45 end
2.預測部分
1 % 對原始資料集居中列平均化 2 for j = 1:m 3 % 注意，此處減去的平均值應該為模型資料集的平均值，而非新資料的平均值 4 X_validate(1:r,j) = X(1:r,j) - mean(X(:,j)); 5 end 6 7 % 計算預測的T 8 for h = 1:num 9 T_est(:,h) = X_validate*W(:,h); 10 X_validate = X_validate - T_est(:,h)*P(:,h)'; 11 end 12 13 % 計算預測的Y 14 for h = 1:num 15 Y_estimated = Y_estimated + B(h)*T_est(:,h)*Q(:,h)'; 16 end 17 for i = 1:p 18 % 注意此處最終的輸出需要加上資料集Y的均值 19 Y_estimated(:,i) = Y_estimated(:,i) + mean(Y(:,i)); 20 end
五、實驗結果
圖1 Matlab PLSR演算法(SIMPLS)和自定義PLSR(NIPALS)方法效果對比