?一 . L0 范數(shù)正則化特征選擇

• 窮舉計(jì)算法

在面對變量選擇時(shí), 我們要進(jìn)行后驗(yàn)表示，我們另 ?rj = 1, 表示第j 個(gè)特征與此后驗(yàn)是相關(guān)的， 其中后驗(yàn)表達(dá)為

? ? ? ? ? ? ? ? ? ? ? ? ??

其中f(r) 為花費(fèi)函數(shù)， .?

假如有樣本 N = 20 ，維數(shù) D ?= 10，進(jìn)行線性回歸模型， 其中數(shù)據(jù)和噪聲為正太分布的，， 我們一般會要 求K 稀疏，表示稀疏的程度。 如K = 5，表示 有5個(gè)w 為非0， 即我們用 w = (0, -1.67, 0.13, 0, 0, 1.19, 0, -0.04, 0.33, 0) ， 并且噪聲方差滿足 ?。

?其實(shí)我們可以枚舉 2^10 = 1024 種模型來計(jì)算p(r|D)， 通過組合所有的特征來枚舉,最終得到最大的八組模型為：

?

進(jìn)行大量的模型枚舉判斷真的不太容易， 我們因而考慮一個(gè)自然而然的后驗(yàn)?zāi)Ｐ?MAP 估計(jì), 似然加上先驗(yàn)，?

? ? ? ??

在實(shí)際計(jì)算時(shí)，通常不會計(jì)算所有情況的概率密度，而是利用中值模型：

?

這就需要計(jì)算后驗(yàn)經(jīng)驗(yàn)概率（marginal inclusion probabilities） 。

我們可以分別得到每個(gè)特征對應(yīng)的后驗(yàn)概率情況：

因而可以通過改變?閾值來進(jìn)行特征的選擇添加過程。

如果上述 MAP 估計(jì)或是?marginal inclusion probabilities 都無效的話，就要考慮算法加速。

?首先先對上述計(jì)算進(jìn)行模型實(shí)例化：

The spike and slab model

后驗(yàn)的公式為：，?即在進(jìn)行MAP后驗(yàn)估計(jì) 是由先驗(yàn)和似然乘積決定。
對于先驗(yàn)，用以下先驗(yàn)對于位（bit:0 or 1）向量表示 ?伯努力： ? ? ? ? ? ? ? ? ? ? ??
其中是相關(guān)特征的概率，??， ||r||0 是 L0 的罰項(xiàng)的模值， 即非零特征的個(gè)數(shù)，
寫出其log 形式，從而化為 更似 l0正則的樣子：?
? ? ? ? ? ? ? ? ? ? ??
其中 lambda 控制模型的稀疏度。 ? ? ? ? ? ? ? ? ? ? ? ??
先驗(yàn)的函數(shù)與l0 掛上鉤了。下面對于似然函數(shù)? 這里對于符號的簡化，我們假設(shè)響應(yīng)y,是中心化過的（ie. ）,因此我們省慮了均值 u.
我們首先討論 p(w | r, ?sigma^2) ?中間那項(xiàng)， ?如果特征項(xiàng)?rj = 0， 則不相關(guān)有 wj = 0. ?如果 rj = 1, 我們希望 wj 是非零的。 ?如果我們對輸入數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化， 我們有一個(gè) 合適的先驗(yàn)（reasonable prior）, 即， ?其中 ?， 控制系數(shù)w與相關(guān)特征（rj =1）變量 的期望的大小（即不為零的wj 為大多呀，這個(gè)期望的值的控制程度）， ?這個(gè)項(xiàng)通過 sigma^2 來進(jìn)行scaled。? ? 因而似然的先驗(yàn)有： 第一項(xiàng)在原點(diǎn)處，即" spike", 當(dāng)? 趨于 無窮時(shí)， 這時(shí) p(wj | rj =1) 接近均勻分布， 即作為"slab", 因而叫做 spike and slab model.
我們可以將 wj = 0時(shí)的系數(shù)為零的情況去掉， 那時(shí)為0 系數(shù)。 這樣我們可化似然為：
? ? ? ? ?

