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