机器学习稀疏之L0正则化

最新推荐文章于 2021-07-26 09:20:17 发布

原创最新推荐文章于 2021-07-26 09:20:17 发布 · 4.8k 阅读

5 ·

本内容遵循CC 4.0 BY-SA版权协议

Machine Learning 专栏收录该内容

9 篇文章

订阅专栏

本文介绍了L0正则化在特征选择中的应用，包括穷举计算法和MAP估计，讨论了正交投影寻踪法（OMP）作为优化手段。通过Prostate数据集的案例分析，展示了L0正则化在实际问题中的运用，并对比了不同选择模型的方法，如BIC准则和交叉验证。最后，指出L0正则化的挑战及其与L1正则化的联系。

一 . L0 范数正则化特征选择

穷举计算法

在面对变量选择时, 我们要进行后验表示，我们另 rj = 1, 表示第j 个特征与此后验是相关的，其中后验表达为

$p(r|D) = \frac{{{e^{ - f(r)}}}}{{\sum\nolimits_{{r^'}} {{e^{ - f(r)}}} }}$

其中f(r) 为花费函数， $f(r) \buildrel \Delta \over = - [logp(D|r) + logp(r)]$ .

假如有样本 N = 20 ，维数 D = 10，进行线性回归模型，其中数据和噪声为正太分布的， ${y_i} \sim N({w^T}{x_i},{\sigma ^2})$ ，我们一般会要求K 稀疏，表示稀疏的程度。如K = 5，表示有5个w 为非0，即我们用 w = (0, -1.67, 0.13, 0, 0, 1.19, 0, -0.04, 0.33, 0) ，并且噪声方差满足 ${\sigma ^2} = 1$ 。

其实我们可以枚举 2^10 = 1024 种模型来计算p(r|D)，通过组合所有的特征来枚举,最终得到最大的八组模型为：

进行大量的模型枚举判断真的不太容易，我们因而考虑一个自然而然的后验模型 MAP 估计, 似然加上先验，

$\hat r = \arg \max p(r|D) = argminf(r)$

在实际计算时，通常不会计算所有情况的概率密度，而是利用中值模型：

$\hat r = \{ j:p({r_j} = 1|D) > 0.5\}$

这就需要计算后验经验概率（marginal inclusion probabilities） $p({r_j} = 1|D)$ 。

我们可以分别得到每个特征对应的后验概率情况：

因而可以通过改变 $p({r_j} = 1|D)$ 阈值来进行特征的选择添加过程。

如果上述 MAP 估计或是 marginal inclusion probabilities 都无效的话，就要考虑算法加速。

首先先对上述计算进行模型实例化：

The spike and slab model

后验的公式为： $p(r|D) \propto p(r)p(D|r)$ ，即在进行MAP后验估计是由先验和似然乘积决定。

对于先验，用以下先验对于位（bit:0 or 1）向量表示伯努力：

$p(r) = \prod\limits_{j = 1}^D {Ber({r_j}|{\pi _0}) = \pi _0^{||r|{|_0}}{{(1 - {\pi _0})}^{D - ||r|{|_0}}}}$

其中 ${\pi _0}$ 是相关特征的概率， $||r|{|_0} = \sum\nolimits_{j = 1}^D {{r_j}}$ ， ||r||0 是 L0 的罚项的模值，即非零特征的个数，

写出其log 形式，从而化为更似 l0正则的样子：

$\begin{align*} {\mathop{\rm logp}\nolimits} (r|{\pi _0}) &= ||r|{|_0}\log {\pi _0} + (D - ||r|{|_0})\log (1 - {\pi _0}) \\ &= {||r|{|_0}(log{\pi _0} - \log (1 - {\pi _0}))} + const \\ &= { - \lambda } ||r|{|_0} + const \end{align*}$

其中 lambda 控制模型的稀疏度。

$\lambda \buildrel \Delta \over = \log \frac{{1 - {\pi _0}}}{{{\pi _0}}}$

先验的函数与l0 挂上钩了。下面对于似然函数

$p(D|r) = p(y|X,r) = \int \int {p(y|X,w,r)p(w|r,{\sigma ^2})p({\sigma ^2})}dwd{\sigma^2}$

这里对于符号的简化，我们假设响应y,是中心化过的（ie. $\bar y = 0$ ）,因此我们省虑了均值 u.

我们首先讨论 p(w | r, sigma^2) 中间那项，如果特征项 rj = 0，则不相关有 wj = 0. 如果 rj = 1, 我们希望 wj 是非零的。如果我们对输入数据进行标准化，我们有一个合适的先验（reasonable prior）, 即 $N(0,{\sigma ^2}\sigma _w^2)$ ，其中 $\sigma _w^2$ ，控制系数w与相关特征（rj =1）变量的期望的大小（即不为零的wj 为大多呀，这个期望的值的控制程度），这个项 $\sigma _w^2$ 通过 sigma^2 来进行scaled。

因而似然的先验有：

$\[p({w_j}|{\sigma ^2},{r_j}) \begin{cases} {\delta _0}({w_j}) & {if }\; {r_j} = 0 \\ N({w_j}|0,{\sigma ^2}\sigma _w^2) &{if }\; {r_j} = 1 \end{cases}$

第一项在原点处，即" spike", 当 $\sigma _w^2$ 趋于无穷时，这时 p(wj | rj =1) 接近均匀分布，即作为"slab", 因而叫做 spike and slab model.

我们可以将 wj = 0时的系数为零的情况去掉，那时为0 系数。这样我们可化似然为：

$p(D|r) = \int \int {N({y}|{X_r}w_r,{\sigma ^2}I_{N})N({{w_j}|0_{{D_r}}},{\sigma ^2}\sigma _w^2{I_{{D_r}}})p({\sigma ^2})}dw_rd{\sigma^2}$

其中 Dr = ||r||0, 是 r 中非0 项的个数，通常简化这个先验形式为: $p(w|r,{\sigma ^2}) = N({w_r}|{0_{{D_r}}},{\sigma ^2}\sum\nolimits_r {} )$ , 其中 $\sum\nolimits_r {}$ 是任意正定矩阵。

在这样的似然先验下，我们可以计算边缘似然（marginal likelihood）。

如果我们知道噪声的方差 sigma, 的话，那么边缘似然为：

$p(D|r,{\sigma ^2}) = \int {N(y|{X_r}{w_r},{\sigma ^2}I)} N({w_r}|{0_{{D_r}}},{\sigma ^2}\sum\nolimits_r {} )d{w_r} = N(y|0,{C_r})$

${C_r} \buildrel \Delta \over = {\sigma ^2}{X_r}\sum\nolimits_r {X_{_r}^T} + {\sigma ^2}{I_N}$

如果不知道噪声方差 sigma,我们可以放一个噪声先验，将噪声项分离出来，通常用 $p({\sigma ^2}) = IG({\sigma ^2}|{a_\sigma },{b_\sigma })$ 。 a, b的确定可参考（kohn et al. 2001），当我们分离出噪声项，这样边缘似然为：

$\begin{align*} p(D|r) & = \int \int {p(y|r,{w_r},{\sigma ^2})p({w_r}|r,{\sigma ^2})p({\sigma ^2})d{w_r}}d{\sigma^2} \\ & \propto |X_{_r}^T{X_r} + \sum\nolimits_r^{ - 1} {} {|^{ - \frac{1}{2}}}|\sum\nolimits_r {} {|^{ - \frac{1}{2}}}{(2{b_\sigma } + S(r))^{ - (2{a_\sigma } + N - 1)/2}} \end{align*}$

其中 S（r）是 RSS , 残差平方和：

$S(r) \buildrel \Delta \over = {y^T}y - {y^T}{X_r}{(X_r^T{X_r} + \sum\nolimits_r^{ - 1} {} )^{ - 1}}X_{_r}^Ty$

边缘似然有的时候不能用这种近似的形式，比如logistic 回归，或者非线性模型，我们就用BIC 进行近似，其形式如下：

${\mathop{\rm logp}\nolimits} (D|r) \approx\ logp(y|X,{{\hat w}_r},\hat \sigma^2 )- \frac{{||r|{|_0}}}{2}\log N$

其中 ${{\hat w}_r}$ 是ML or MAP 估计基于 Xr, 这时 ||r||0 是这个模型的自由度“degrees of freedom”.

这时添加上前面的先验，得到目标函数为：

${\mathop{\rm logp}\nolimits} (r|D) \approx\ logp(y|X,{{\hat w}_r},\hat \sigma^2 )- \frac{{||r|{|_0}}}{2}\log N -\lambda{{||r|{|_0}}}+const$

其中两个lambda 可以结合为一个，这样就得到了l0 正则的形式 .

2.From the Bernoulli-Gaussian mode to L0 regularization

这个模型是通常使用的模型，其中参数分布为：

$\begin{align*} [y|{x_r},w,r,{\sigma ^2} &\sim N(\sum\limits_j {{r_j}{w_j}{x_{ij}},} {\sigma ^2}) \\ {r_j} & \sim Ber({\pi _0})\\ {w_j}& \sim N(0,\sigma _{_w}^2) \end{align*}$

也叫做二值掩码模型，因为可用 rj 变量来掩除 wj 权重值，使其为0.

不像 1. 中尖峰和平坡模型中，我们不能解释那些不相关的系数，这些系数是实际存在的，二值掩码模型是通过 rj->y <-wj, 而尖峰与平坡模型是 rj->wj->y, 在此二值掩码模型中，只有内积项rjwj 能通过似然估计出来。

通过大量用自己的选择来推得目标函数。

首先联合先验的形式为:

$p(r,w) \propto N(w|0,\sigma _{_w}^2I)\pi _{_0}^{||r|{|_0}}{(1 - {\pi _0})^{D - ||r|{|_0}}}$

因此非正太 -log 后验形式为：

$f(r,w) \buildrel \Delta \over = - 2{\sigma ^2}\log p(r,w,y|X) = ||y - X(r.*w)|{|^2} + \frac{{{\sigma ^2}}}{{\sigma _w^2}}||w|{|^2} +\lambda||r||{}_0+\text{const}$

其中

$\lambda \buildrel \Delta \over = 2{\sigma ^2}\log (\frac{{1 - {\pi _0}}}{{{\pi _0}}})$ ，我们需要对w 进行分拨离析，即 ${w_{ - r}}$ 和 ${w_{ r}}$ ，分别对应 r 的零项和非零项，因为 X(r.*w)=XrWr, 因此我们使 ${w_{ - r}}$ = 0.

当 $\sigma _w^2$ 趋于无穷时，。我们不需要考虑非零的w 的正则项，(这时因为没有复杂的罚来自边缘似然或BIC近似)，这时目标函数变为:

$f(r,w) = ||y - {X_r}{w_r}|{|^2} + \lambda ||r||{}_0$

我们现在不跟踪每个r ，用一个替换的相关变量来支持，或者w的非零项，然后我们可以重写目标函数为：

$f(w) = ||y - Xw|{|^2} + \lambda ||w||{}_0$

这就是l0正则了，因为 $r \in {\{ 0,1\} ^D}$ , 因此此函数 f(w) 非常的不平滑，因此很难去优化，我们需要找方法来优化。

加速优化方法

得到了L0正则化的目标函数，我们利用最小二乘以及贪心前进寻找优化解：

1. Single best replacement(最佳单特征替换)

也叫爬山法，开始为空，每日有选择特征 r = 0, 然后单独对每个特征进行计算，找到一个最优的，之后对当前的模型的邻域为翻转r 的一个比特，如果邻域比特不在模型中，将其加入计算，若在模型中，将其删除，计算，这样一直添加、删除，直到模型没有了提高。

2. 正交最小二乘法

如果SBR 不允许删除，就退退化到此方法，过程如图;

这个算法的缺点就是需要的计算量很大！

3. Orthogonal matching persuits(正交投影寻踪法)

这个算法一般在L0正则化时会经常应用，因此这里对这个算法详细的介绍了解下。

首先引入 MP 没有正交的情况

在一个完备集合空间中 D = {xi}, 构成哦 Hilbert 空间中，其中 V = Span{xn}, v 是D 张成的那个空间。这里这个D中的各个方向xi，并非相互垂直的。其中 $w = {v^ \bot }$ (in H)，这里我们假设这里的基底 xi，都是已经单位向量即||xn||= 1.

Matching pursuit 即做投影，然后匹配，求误差，添加基底利于表示，然后再投影.....

${P_v}f = \sum\limits_n {{a_n}{x_n}}$

其中 f 投影之后的表示为 : 近似 + 残差

$f = \sum\limits_{i = 1}^k {{a_i}} {x_{{n_i}}} + {R_k}f = {f_k} + {R_k}f$

fk 是现在的近似， Rkf 是现在的残差，在我们进行迭代的时候，我们首先要进行初始化， ${R_0}f = f,{f_0} = 0$ 并且 k =1 ，开始进行迭代

（I） 首先计算残差与剩余各项的残差，然后找到最大的进行添加进来

${\{ \left\langle {{R_k}f,{x_n}} \right\rangle \} _n}$

（II） 找到最大的内积项，添加进来：

$\begin{array}{l} \left| {\left\langle {{R_k}f,{x_{{n_{k + 1}}}}} \right\rangle } \right| \ge \alpha \mathop {\sup }\limits_j \left| {\left\langle {{R_k}f,{x_{{n_j}}}} \right\rangle } \right|\\ where{\kern 1pt} {\kern 1pt} 0 < \alpha \le 1 \end{array}$

（III） 添加进来之后，我们就利用这新添加的特征向量，来进行更新：将此残差最大的那项添加到 f 的近似中，这样 f ，将更加接近实际值，而残差在这一新添加项的的部分将被拿掉，给了 fk 作为近似使用了。

$\begin{array}{l} {f_{k + 1}} = {f_k} + \left\langle {{R_k}f,{x_{{n_{k + 1}}}}} \right\rangle {x_{{n_{k + 1}}}}\\ {R_{k + 1}}f = {R_k}f - \left\langle {{R_k}f,{x_{{n_{k + 1}}}}} \right\rangle {x_{{n_{k + 1}}}} \end{array}$

（IV） 更新k , k<-k+1, 重复上述（I）-(IV)，直到收敛。

收敛性证明：很重要的一个条件是满足

$\left\langle {{R_{k + 1}}f,{x_{{n_{k + 1}}}}} \right\rangle = 0$ 即在原来的Rkf 残差中，将最大的那个 ${{x_{{n_{k + 1}}}}}$ 内积已经去掉，因而这两项现在满足正交的结果，因而由能量守恒有：

$||{R_k}f|{|^2} = ||{R_{k + 1}}f|{|^2} + |\left\langle {{R_k}f,{x_{{n_{k + 1}}}}} \right\rangle {|^2}$

这样我们每次迭代，都是残差减小，直到收敛。

缺点: MP 缺点即需要进行迭代的次数很多，因而进行缓慢，

每次迭代添加的项构成了我们的f , ${f_N} = \sum\limits_{k = 0}^{N - 1} {\left\langle {{R_k}f,{x_{{n_{k + 1}}}}} \right\rangle } {x_{{n_{k + 1}}}}$ , 这样虽然 xn1, xn2.....xnn 是一个张成的空间， fN 在这个上面的投影来进行不断的添加，但是只有当 ${R_N}f \in V_N^ \bot$ ，满足是我们的每次迭代才是最优近似的迭代，而我们现在只满足： ${R_k}f \bot {x_{{n_N}}}$ ，即指正交与新添加项，这样

我们下次迭代式还可能对已经出现过的项进行重新迭代，这样效率大大降低了。因而我们得到的结果是次优的。

举个栗子：

如果空间D = {x1, x2}，且有图表示 f 与 x1, x2 的关系为：

MP算法迭代会发现总是在 x1 和 x2 上反复迭代，即 $f = {a_1}{x_1} + {a_2}{x_2} + {a_3}{x_1} + {a_4}{x_2} + ....$ ，这个就是信号(残值)在已选择的原子进行垂直投影的非正交性导致的。

这样就引进了 OMP 正交投影寻踪法。

OMP 娓娓道来了：

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

每次迭代都保证残差与现有的空间正交： ${R_N}f \in V_N^ \bot$ ，这里OMP 准确的是两个迭代，保证向后的计算是在正交情况下进行的。

首先我们假设已经有了第 k 阶的结果：

$f = \sum\limits_{n = 1}^k {a_{_n}^k{x_n} + {R_k}f,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} with} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\langle {{R_k}f,{x_n}} \right\rangle = 0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} n = 1,...,k.$

我们要做的是进行更新，k->k+1，注意的是这里的 $\left\langle {{R_k}f,{x_n}} \right\rangle = 0$ ，第 k 次的残差与已添加所有的特征向量都正交。

k+1 :

$f = \sum\limits_{n = 1}^{k + 1} {a_{_n}^{k + 1}{x_n} + {R_{k + 1}}f,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} with} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\langle {{R_{k + 1}}f,{x_n}} \right\rangle = 0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} n = 1,...,k + 1.$

这里的 k+1 次迭代情况，与上述第 k 次满足同样的条件。

两者做差有：

$\sum\limits_{n = 1}^{k + 1} {(a_n^{k + 1} - a_n^k){x_n}} + a_{k + 1}^{k + 1}{x_{k + 1}} + {R_{k + 1}}f - {R_k}f$

因为在D 中的各个特征向量不是正交的，因此我们构造了一个辅助方程来表示 ${x_{k + 1}}$ , 与前 xn 的关系， ${x_{k + 1}}$ 在Span(x1, x2..... ,xn)上的正交投影操作，后面的一项是残差。。（对的不要理解错误了，是残差与已经添加的特征向量正交，而向量之间不正交。）

${x_{k + 1}} = \sum\limits_{n = 1}^k {b_n^k{x_n}} + {r_k},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} with{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\langle {{r_k},{x_n}} \right\rangle = 0,{\kern 1pt} {\kern 1pt} {\kern 1pt} n = 1,...,k$

关系表示： $\sum\limits_{n = 1}^k {b_n^k{x_n}} = {P_{{V_k}}}{x_{k + 1}},and{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {r_k} = {P_{V_k^ \bot }}{x_{k + 1}}$

将 ${x_{k + 1}}$ 表示的式子，带入上面的做差的式子：

$\sum\limits_{n = 1}^k {(a_n^{k + 1} - a_n^k + a_{k + 1}^{k + 1}b_n^k){x_n}} + (a_{k + 1}^{k + 1}{r_k} + {R_{k + 1}}f - {R_k}f) = 0$

如果此式的系数项和常数项分别满足，则此公式一定满足。

$\begin{align*} {a_n^{k + 1} - a_n^k + a_{k + 1}^{k + 1}b_n^k} & =0 \\ a_{k + 1}^{k + 1}{r_k} + {R_{k + 1}}f - {R_k}f &= 0 \end{align*}$

我们令： $a_{k + 1}^{k + 1} = {a_k}$ ，有上两式变为

$a_n^{k + 1} = a_n^k + {a_k}b_n^k,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} n = 1,....,k$

后面的式子，两边都对 ${x_{k + 1}}$ 做内积有：

${a_k} = \frac{{\left\langle {{R_k}f,{x_{k + 1}}} \right\rangle }}{{\left\langle {{r_k},{x_{k + 1}}} \right\rangle }} = \frac{{\left\langle {{R_k}f,{x_{k + 1}}} \right\rangle }}{{||{r_k}|{|^2}}} = \frac{{\left\langle {{R_k}f,{x_{k + 1}}} \right\rangle }}{{||{x_{k + 1}}|{|^2} - \sum\nolimits_{n = 1}^k {b_n^k\left\langle {{x_k},{x_{k + 1}}} \right\rangle } }}$

这样就完成了得带的参数求解的过程，其中 $b_n^k$ , 可通过计算得出，后续参考中的《OMP_Krishnaprasad.pdf》给出了其计算过程，rk 可以通过求出的 $b_n^k$ ，来求残差。这样可以得到第k+1 步迭代的参数。

收敛性证明：因为满足 ${R_k}f = {R_{k + 1}}f + {a_k}{r_k}$

又因为 ${r_k}$ 与 $R_{k+1}f$ 相互正交，所以由残差构成哦勾股定理有

$||{R_k}f|{|^2} = ||{R_{k + 1}}f|{|^2} + \frac{{|\left\langle {{R_k}f,{x_{{n_{k + 1}}}}} \right\rangle {|^2}}}{{||{r_k}|{|^2}}}$ ，这样只要经过 N步（方向向量的个数），就可全部迭代完毕。

算法的过程：

以上的算法的理解是建立在参数的理解过程的。另外更简单的方法是直接选特征（最大内积），然后用最小二乘直接计算参数，从而方面很多。

二、L0正则化实例

Prostate数据集上回归分析

数据集介绍

该数据集包括了97位准备做前列腺根治手术病人的前列腺特殊抗原（lpsa: log prostate specific antigen）和8个临床指标：

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. 现在要做用L0正则做线性回归：

用L0 就是对自己进行选择的过程，

方法一：穷举搜索：

所有可能的模型共有 ${2^D} = {2^8} = 256$ 个。对每个模型计算其最小二乘结果，并比较每个模型的训练误差和测试误差。

首先读取数据集：

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

读取数据，然后送去穷举：

k = load('prostate.mat') % from prostateDataMake

[n d] = size(Xtrain);
[w, mseTrain, mseTest, sz, members] =  allSubsetsRegression(Xtrain, ytrain, Xtest, ytest, 0:d, 1);

穷举的过程，首先是标准化数据：

if doStandardize   %对训练集x进行标准化
  [Xtrain, mu]  = center(Xtrain);
  [Xtrain, s]  = mkUnitVariance(Xtrain);
  if ~isempty(Xtest)
    Xtest = center(Xtest, mu);
    Xtest = mkUnitVariance(Xtest, s);
  end
end

之后标准化结束后，要得到所有的组合 256种，

ndx = ind2subv( 2*ones(1,d),1:(2^d) ) -1;   %256×8 是个矩阵,子矩阵的排列组合。
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))
    continue
  else
    %fprintf('including %s\n', sprintf('%d,', include));
  end

  members{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);
  else
    mseTest = [];
  end
  
  j = j + 1;