其中 Dr = ||r||0, 是 r 中非0 項(xiàng)的個(gè)數(shù)， 通常簡化這個(gè)先驗(yàn)形式為:?, ?其中? 是任意正定矩陣。
在這樣的似然先驗(yàn)下， 我們可以計(jì)算 邊緣似然（marginal likelihood）。
如果我們知道噪聲的方差 sigma, 的話，那么邊緣似然為： ? ? ? ? ??? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??


如果不知道噪聲方差 sigma,我們可以放一個(gè) 噪聲先驗(yàn)，將噪聲項(xiàng)分離出來， 通常用?。 ?a, b的確定可參考（kohn et al. 2001） ， 當(dāng)我們分離出噪聲項(xiàng)，這樣邊緣似然為： ? ? ? ? ? ?

其中 S（r） 是 RSS , 殘差平方和： ? ? ? ? ? ??
邊緣似然 有的時(shí)候不能用這種近似的形式，比如logistic 回歸，或者非線性模型， 我們就用BIC 進(jìn)行近似，其形式如下： ? ? ? ? ? ? ?
其中??是ML or MAP 估計(jì) 基于 Xr, 這時(shí) ||r||0 是 這個(gè)模型的自由度“degrees of freedom”. 這時(shí)添加上前面的先驗(yàn)，得到目標(biāo)函數(shù)為： ? ? ? ? ? ? ?
其中兩個(gè)lambda 可以結(jié)合為一個(gè)，這樣就得到了l0 正則的形式 .
? ? ? ? ? 2.From the Bernoulli-Gaussian mode to L0 regularization 這個(gè)模型是通常使用的模型，其中參數(shù)分布為： ? ? ? ? ? ? ? ?

也叫做二值掩碼模型， 因?yàn)榭捎?rj 變量 來 掩除 wj 權(quán)重值， 使其為0. ?不像 1. 中 尖峰和 平坡模型中， 我們不能解釋那些不相關(guān)的系數(shù)， 這些系數(shù)是實(shí)際存在的， 二值掩碼模型 是通過 rj->y <-wj, ?而尖峰與平坡模型是 rj->wj->y, ?在此二值掩碼模型中， 只有 內(nèi)積項(xiàng)rjwj 能通過似然 估計(jì)出來。 通過大量用自己的選擇來推得目標(biāo)函數(shù)。? 首先聯(lián)合先驗(yàn)的形式為:
因此 非正太 -log 后驗(yàn)形式為：

其中? ， 我們需要 對w 進(jìn)行分撥離析， 即 和， 分別對應(yīng) r 的零項(xiàng)和非零項(xiàng)， 因?yàn)?X(r.*w)=XrWr, 因此我們使 = 0. 當(dāng)? 趨于 無窮時(shí) ， 。我們不需要考慮非零的w 的正則項(xiàng)，(這時(shí)因?yàn)闆]有復(fù)雜的罰來自邊緣似然或BIC近似)，這時(shí)目標(biāo)函數(shù)變?yōu)?
? ? ? ? ? ?? 我們現(xiàn)在不跟蹤每個(gè)r ， 用一個(gè)替換的相關(guān)變量來支持，或者w的非零項(xiàng)， 然后我們可以重寫目標(biāo)函數(shù)為： ? ? ? ? ? ??
這就是l0正則了， ?因?yàn)?, 因此此函數(shù) f(w) 非常的不平滑， 因此很難去優(yōu)化，我們需要找方法來優(yōu)化。

• 加速優(yōu)化方法

得到了L0正則化的目標(biāo)函數(shù)， 我們利用最小二乘以及貪心前進(jìn)尋找優(yōu)化解： 1. Single best replacement(最佳單特征替換) 也叫爬山法， ?開始為空， 每日有選擇特征 r = 0, ?然后單獨(dú)對每個(gè)特征進(jìn)行計(jì)算，找到一個(gè)最優(yōu)的， 之后對當(dāng)前的模型的鄰域?yàn)榉D(zhuǎn)r 的一個(gè)比特， 如果鄰域比特不在模型中，將其加入計(jì)算，若在模型中，將其刪除，計(jì)算， 這樣一直添加、刪除，直到模型沒有了提高。
2. 正交最小二乘法? ?如果SBR 不允許刪除，就退退化到此方法，過程如圖;

