簡介KPLS核偏最小二乘回歸的MATLAB實(shí)現(xiàn)程序包面向需處理高維、非線性數(shù)據(jù)關(guān)聯(lián)的科研人員與工程師解決難以用線性PLS建模的問題。壓縮包共6個(gè)文件以MATLAB的m源碼和mat格式數(shù)據(jù)文件為主另含asv自動(dòng)備份文件整體僅17KB輕量易用。其中m文件實(shí)現(xiàn)核矩陣計(jì)算、模型訓(xùn)練與預(yù)測的核心流程mat文件內(nèi)置完整示例數(shù)據(jù)集用戶可直接運(yùn)行驗(yàn)證算法效果并借助代碼理解高斯核、多項(xiàng)式核等的選擇邏輯——無需預(yù)先設(shè)定數(shù)據(jù)的非線性結(jié)構(gòu)。已有1044人學(xué)習(xí)下載說明該程序包在非線性回歸與模式識別場景中具備一定參考價(jià)值。借助該資源讀者既能拆解KPLS的算法實(shí)現(xiàn)細(xì)節(jié)又能將示例數(shù)據(jù)替換為自己的數(shù)據(jù)開展建模分析省去從零編碼的繁瑣過程。1. 拿到KPLS程序后先想清楚它要解決什么核偏最小二乘問題KPLS核偏最小二乘程序在Matlab里一直是個(gè)有點(diǎn)“玄學(xué)”的存在代碼看著不長跑起來卻容易翻車。你在近紅外光譜定量分析、化工過程軟測量或者多元回歸任務(wù)里一旦遇到“線性PLS的預(yù)測殘差壓不下去”的情況通常就該把KPLS拿出來試了——它用核函數(shù)把原始數(shù)據(jù)映射到高維特征空間再在這個(gè)空間里做偏最小二乘回歸本質(zhì)上是給PLS加了一層“非線性外掛”。這篇筆記我就按自己實(shí)際調(diào)這套程序的路徑把原理、Matlab實(shí)現(xiàn)、參數(shù)整定和踩過的坑完整講一遍。新手跟著代碼能跑通熟手可以重點(diǎn)看第5章的邊界條件。2. 從線性PLS到KPLS核矩陣中心化是第一個(gè)分水嶺2.1 線性PLS為什么在非線性數(shù)據(jù)上失靈偏最小二乘回歸PLS的核心操作是在輸入矩陣X和輸出矩陣Y之間迭代提取潛變量每一步找一組權(quán)重w和c讓X的得分t和Y的得分u之間的協(xié)方差最大化然后做回歸、降維重復(fù)直到提取完設(shè)定數(shù)量的潛變量。這個(gè)過程的每一步都是線性投影所以當(dāng)X和Y之間存在明顯的非線性關(guān)系時(shí)——比如光譜吸收隨濃度變化呈非線性、化學(xué)反應(yīng)速率隨溫度呈指數(shù)變化——線性PLS的潛變量很難用一個(gè)平面把數(shù)據(jù)云壓平殘差里會殘留系統(tǒng)性結(jié)構(gòu)。常見的補(bǔ)救辦法是先把X做非線性擴(kuò)維加入x2、xy、log(x)這類手工特征但特征工程需要靠經(jīng)驗(yàn)試而且維度爆炸快。KPLS的思路更直接不顯式構(gòu)造高維特征而是用核函數(shù)計(jì)算樣本兩兩之間的內(nèi)積相當(dāng)于在隱含的高維空間里做PLS。這個(gè)空間里的線性關(guān)系對應(yīng)原始空間的非線性關(guān)系。2.2 核矩陣K和M矩陣KPLS的兩根支柱在Matlab實(shí)現(xiàn)KPLS之前要先理解兩個(gè)矩陣。第一個(gè)是核矩陣K維度是n×nn為樣本數(shù)K(i,j) k(x_i, x_j)k是核函數(shù)。最常用的RBF核形式為k(x_i, x_j) exp(-||x_i - x_j||2 / (2σ2))其中σ是核寬度也叫尺度參數(shù)。第二個(gè)是中心化后的核矩陣Kc它必須經(jīng)過以下變換才能用于回歸Kc K - 1_n·K - K·1_n 1_n·K·1_n其中1_n是n×n、所有元素為1/n的矩陣。這一步做的是“在特征空間里把數(shù)據(jù)云的中心移到原點(diǎn)”。很多初版KPLS代碼跑出來的結(jié)果比線性PLS還差往往就是漏了這一項(xiàng)。線性PLS對X做的是均值中心化而核空間里沒法直接對高維特征做均值只能通過核矩陣的代數(shù)操作完成漏掉等于拿一個(gè)帶偏移量的K去做迭代權(quán)重w的每個(gè)內(nèi)積項(xiàng)都帶常數(shù)偏差預(yù)測自然漂移。2.3 KPLS與線性PLS的迭代差異線性PLS每次迭代更新的是X和Y本身做完一個(gè)潛變量后用X減去得分t對載荷p的投影Y減去得分t對回歸系數(shù)的投影。KPLS迭代時(shí)更新的是核矩陣而不是原始數(shù)據(jù)——因?yàn)樵紨?shù)據(jù)在高維空間里不可見所有運(yùn)算都發(fā)生在核矩陣上。具體迭代時(shí)對Kc做潛變量提取每一步都要用殘差矩陣替換當(dāng)前Kc殘差的計(jì)算和線性PLS類似只是矩陣維度始終是n×n與原始X的變量數(shù)無關(guān)。這一點(diǎn)讓KPLS對“變量數(shù)遠(yuǎn)大于樣本數(shù)”的光譜數(shù)據(jù)特別友好也是它能處理高維輸入的核心原因。3. 在Matlab里把KPLS程序?qū)懲ê诵牡c參數(shù)定義3.1 訓(xùn)練階段RBF核矩陣、中心化與潛變量提取的完整代碼下面這段代碼是我自己整理后的KPLS訓(xùn)練主流程輸入Xtr是訓(xùn)練集特征矩陣每行一個(gè)樣本Ytr是對應(yīng)的輸出矩陣單輸出時(shí)是n×1向量。核函數(shù)先用RBF參數(shù)sigma控制核寬度。需要Matlab 2016a以上版本即可不需要額外工具箱。function [Kc, T, W, P, B, sigma] kpls_train(Xtr, Ytr, A, sigma) % KPLS_train: 核偏最小二乘訓(xùn)練 % 輸入: % Xtr - 訓(xùn)練輸入矩陣n×d每行一個(gè)樣本 % Ytr - 訓(xùn)練輸出矩陣n×m單輸出時(shí) m1 % A - 潛變量個(gè)數(shù)正整數(shù) % sigma - RBF核寬度標(biāo)量正數(shù) % 輸出: % Kc - 中心化核矩陣 n×n % T - 潛變量得分矩陣 n×A % W - 權(quán)重系數(shù)矩陣 n×A % P - 載荷矩陣 n×A % B - 回歸系數(shù)矩陣 n×m用于新樣本預(yù)測 % 注意: 調(diào)用前需對Xtr和Ytr做標(biāo)準(zhǔn)化否則sigma含義會失真 n size(Xtr, 1); % 1. 計(jì)算原始核矩陣RBF核 K zeros(n, n); for i 1:n for j 1:n diff Xtr(i,:) - Xtr(j,:); K(i,j) exp(-(diff * diff) / (2 * sigma^2)); end end % 2. 核矩陣中心化關(guān)鍵步驟漏掉這一項(xiàng)模型基本不能用 OneN ones(n, n) / n; Kc K - OneN * K - K * OneN OneN * K * OneN; % 3. 潛變量提取迭代 T zeros(n, A); W zeros(n, A); P zeros(n, A); B zeros(n, A); U zeros(n, A); Kres Kc; % 殘差核矩陣每迭代一次更新一次 for a 1:A % 3.1 對當(dāng)前殘差核矩陣做特征值分解取最大特征值對應(yīng)特征向量作為權(quán)重 % 注意: Kres 可能非對稱需轉(zhuǎn)成對稱矩陣處理 Ksym (Kres Kres) / 2; [V, D] eig(Ksym); [~, idx] max(diag(D)); t V(:, idx); % 得分向量 % 3.2 計(jì)算該潛變量的載荷和回歸系數(shù) p Kres * t / (t * t); u Ytr * t / (t * t); w p / norm(p); % 3.3 保存結(jié)果 T(:, a) t; W(:, a) w; P(:, a) p; U(:, a) u; B(:, a) w; % 3.4 更新殘差核矩陣減去當(dāng)前潛變量解釋的部分 Kres Kres - (p * p) * (t * t) / (t * t); end % 4. 計(jì)算回歸系數(shù)矩陣n×m B T * (T * Kc * T) \ T * Ytr; end這段代碼的核心思路是先算RBF核矩陣并做中心化然后在中心化核矩陣上逐次提取潛變量。每次迭代用eig分解找最大特征向量——注意這里對Kres做了對稱化處理因?yàn)閿?shù)值誤差會讓核矩陣出現(xiàn)極小非對稱量直接eig可能拿到復(fù)數(shù)特征向量。潛變量數(shù)A需要預(yù)先指定可以用交叉驗(yàn)證來選我一般從2開始試到15步長1。sigma的取值則和數(shù)據(jù)尺度強(qiáng)相關(guān)標(biāo)準(zhǔn)化后的數(shù)據(jù)通常從1.0開始按對數(shù)網(wǎng)格向兩端搜索。3.2 預(yù)測階段新樣本的核行與中心化不能直接套訓(xùn)練代碼預(yù)測新樣本時(shí)要構(gòu)造新樣本與所有訓(xùn)練樣本的核行向量然后做“行向中心化”——這是最容易寫錯(cuò)的地方。很多人在預(yù)測階段直接拿原始核行去和訓(xùn)練得到的B相乘結(jié)果預(yù)測曲線嚴(yán)重偏離原因就是沒有對核行做和訓(xùn)練階段一致的中心化操作。function Ypred kpls_predict(Xtr, Xte, Kc, B, sigma) % kpls_predict: KPLS預(yù)測 % 輸入: % Xtr - 訓(xùn)練輸入矩陣 n×d % Xte - 測試輸入矩陣 m×d每行一個(gè)新樣本 % Kc - 訓(xùn)練階段輸出的中心化核矩陣 n×n % B - 訓(xùn)練階段輸出的回歸系數(shù) n×m % sigma - 與訓(xùn)練階段相同的核寬度 % 輸出: % Ypred - 預(yù)測輸出 m×m n size(Xtr, 1); mte size(Xte, 1); % 1. 構(gòu)造測試核矩陣每個(gè)測試樣本與所有訓(xùn)練樣本的RBF核 Kte zeros(mte, n); for i 1:mte for j 1:n diff Xte(i,:) - Xtr(j,:); Kte(i,j) exp(-(diff * diff) / (2 * sigma^2)); end end % 2. 測試核矩陣行向中心化必須與訓(xùn)練階段的中心化公式呼應(yīng) OneN ones(n, n) / n; OneM ones(mte, n) / n; Kte_c Kte - OneM * K - OneM * K * OneN OneM * Kc; % 注意這里的K需要原始核矩陣 % 3. 預(yù)測 Ypred Kte_c * B; end但這里有一個(gè)隱患kpls_predict函數(shù)需要用到原始未中心化的訓(xùn)練核矩陣K而kpls_train的輸出里沒有包含K。我實(shí)際使用時(shí)會改一下訓(xùn)練函數(shù)把K也作為輸出返回或者干脆在預(yù)測函數(shù)里重新算一遍Kte和K。兩種做法都行重點(diǎn)是保持中心化公式一致。行向中心化的數(shù)學(xué)含義是把新樣本的核行投影到訓(xùn)練樣本張成的特征空間里然后減去訓(xùn)練核矩陣的列均值、加上全局均值這樣才能保證新樣本和訓(xùn)練樣本在同一坐標(biāo)系下比較。3.3 數(shù)據(jù)標(biāo)準(zhǔn)化不是可選項(xiàng)是sigma生效的前提我在3.1代碼注釋里寫了“調(diào)用前需對Xtr和Ytr做標(biāo)準(zhǔn)化”這里展開說為什么。RBF核的距離項(xiàng)||x_i-x_j||2對特征的尺度極其敏感如果特征A取值范圍0~1特征B取值范圍0~10000那么距離幾乎完全由特征B主導(dǎo)核矩陣會退化成只有單個(gè)特征起作用的“退化核”。最直接的后果是sigma不管怎么調(diào)模型都只盯著量級最大的特征KPLS的非線性擬合能力發(fā)揮不出來。常見做法是對每個(gè)特征做z-score標(biāo)準(zhǔn)化即減去均值除以標(biāo)準(zhǔn)差讓所有特征同尺度。Ytr的標(biāo)準(zhǔn)化同樣重要否則預(yù)測結(jié)果會被Ytr中量級大的輸出變量主導(dǎo)。標(biāo)準(zhǔn)化參數(shù)均值和標(biāo)準(zhǔn)差必須保存下來預(yù)測新樣本時(shí)用同一套參數(shù)轉(zhuǎn)換不能重新計(jì)算測試集的均值。4. KPLS程序的三組關(guān)鍵參數(shù)核寬度、潛變量數(shù)與核函數(shù)選型4.1 核寬度sigma的含義與搜索范圍sigma在RBF核中控制的是“局部性”sigma越小核值隨距離衰減越快模型越傾向于記住訓(xùn)練樣本附近的細(xì)節(jié)sigma越大核值衰減越平緩模型越接近線性核。極端情況下sigma趨近無窮大RBF核退化為常數(shù)核KPLS退化成類似線性PLS的效果。反過來sigma趨近0核矩陣趨近單位陣模型變成對訓(xùn)練樣本的逐點(diǎn)記憶泛化能力清零。我一般在標(biāo)準(zhǔn)化數(shù)據(jù)上從sigma0.1到sigma10按對數(shù)等間隔取10個(gè)值配合交叉驗(yàn)證一起搜。光譜數(shù)據(jù)我自己的經(jīng)驗(yàn)是sigma在0.5~2之間表現(xiàn)穩(wěn)定超出這個(gè)區(qū)間要么過擬合要么欠擬合。4.2 潛變量數(shù)A的選擇不是越多越好潛變量數(shù)A對應(yīng)的是“模型用幾個(gè)隱成分來解釋X-Y的關(guān)系”。A太小模型欠擬合殘差里還有結(jié)構(gòu)A太大模型開始擬合噪聲測試集誤差反彈。判斷A的合適范圍有個(gè)簡單手段畫出訓(xùn)練集和驗(yàn)證集的R2隨A的變化曲線兩條線一起上升時(shí)說明還能加潛變量訓(xùn)練集繼續(xù)上升而驗(yàn)證集開始下降的拐點(diǎn)就是合適的A。光譜數(shù)據(jù)里A5~10比較常見但這不絕對。另外注意KPLS的潛變量數(shù)并不等價(jià)于線性PLS的潛變量數(shù)因?yàn)楹丝臻g維度高得多前幾個(gè)潛變量解釋的方差比例通常也比線性PLS高得多。4.3 核函數(shù)選型對比RBF、線性核與多項(xiàng)式核核函數(shù)表達(dá)式適用場景需要注意的問題線性核x_i·x_j數(shù)據(jù)本身接近線性與線性PLS結(jié)果幾乎一致不需要用KPLSRBF核高斯核exp(-‖x_i-x_j‖2/2σ2)大多數(shù)非線性回歸問題光譜、軟測量需要調(diào)sigma對數(shù)據(jù)標(biāo)準(zhǔn)化敏感多項(xiàng)式核(x_i·x_j c)^d數(shù)據(jù)有明確的多項(xiàng)式關(guān)系階數(shù)d不宜超過3否則數(shù)值溢出Sigmoid核tanh(a·x_i·x_j b)偶爾用于某些特定分布參數(shù)a、b難調(diào)不是每個(gè)數(shù)據(jù)集都正定我平時(shí)90%的情況直接用RBF核原因是它只有一個(gè)參數(shù)sigma要調(diào)而且理論上可以逼近任意連續(xù)函數(shù)搜索結(jié)果一般不會太差。多項(xiàng)式核對階數(shù)敏感d2或3還行再高核矩陣就容易出現(xiàn)數(shù)值爆炸。線性核則沒必要用KPLS包裝直接跑線性PLS更省事。5. KPLS程序常見的坑從奇異核矩陣到預(yù)測漂移的排查5.1 漏掉核矩陣中心化所有預(yù)測值整體偏移現(xiàn)象訓(xùn)練集擬合很好但測試集預(yù)測值整體偏高或偏低一個(gè)常數(shù)殘差圖里出現(xiàn)明顯的水平偏移帶。原因訓(xùn)練階段沒有對核矩陣做中心化回歸系數(shù)B里包含了數(shù)據(jù)均值偏移量預(yù)測自然跟著偏。解決檢查訓(xùn)練代碼中是否存在Kc K - OneNK - KOneN OneNKOneN這一行不存在就是漏了。加上后重新訓(xùn)練預(yù)測就會回到正確水平。5.2 sigma調(diào)得過小訓(xùn)練集R2接近1但驗(yàn)證集慘不忍睹現(xiàn)象sigma從1.0改成0.1后訓(xùn)練集R2從0.85跳到0.98驗(yàn)證集R2反而從0.80跌到0.55。原因sigma太小核矩陣對角線遠(yuǎn)大于非對角線KPLS實(shí)際上在記憶訓(xùn)練樣本核空間里的局部鄰域小到只剩樣本自身。解決把sigma往大方向調(diào)觀察訓(xùn)練集與驗(yàn)證集R2的差距是否縮小。如果差距仍然很大考慮是不是數(shù)據(jù)本身噪聲太高需要先做平滑或變量篩選。5.3 預(yù)測階段核行中心化公式寫錯(cuò)在線預(yù)測結(jié)果抖動(dòng)現(xiàn)象離線測試集效果正常但換到在線預(yù)測場景每次只來一個(gè)新樣本結(jié)果偶爾跳變且偏差方向不固定。原因預(yù)測時(shí)核行中心化用了錯(cuò)誤的列均值。新樣本的核行向量長度等于訓(xùn)練樣本數(shù)n中心化需要同時(shí)減去訓(xùn)練核矩陣的列均值向量每個(gè)元素是所有訓(xùn)練樣本與該列樣本的核均值和行向量均值再補(bǔ)回全局均值。漏掉其中任何一項(xiàng)都會導(dǎo)致預(yù)測對訓(xùn)練樣本的索引順序敏感。解決把訓(xùn)練時(shí)的K保存下來預(yù)測時(shí)嚴(yán)格按公式Kte_c Kte - repmat(mean(K,1),mte,1) - repmat(mean(Kte,2),1,n) mean(K(:))/n計(jì)算其中mean(K,1)是訓(xùn)練核矩陣的列均值。5.4 原始數(shù)據(jù)有離群點(diǎn)前兩個(gè)潛變量被帶偏現(xiàn)象刪除兩個(gè)離群樣本后模型預(yù)測能力明顯提升殘差大幅下降。原因KPLS迭代第一步提取的是最大特征值對應(yīng)的特征向量離群樣本在核空間里遠(yuǎn)離數(shù)據(jù)云主體會顯著拉大核矩陣的距離項(xiàng)使得第一個(gè)潛變量花大量方差去擬合這個(gè)“遠(yuǎn)點(diǎn)”。解決訓(xùn)練前用簡單的3σ準(zhǔn)則或箱線圖剔除離群樣本更穩(wěn)妥的做法是訓(xùn)練后用HAT矩陣核矩陣的帽子矩陣診斷樣本杠桿值刪除杠桿值過大的樣本。5.5 訓(xùn)練集和測試集標(biāo)準(zhǔn)化參數(shù)不一致現(xiàn)象把訓(xùn)練集代碼復(fù)制一份改路徑后直接跑測試集結(jié)果預(yù)測殘差比預(yù)期大很多量級甚至超出訓(xùn)練集范圍。原因測試集單獨(dú)做了標(biāo)準(zhǔn)化均值和標(biāo)準(zhǔn)差是測試集自己的值和訓(xùn)練集不匹配。解決訓(xùn)練階段保存Xtr_mean、Xtr_std、Ytr_mean、Ytr_std四個(gè)標(biāo)量預(yù)測階段用同一組參數(shù)對Xte做標(biāo)準(zhǔn)化。6. 用交叉驗(yàn)證鎖定sigma和A再看殘差曲線確認(rèn)模型真沒白做KPLS程序跑通只是第一步參數(shù)整定才是決定模型好壞的關(guān)鍵。我自己的習(xí)慣是做兩層交叉驗(yàn)證外層10折評估泛化性能內(nèi)層在每一折里對sigma和A做網(wǎng)格搜索。A候選1~15sigma候選[0.1, 0.3, 1, 3, 10]每一組參數(shù)計(jì)算內(nèi)層驗(yàn)證集的均方根誤差RMSE選最小的一組。網(wǎng)格搜索的代碼框架非常簡單sigma_list [0.1, 0.3, 1, 3, 10]; A_list 1:15; best_sigma 1; best_A 5; best_rmse inf; for s sigma_list for a A_list % 這里嵌套5折交叉驗(yàn)證計(jì)算平均RMSE rmse cross_val_rmse(Xtr, Ytr, a, s); if rmse best_rmse best_rmse rmse; best_sigma s; best_A a; end end end參數(shù)鎖定的結(jié)果不能只看RMSE還要畫一條預(yù)測值與真實(shí)值的散點(diǎn)圖確認(rèn)點(diǎn)在45度線附近且沒有明顯的彎曲——如果散點(diǎn)在低值段聚攏、高值段發(fā)散說明模型對高值樣本的系統(tǒng)性偏差沒壓住這時(shí)優(yōu)先檢查是不是數(shù)據(jù)標(biāo)準(zhǔn)化時(shí)把高值段壓得太扁或者sigma偏小讓高值樣本落在核函數(shù)尾端。我在實(shí)際項(xiàng)目里還會畫潛變量得分T的前三列的三維散點(diǎn)看樣本在潛變量空間中是否有明顯的分簇趨勢這能幫助判斷是否還存在未建模的非線性結(jié)構(gòu)。最后說一個(gè)我養(yǎng)成的習(xí)慣每次訓(xùn)練完都把sigma、A、標(biāo)準(zhǔn)化參數(shù)、核矩陣K、回歸系數(shù)B一起存成MAT文件文件名帶數(shù)據(jù)版本號。KPLS調(diào)參的復(fù)現(xiàn)性比線性PLS敏感很多換一個(gè)sigma結(jié)果可能完全不同不保存現(xiàn)場就等于白跑。希望這篇筆記能幫你少走一段彎路——KPLS沒有想象中那么玄但每一步都值得確認(rèn)清楚。本文還有配套的精品資源點(diǎn)擊獲取