end

分别找到不同个数特征下最好的情况，然后找到测试结果最好的情况。

for i=0:d
  ndx = 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
bestTestMSE

figure(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-.');   %将顶点连线画出来

xlabel('subset size')
legend('best RSS on training set', 'best RSS on testing set');
title('all subsets on prostate cancer')

效果图：

这里的测试误差是模型在30个测试样本上的真正测试误差，不是用交叉验证估计的结果。训练误差随着模型复杂度的增加，一直下降；但测试误差开始随模型复杂度增加而下降，当模型达到一定复杂度（特征数目为3）时，测试误差随着模型复杂度的增加而增加。所以当模型复杂度高时，模型的训练误差很小，但测试误差较大，这就发生了过拟合现象。如果以真实测试误差为标准，则我们选择的最佳模型为子集大小为k=3。

最佳结果模型的参数表为：

该模型的训练误差为0.6309，在测试集上的真实测试误差为0.3828。

方法二： BIC 准则

如果不是用真正的测试误差，而是用BIC准则作为测试误差的估计的话，得到的最佳模型为子集大小为k=3。其中BIC 判别分数为：

$BIC(M) = \frac{1}{N}\sum\limits_{i = 1}^N {{{({y_i} - {{\hat w}^T}{x_i})}^2} + (logN)\frac{{dof(M)}}{N}{\sigma ^2}}$

这里若噪声的方差未知，可用其极大似然估计 ${{\hat \sigma }^2}$ 代替， dof(M) 为模型M 中的自由参数的数目，在线性回归中为特征数目。

首先读取数据：

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

标准化：

if doStandardize
  [Xtrain, mu]  = center(Xtrain);
  [Xtrain, s]  = mkUnitVariance(Xtrain);
  if ~isempty(Xtest)
    Xtest = center(Xtest, mu);
    Xtest = mkUnitVariance(Xtest, s);
  end
end

然后似然估计噪声方差

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

特征组合并穷举记录结果：

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))
    continue
  else
    %fprintf('including %s\n', sprintf('%d,', include));
  end
  
  members{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 mseTest
  mseTest_BIC(j) = mseTrain(j) + log(n)*sz(j)*sig2_hat/n;   % BIC; 
  
  %test
  if ~isempty(Xtest)
    ntest = size(Xtest, 1);
    ypredTest = [ones(ntest,1) Xtest(:,include)]*w{j};
    mseTest(j) = mean((ypredTest-ytest).^2);
  else
    mseTest = [];
  end
  
  j = j + 1;
end

选出最小误差结果，并画出来：

for i=0:d
  ndx = 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 的结果
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
bestmseTest

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

结果图：

其中参数结果为：

该模型的训练误差为0.5210，在测试集上的真实测试误差为0.4815。

方法三：交叉验证

如果用10折交叉验证测试误差估计，并采用1倍标准差原则的话，得到的最佳模型为子集大小为k=2。

首先读取数据：

load('prostate.mat') % from prostateDataMake

[n d] = size(Xtrain);
[w, mseTrain, mseTest,model] =  allSubsetsRegression_CV(Xtrain, ytrain, Xtest, ytest, 0:d, 1);

标准化数据：

%标准化数据
if doStandardize
  [Xtrain, mu]  = center(Xtrain);
  [Xtrain, s]  = mkUnitVariance(Xtrain);
  if ~isempty(Xtest)
    Xtest = center(Xtest, mu);
    Xtest = mkUnitVariance(Xtest, s);
  end
end

将数据分成10折，比如100 个数据，第一次[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 交叉验证

分别选每折数据进行，组合穷举：

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

算出每个折的一倍标准差，以及每折得到的对应的个数求出的均方差的平均。

%mseCVTest  这个得到的是那个最好的测试的结果，就是每折之后， 选出的[], 1，2 ，3 ...,8 枚举过的最小的均方误差。
mseCVMean = mean(mseCVTest,1);  %每一折的均方差的平均 ， 即列的平均啊， 最小的那个测试标准偏差。
mseCVse = std(mseCVTest,[],1)/sqrt(Nfolds);   %一倍标准差， 标准偏差， 1×9

然后是其一倍标准差表示，以及找到在变量最少（模型最简单）的情况下，满足最小均方差一倍标准差范围的。

figure;                    %误差带的上下各为 mseCVse
errorbar(df, mseCVMean, mseCVse);   % confidence intervals of data or the deviation along a curve.
kstar = oneStdErrorRule(mseCVMean, mseCVse);  %找出根据均值和一倍标准差， 最好的一个满足的。， 因为最小到大其特征数在增加，因而其复杂度在增加。
hold on
vline(kstar,'r-.');
xlabel('size of subset')
ylabel('cv error')
title(sprintf('all subsets CV'))

通过得到这个最少变量得个数，在这个最小变量个数下能表示此模型的情况。

% Refit using all training data to find set of chosen size
[w, mseTrain, junk, junk2, members] = ...
    allSubsetsRegression(Xtrain, ytrain, [], [], sizes(kstar), 0);  %所以找到两个的最佳的那个w， 
[junk,best]=min(mseTrain);
w = w(best);
model = members{best};
mseTrain = mseTrain(best);

model

if ~isempty(Xtest)
    ntest = size(Xtest, 1);
    ypredTest = [ones(ntest,1) Xtest(:,model)]* cell2mat(w);
    mseTest = mean((ypredTest-ytest).^2);
  else
    mseTest = 0;
end
  
mseTrain
mseTest

此时模型的参数结果为：

该模型的训练误差为0. 5536，在测试集上的真实测试误差为0.5737。

方法四：正交标准寻踪

我们以BIC准则估计模型的测试误差，采用正交投影寻踪每次增加一个与当前残差最相关的特征，然后更新模型参数。前向逐步回归选择的模型为子集大小为4；当特征数目变成5时，BIC分数增加，迭代终止。

读取数据：

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

数据标准化：

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;

当数据为空时，最开始的初始化情况有：

% 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);
else
    mseTest = [];
end

ndx_star = j;

然后进入OMP 逐渐的迭代过程，利用BIC 分数作为判别的标准，当BIC 分数增加时（误差变大）迭代终止：

j = j+1;
for i = 1 : d
  %当前模型的残差
  residual = ytrain-ypredTrain;
  %残差与特征之间的相关性（投影）
  proj = Xtrain' * residual;
  %选择与残差最相关的特征
  [maxval, bestk] = max(abs(proj));
  
  include = [include bestk];        %add the mink feature
  members{j} = include;        % add best k th feature temp
  sz(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 mseTest
  mseTest_BIC(j) = mseTrain(j) + log(n)*sz(j)*sig2_hat/n;   % BIC;  利用此分数来进行终止条件判断。
      
  %test
  if ~isempty(Xtest)
    ntest = size(Xtest, 1);
    ypredTest = [ones(ntest,1) Xtest(:,members{j})]* w{j} ;
    mseTest(j) = mean((ypredTest-ytest).^2);
  else
    mseTest = [];
  end

  
  if mseTest_BIC(j) > bestBIC     % the BIC tend to increase
      break;
  else
      bestBIC = mseTest_BIC(j);
  end
  
  j = j + 1;
end

结果：