這個(gè)算法的缺點(diǎn)就是需要的計(jì)算量很大！ 3. Orthogonal matching persuits(正交投影尋蹤法) 這個(gè)算法一般在L0正則化時(shí)會經(jīng)常應(yīng)用，因此這里對這個(gè)算法詳細(xì)的介紹了解下。?

• 首先引入 MP 沒有正交的情況

在一個(gè)完備集合空間中 D = {xi}, 構(gòu)成哦 Hilbert 空間中， 其中 V = Span{xn}, v 是D 張成的那個(gè)空間。 這里這個(gè)D中的各個(gè)方向xi， 并非相互垂直的。 其中??(in H)， ?這里我們假設(shè)這里的 基底 xi， 都是已經(jīng)單位向量即||xn||= 1. Matching pursuit ?即做投影， 然后匹配， 求誤差， 添加 基底 利于表示，然后再投影..... ? ? ? ? ? ?
其中 f 投影之后的表示 為 : 近似 + 殘差? ? fk 是現(xiàn)在的近似， Rkf 是 現(xiàn)在的殘差， 在我們進(jìn)行迭代的時(shí)候，我們首先要進(jìn)行初始化， ?并且 k =1 ，開始進(jìn)行迭代
（I） ?首先計(jì)算 殘差與剩余各項(xiàng)的殘差，然后找到最大的進(jìn)行添加進(jìn)來 ? （II） 找到最大的內(nèi)積項(xiàng)，添加進(jìn)來 ：
（III） 添加進(jìn)來之后，我們就利用這新添加的特征向量，來進(jìn)行更新： 將此殘差最大的那項(xiàng)添加到 f 的近似中，這樣 f ，將更加接近實(shí)際值， 而 殘差在這一新添加項(xiàng)的的部分將被拿掉，給了 fk 作為近似使用了。 ? ? ? ? ?
（IV） 更新k , k<-k+1, 重復(fù)上述（I）-(IV)，直到收斂。
收斂性證明： 很重要的一個(gè)條件是滿足? ? ? ??? ?即在原來的Rkf 殘差中，將最大的那個(gè)?內(nèi)積已經(jīng)去掉，因而 這兩項(xiàng)現(xiàn)在滿足正交的結(jié)果， 因而由能量守恒有：
這樣我們每次迭代，都是殘差減小， 直到收斂。?
缺點(diǎn): MP 缺點(diǎn)即 需要進(jìn)行迭代的次數(shù)很多，因而進(jìn)行緩慢，?
每次迭代添加的項(xiàng)構(gòu)成了我們的f , ??, ?這樣雖然 xn1, xn2.....xnn 是一個(gè)張成的空間， fN 在這個(gè)上面的投影來進(jìn)行不斷的添加， 但是 只有當(dāng)?，滿足是我們的每次迭代才是 最優(yōu)近似的迭代， 而我們現(xiàn)在只滿足 ：， 即指正交與新添加項(xiàng)， 這樣
我們下次迭代式還可能對已經(jīng)出現(xiàn)過的項(xiàng)進(jìn)行重新迭代，這樣效率大大降低了。 因而我們得到的結(jié)果是次優(yōu)的。

舉個(gè)栗子：

如果 ?空間D = {x1, x2}，?且有 圖表示 f 與 x1, x2 的關(guān)系為：

?MP算法迭代會發(fā)現(xiàn)總是在 x1 和 x2 上反復(fù)迭代，即，這個(gè)就是信號(殘值)在已選擇的原子進(jìn)行垂直投影的非正交性導(dǎo)致的。?

這樣就引進(jìn)了 OMP 正交投影尋蹤法。

• OMP 娓娓道來了：

有了上面 MP 的了解之后， 這里就容易多了：

?每次迭代都保證 殘差與 現(xiàn)有的空間正交：?， 這里OMP 準(zhǔn)確的是兩個(gè)迭代，保證向后的計(jì)算是在正交情況下進(jìn)行的。

首先我們假設(shè)已經(jīng)有了 第 k 階 的結(jié)果 ：?

? ? ? ? ? ??

我們要做的是 進(jìn)行更新，k->k+1， ?注意的是 這里的?， ?第 k 次的殘差與 已添加所有的特征向量都正交。

k+1 : ? ?

? ? ? ? ? ???

這里的 k+1 次迭代情況， 與上述第 k 次 滿足同樣的條件。

兩者做差有：?

因?yàn)樵贒 中的各個(gè)特征向量 不是正交的， 因此我們構(gòu)造了一個(gè)輔助方程 來表示?, 與前 xn 的關(guān)系，?在Span(x1, x2..... ,xn)上的正交投影操作， 后面的一項(xiàng)是殘差。?。 （對的不要理解錯(cuò)誤了，是殘差與已經(jīng)添加的特征向量正交，而向量之間不正交。）?

? ? ? ? ? ? ??

?關(guān)系表示：?

將??表示的式子 ，帶入上面的做差的式子：?

? ? ? ? ? ? ??

如果此式的系數(shù)項(xiàng)和常數(shù)項(xiàng)分別滿足 ，則此公式一定滿足 。

? ? ? ? ? ? ? ?

我們令：?， 有 上兩式變?yōu)?/span>

后面的式子 ，兩邊都對??做內(nèi)積有： ? ? ? ? ? ??

? ? ? ? ?

這樣就完成了得帶的參數(shù)求解的過程， 其中??, 可通過 計(jì)算得出，后續(xù)參考中的《OMP_Krishnaprasad.pdf》給出了 其計(jì)算過程，rk 可以通過求出的?，來求殘差。 這樣可以得到第k+1 步迭代的參數(shù)。

? ? ? ? ?收斂性證明： 因?yàn)闈M足?

又因?yàn)?與?相互正交，所以由殘差構(gòu)成哦勾股定理有?

? ? ? ?， 這樣只要經(jīng)過 N步（方向向量的個(gè)數(shù)），就可全部迭代完畢。

算法的過程：

以上的算法的理解是建立在參數(shù)的理解過程的。 另外更簡單的方法是 直接選特征（最大內(nèi)積），然后用最小二乘直接計(jì)算參數(shù)，從而方面很多。

二、L0正則化實(shí)例

Prostate數(shù)據(jù)集上回歸分析

數(shù)據(jù)集介紹

該數(shù)據(jù)集包括了97位準(zhǔn)備做前列腺根治手術(shù)病人的前列腺特殊抗原（lpsa:?log?prostate?specific?antigen）和8個(gè)臨床指標(biāo)：

lcavol：log?cancel?volume?（腫瘤體積）

lweight：log?prostate?weight?（前列腺重量）

age：（年齡）

lbph：log?bengin?prostatic?hypcrplasia?（良性前列腺增生量）

svi：seminal?vesicle?invasion?（精囊浸潤）

lcp：log?of?capsular?penetration?（包膜穿透）

gleason：gleason?score?（Gleason積分）

pgg45：percent?of?Gleason?scores?4?or?5??（?Gleason4/5所占百分比?）。

2. 現(xiàn)在要做用L0正則做線性回歸：

用L0 就是對自己進(jìn)行選擇的過程，

?方法一： 窮舉搜索：

所有可能的模型共有 個(gè)。對每個(gè)模型計(jì)算其最小二乘結(jié)果，并比較每個(gè)模型的訓(xùn)練誤差和測試誤差。

首先讀取數(shù)據(jù)集：?

? ? istrain: [1x97 logical]
? ? ? ? ? y: [97x1 double]
? ? ? names: {1x9 cell}
? ? ? ? ? X: [97x8 double]
? ? ? ytest: [30x1 double]
? ? ?Xtrain: [67x8 double]
? ? ?ytrain: [67x1 double]
? ? ? Xtest: [30x8 double]

讀取數(shù)據(jù)，然后送去窮舉：

k = load('prostate.mat') % from prostateDataMake[n d] = size(Xtrain); [w, mseTrain, mseTest, sz, members] = allSubsetsRegression(Xtrain, ytrain, Xtest, ytest, 0:d, 1);
窮舉的過程，首先是標(biāo)準(zhǔn)化數(shù)據(jù)：

if doStandardize %對訓(xùn)練集x進(jìn)行標(biāo)準(zhǔn)化[Xtrain, mu] = center(Xtrain);[Xtrain, s] = mkUnitVariance(Xtrain);if ~isempty(Xtest)Xtest = center(Xtest, mu);Xtest = mkUnitVariance(Xtest, s);end end
之后標(biāo)準(zhǔn)化結(jié)束后，要得到所有的組合 256種，

ndx = ind2subv( 2*ones(1,d),1:(2^d) ) -1; %256×8 是個(gè)矩陣,子矩陣的排列組合。 NN = size(ndx,1);
之后開始進(jìn)行枚舉， 并記錄下每一次的計(jì)算的結(jié)果：

j = 1; for i=1:size(ndx,1) include = find(ndx(i,:));if ~ismember(length(include), allowableSizes)%fprintf('skipping size %d\n', length(include))continueelse%fprintf('including %s\n', sprintf('%d,', include));endmembers{j} = include;sz(j) = length(include);w{j} = [ones(n,1) Xtrain(:,include)]\ytrain; ypredTrain = [ones(n,1) Xtrain(:,include)]*w{j};RSS = sum((ytrain-ypredTrain).^2);mseTrain(j) = RSS/n;if ~isempty(Xtest)ntest = size(Xtest, 1);ypredTest = [ones(ntest,1) Xtest(:,include)]*w{j};mseTest(j) = mean((ypredTest-ytest).^2);elsemseTest = [];endj = j + 1; end
分別找到不同個(gè)數(shù)特征下最好的情況，然后找到測試結(jié)果最好的情況。

for i=0:dndx = find(sz==i);[besTrainScore(i+1) bestTrainSet] = min(mseTrain(ndx));[besTestScore(i+1) bestTestSet] = min(mseTest(ndx)); end[bestTestMSE ndx_star] = min(mseTest); w_star = w(ndx_star); model_star = members(ndx_star); bestTrainMSE = mseTrain(ndx_star); bestTrainMSE bestTestMSEfigure(1);clf plot(0:d,besTrainScore, 'bo-') hold on plot(0:d,besTestScore, 'ro-')%the selected/best model [bestTestMSE kstar] = min(besTestScore); vline(kstar-1,'r-.'); %將頂點(diǎn)連線畫出來xlabel('subset size') legend('best RSS on training set', 'best RSS on testing set'); title('all subsets on prostate cancer')
效果圖：

這里的測試誤差是模型在30個(gè)測試樣本上的真正測試誤差，不是用交叉驗(yàn)證估計(jì)的結(jié)果。訓(xùn)練誤差隨著模型復(fù)雜度的增加，一直下降；但測試誤差開始隨模型復(fù)雜度增加而下降，當(dāng)模型達(dá)到一定復(fù)雜度（特征數(shù)目為3）時(shí)，測試誤差隨著模型復(fù)雜度的增加而增加。所以當(dāng)模型復(fù)雜度高時(shí)，模型的訓(xùn)練誤差很小，但測試誤差較大，這就發(fā)生了過擬合現(xiàn)象。如果以真實(shí)測試誤差為標(biāo)準(zhǔn)，則我們選擇的最佳模型為子集大小為k=3。

最佳結(jié)果模型的參數(shù)表為：

該模型的訓(xùn)練誤差為0.6309，在測試集上的真實(shí)測試誤差為0.3828。

方法二： BIC 準(zhǔn)則

如果不是用真正的測試誤差，而是用BIC準(zhǔn)則作為測試誤差的估計(jì)的話，得到的最佳模型為子集大小為k=3。 其中BIC 判別分?jǐn)?shù)為：

? ? ? ? ?

這里若噪聲的方差未知，可用其極大似然估計(jì)代替， dof(M) 為模型M 中的自由參數(shù)的數(shù)目， 在線性回歸中為特征數(shù)目。

首先讀取數(shù)據(jù)：

load('prostate.mat') % from prostateDataMake[n d] = size(Xtrain); [w, mseTrain, mseTest, mseTest_BIC, sz, members] = allSubsetsRegression_BIC(Xtrain, ytrain, Xtest, ytest, 0:d, 1); 標(biāo)準(zhǔn)化：

if doStandardize[Xtrain, mu] = center(Xtrain);[Xtrain, s] = mkUnitVariance(Xtrain);if ~isempty(Xtest)Xtest = center(Xtest, mu);Xtest = mkUnitVariance(Xtest, s);end end
然后似然估計(jì) 噪聲方差

%OLS = MLE, estimate the variance of the noise, which will be used in compute %BIC beta_LS = [ones(n,1) Xtrain]\ytrain; ypredTrain = [ones(n,1) Xtrain]*beta_LS; % Predicted responses at each data point. RSS = sum((ytrain-ypredTrain).^2); % Residuals. sig2_hat = RSS/(n-d-1);
特征組合 并窮舉記錄結(jié)果：

ndx = ind2subv( 2*ones(1,d),1:(2^d) ) -1; NN = size(ndx,1); j = 1; for i=1:size(ndx,1)include = find(ndx(i,:));if ~ismember(length(include), allowableSizes)%fprintf('skipping size %d\n', length(include))continueelse%fprintf('including %s\n', sprintf('%d,', include));endmembers{j} = include;sz(j) = length(include);w{j} = [ones(n,1) Xtrain(:,include)]\ytrain; ypredTrain = [ones(n,1) Xtrain(:,include)]*w{j};RSS = sum((ytrain-ypredTrain).^2);mseTrain(j) = RSS/n;%BIC estimation of mseTestmseTest_BIC(j) = mseTrain(j) + log(n)*sz(j)*sig2_hat/n; % BIC; %testif ~isempty(Xtest)ntest = size(Xtest, 1);ypredTest = [ones(ntest,1) Xtest(:,include)]*w{j};mseTest(j) = mean((ypredTest-ytest).^2);elsemseTest = [];endj = j + 1; end
選出最小誤差結(jié)果，并畫出來：

for i=0:dndx = find(sz==i);[besTrainScore(i+1) bestTrainSet] = min(mseTrain(ndx));[besTestScore(i+1) bestTestSet] = min(mseTest(ndx));[besTestScore_BIC(i+1) bestTestSet_BIC] = min(mseTest_BIC(ndx)); %最小的BIC 的結(jié)果 end[bestBIC ndx_star] = min(mseTest_BIC); w_star = w(ndx_star); model_star = members(ndx_star); bestmseTrain = mseTrain(ndx_star); bestmseTest = mseTest(ndx_star); bestBIC bestmseTrain bestmseTestfigure clf plot(0:d,besTestScore_BIC, 'mo-') hold on plot(0:d,besTrainScore, 'bo-') hold on plot(0:d,besTestScore, 'ro-')%the selected/best model [bestBIC kstar] = min(besTestScore_BIC);hold on vline(kstar-1,'m-.');xlabel('subset size') ylabel('error')legend('best BIC score','best RSS on training set', 'best RSS on testing set'); title('all subsets on prostate cancer BIC')
結(jié)果圖：

其中參數(shù)結(jié)果為：

該模型的訓(xùn)練誤差為0.5210，在測試集上的真實(shí)測試誤差為0.4815。

方法三： 交叉驗(yàn)證

如果用10折交叉驗(yàn)證測試誤差估計(jì)，并采用1倍標(biāo)準(zhǔn)差原則的話，得到的最佳模型為子集大小為k=2。

首先讀取數(shù)據(jù)：

load('prostate.mat') % from prostateDataMake[n d] = size(Xtrain); [w, mseTrain, mseTest,model] = allSubsetsRegression_CV(Xtrain, ytrain, Xtest, ytest, 0:d, 1);
標(biāo)準(zhǔn)化數(shù)據(jù)：

%標(biāo)準(zhǔn)化數(shù)據(jù) if doStandardize[Xtrain, mu] = center(Xtrain);[Xtrain, s] = mkUnitVariance(Xtrain);if ~isempty(Xtest)Xtest = center(Xtest, mu);Xtest = mkUnitVariance(Xtest, s);end end
將數(shù)據(jù)分成10折， 比如100 個(gè)數(shù)據(jù)，第一次[0-10] test, [11-100] train, 第二次[11-20] test, [0-10],[21-100] train .... 等等

%cross validation Nfolds = 10; seed= 0; rand('state', seed);sizes = 0:d; perm = randperm(n); xtrainAll = Xtrain(perm,:); ytrainAll = ytrain(perm); [trainfolds, testfolds] = Kfold(size(xtrainAll,1), Nfolds); %10 folds 交叉驗(yàn)證

分別選每折數(shù)據(jù)進(jìn)行， 組合 窮舉： clear mseCVTrain mseCVTest for f=1:length(trainfolds)xtrainCV = xtrainAll(trainfolds{f},:);ytrainCV = ytrainAll(trainfolds{f});xtestCV = xtrainAll(testfolds{f},:);ytestCV = ytrainAll(testfolds{f});[w, mseCVTrainAll(f,:), mseCVTestAll(f,:), sz, members, df, mseCVTest(f,:)] = ...allSubsetsRegression(xtrainCV, ytrainCV, xtestCV, ytestCV, sizes, 0); end
算出每個(gè)折的一倍標(biāo)準(zhǔn)差，以及每折得到的對應(yīng)的個(gè)數(shù)求出的均方差的平均。 %mseCVTest 這個(gè)得到的是那個(gè)最好的測試的結(jié)果，就是每折之后， 選出的[], 1，2 ，3 ...,8 枚舉過的最小的均方誤差。 mseCVMean = mean(mseCVTest,1); %每一折的均方差的平均 ， 即列的平均啊， 最小的那個(gè)測試標(biāo)準(zhǔn)偏差。 mseCVse = std(mseCVTest,[],1)/sqrt(Nfolds); %一倍標(biāo)準(zhǔn)差， 標(biāo)準(zhǔn)偏差， 1×9
然后是其一倍標(biāo)準(zhǔn)差表示， 以及找到在變量最少（模型最簡單）的情況下， 滿足最小均方差一倍標(biāo)準(zhǔn)差范圍的。 figure; %誤差帶的上下各為 mseCVse errorbar(df, mseCVMean, mseCVse); % confidence intervals of data or the deviation along a curve. kstar = oneStdErrorRule(mseCVMean, mseCVse); %找出根據(jù)均值和一倍標(biāo)準(zhǔn)差， 最好的一個(gè)滿足的。， 因?yàn)樽钚〉酱笃涮卣鲾?shù)在增加，因而其復(fù)雜度在增加。 hold on vline(kstar,'r-.'); xlabel('size of subset') ylabel('cv error') title(sprintf('all subsets CV'))
通過得到這個(gè)最少變量得個(gè)數(shù)，在這個(gè)最小變量個(gè)數(shù)下能表示 此模型的情況。 % Refit using all training data to find set of chosen size [w, mseTrain, junk, junk2, members] = ...allSubsetsRegression(Xtrain, ytrain, [], [], sizes(kstar), 0); %所以找到兩個(gè)的最佳的那個(gè)w， [junk,best]=min(mseTrain); w = w(best); model = members{best}; mseTrain = mseTrain(best);modelif ~isempty(Xtest)ntest = size(Xtest, 1);ypredTest = [ones(ntest,1) Xtest(:,model)]* cell2mat(w);mseTest = mean((ypredTest-ytest).^2);elsemseTest = 0; endmseTrain mseTest



此時(shí)模型的參數(shù)結(jié)果為：

該模型的訓(xùn)練誤差為0.?5536，在測試集上的真實(shí)測試誤差為0.5737。

方法四：正交標(biāo)準(zhǔn)尋蹤

我們以BIC準(zhǔn)則估計(jì)模型的測試誤差，采用正交投影尋蹤每次增加一個(gè)與當(dāng)前殘差最相關(guān)的特征，然后更新模型參數(shù)。前向逐步回歸選擇的模型為子集大小為4；當(dāng) 特征數(shù)目變成5時(shí)，BIC分?jǐn)?shù)增加，迭代終止。

讀取數(shù)據(jù)：

load('prostate.mat') % from prostateDataMake [n d] = size(Xtrain); [w, mseTrain, mseTest, mseTest_BIC, sz, members] = OMP_BIC(Xtrain, ytrain, Xtest, ytest, 0:d, 1);
數(shù)據(jù)標(biāo)準(zhǔn)化：

if doNormalized[Xtrain, mu] = center(Xtrain);[Xtrain, norm2] = normalize(Xtrain);if ~isempty(Xtest)Xtest = center(Xtest, mu);Xtest = normalize(Xtest, norm2);end end
然后最大似然求方差?

%OLS, for estimate the variance of the noise, which will be used in compute %BIC beta_LS = [ones(n,1) Xtrain]\ytrain; ypredTrain = [ones(n,1) Xtrain]*beta_LS; % Predicted responses at each data point. RSS = sum((ytrain-ypredTrain).^2); % Residuals. sig2_hat = RSS/n;
當(dāng)數(shù)據(jù)為空時(shí)， 最開始的初始化情況有：

% the empty model include = []; %包含的特征 j = 1; members{j} = include; sz(j) = length(include); w0 = mean(ytrain); w{j} = w0; ypredTrain = [ones(n,1) Xtrain(:,include)]* w{j} ; RSS = sum((ytrain-ypredTrain).^2); mseTrain(j) = RSS/n;%BIC estimation of mseTest mseTest_BIC(j) = mseTrain(j) + log(n)*sz(j)*sig2_hat/n; % BIC; bestBIC = mseTest_BIC(j);%test if ~isempty(Xtest)ntest = size(Xtest, 1);ypredTest = [ones(ntest,1) Xtest(:,include)]*w{j};mseTest(j) = mean((ypredTest-ytest).^2); elsemseTest = []; endndx_star = j;
然后進(jìn)入OMP 逐漸的迭代過程，利用BIC 分?jǐn)?shù)作為判別的標(biāo)準(zhǔn)，當(dāng)BIC 分?jǐn)?shù)增加時(shí)（誤差變大）迭代終止：

j = j+1; for i = 1 : d%當(dāng)前模型的殘差residual = ytrain-ypredTrain;%殘差與特征之間的相關(guān)性（投影）proj = Xtrain' * residual;%選擇與殘差最相關(guān)的特征[maxval, bestk] = max(abs(proj));include = [include bestk]; %add the mink featuremembers{j} = include; % add best k th feature tempsz(j) = length(members{j}); %least squares w{j} = [ones(n,1) Xtrain(:,members{j})]\ytrain; ypredTrain = [ones(n,1) Xtrain(:,members{j})]* w{j} ;RSS = sum((ytrain-ypredTrain).^2);mseTrain(j) = RSS/n;%BIC estimation of mseTestmseTest_BIC(j) = mseTrain(j) + log(n)*sz(j)*sig2_hat/n; % BIC; 利用此分?jǐn)?shù)來進(jìn)行終止條件判斷。%testif ~isempty(Xtest)ntest = size(Xtest, 1);ypredTest = [ones(ntest,1) Xtest(:,members{j})]* w{j} ;mseTest(j) = mean((ypredTest-ytest).^2);elsemseTest = [];endif mseTest_BIC(j) > bestBIC % the BIC tend to increasebreak;elsebestBIC = mseTest_BIC(j);endj = j + 1; end
結(jié)果：

三、總結(jié)

通過L0 正則化進(jìn)行特征選擇， 但是一般情況L0 正則導(dǎo)致函數(shù)不光滑， 優(yōu)化方法不容易計(jì)算， 因此，我們將繼續(xù)研究L1正則化， 作為L0 的 “最優(yōu)突近似”。?

四、參考文獻(xiàn)

1. <<?Machine Learning?A Probabilistic Perspective >>?

2.?Orthogonal Matching Pursuit:Recursive Function Approximat ion with Applications to Wavelet?Decomposition?

3.?機(jī)器學(xué)習(xí)中的范數(shù)規(guī)則化之（一）L0、L1與L2范數(shù)?-?http://blog.csdn.net/zouxy09/article/details/24971995

總結(jié)

以上是生活随笔為你收集整理的机器学习稀疏之L0正则化的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。