Lasso、Ridge Regression(岭回归)与线性回归模型(含Matlab代码)

引言

LASSO是由1996年Robert Tibshirani首次提出，全称Least absolute shrinkage and selection operator。该方法是一种压缩估计。它通过构造一个惩罚函数得到一个较为精炼的模型，使得它压缩一些回归系数，即强制系数绝对值之和小于某个固定值；同时设定一些回归系数为零。因此保留了子集收缩的优点，是一种处理具有复共线性数据的有偏估计。¹
本文立足于普通线性回归模型，首先提及普通最小二乘法(OLS)，并说明Lasso与Ridge Regression(岭回归)为OLS的一个延申方法，用于处理OLS无法解决的问题。
本文假设读者已有概率论与数理统计、线性代数以及运筹学的先验知识，并对线性回归模型有着一定的了解。

一、从线性回归说起

线性回归模型是统计学中最基础也是最实用的模型，他旨在给出一个形如
$$y={\beta }_0+{\mathbf{\beta }}^T\mathbf{x}$$
的线性表达式，以表示出自变量$y$与因变量$\mathbf{x}$之间的线性关系
(其中$\mathbf{\beta }$与$\mathbf{x}$都是一个$d$维的列向量)

1.1 模型假设

本文始终按照如下假设，后文不再进行叙述：
现 有 n 组 数 据 样 本 点 { ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( x n , y n ) } 现有n组数据样本点\{(\boldsymbol{x}_1,y_1),(\boldsymbol{x}_2,y_2),...,(\boldsymbol{x}_n,y_n)\} 现有n组数据样本点{(x1​,y1​),(x2​,y2​),...,(xn​,yn​)}
$$\begin{gathered} 其中{\mathbf{x}}_i\in R^d是一个d维的列向量，y_i\in R是一个实数 \\ 也就是说现在有d个自变量与1个因变量，这种情况被称为MutipleLinearRegression(多元线性回归) \\ 需要警惕不能与MultivariateLinearRegression(因变量y亦为向量)混淆 \end{gathered}$$

1.2 普通最小二乘法

1.2.1 历史由来

1801年，意大利天文学家朱赛普·皮亚齐发现了第一颗小行星谷神星。经过40天的跟踪观测后，由于谷神星运行至太阳背后，使得皮亚齐失去了谷神星的位置。随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星，但是根据大多数人计算的结果来寻找谷神星都没有结果。时年24岁的高斯也计算了谷神星的轨道。奥地利天文学家海因里希·奥尔伯斯根据高斯计算出来的轨道重新发现了谷神星。
高斯使用的最小二乘法的方法发表于1809年他的著作《天体运动论》中。
法国科学家勒让德于1806年独立发明“最小二乘法”，但因不为世人所知而默默无闻。勒让德也曾与高斯为谁最早创立最小二乘法原理发生争执。
1829年，高斯提供了最小二乘法的优化效果强于其他方法的证明，因此被称为高斯-马尔可夫定理(OLS估计必定为BLUE)。²

1.2.2 最小二乘法的理论依据

相信大家对最小二乘法已经是非常熟悉了，其本质是一个优化模型，找到最优的回归系数${\beta }_0$与$\mathbf{\beta }$使得以下这个式子取到最小值
$$SSE=\sum _{i=1}^n{\left(y_i-{\beta }_0-{\mathbf{\beta }}^T{\mathbf{x}}_i\right)}^2$$
上面这个式子称为残差平方和，它在标准误差的估计以及线性回归的检验中还有妙用，不过本文仅仅考虑最小二乘法本身而不深究后续的统计推断，故仅需要知道最小二乘法的本质是找到使残差平方和最小的回归系数。
由于样本量$n$也会对残差平方和的大小产生影响，我们不妨以平均残差平方和为探究目标，即以
$$Q({\beta }_0,\mathbf{\beta })=\frac{1}{n}\sum _{i=1}^n{\left(y_i-{\beta }_0-{\mathbf{\beta }}^T{\mathbf{x}}_i\right)}^2$$
为目标函数，故最小二乘法的理论模型可以写为：
$$({\beta }_0^{\ast },{\mathbf{\beta }}^{\ast })=argmi{n}_{({\beta }_0,\mathbf{\beta })}Q({\beta }_0,\mathbf{\beta })$$

1.2.3 模型的化简

实际上，我们可以假设${\mathbf{x}}_i$与$y_i$的均值皆为0，即
$$\sum _{i=1}^n{\mathbf{x}}_i=0,\sum _{i=1}^n{y}_i=0$$
以便化简模型，下面我们来证明这个假设是可行的
$证明$
$$\begin{gathered} 由于原模型是一个无约束优化问题，且其目标函数是个处处可导函数，取到极小值必定有各偏导为0，那我们不妨尝 \\ 试对常数项系数{\beta }_0进行求导并令其为0 \end{gathered}$$
$$令\frac{\partial Q({\beta }_0,\mathbf{\beta })}{\partial {\beta }_{)}}=-\frac{2}{n}\sum _{i=1}^n\left(y_i-{\beta }_0-{\mathbf{\beta }}^T{\mathbf{x}}_i\right)=0$$
$$\begin{gathered} 解得{\beta }_0=\bar{y}-{\mathbf{\beta }}^T\bar{\mathbf{x}} \\ 将其代回Q({\beta }_0,\mathbf{\beta })，可得 \\ Q(\mathbf{\beta })=\frac{1}{n}\sum _{i=1}^n{\left((y_i-\bar{y})-{\mathbf{\beta }}^T({\mathbf{x}}_i-\bar{\mathbf{x}})\right)}^2 \\ 故不妨假设\sum _{i=1}^n{\mathbf{x}}_i=0,\sum _{i=1}^n{y}_i=0，否则取{\mathbf{x}}_i^{\ast }={\mathbf{x}}_i-\bar{\mathbf{x}},y_i^{\ast }=y_i-\bar{y} \\ 此时仅需研究一个不含常数项的函数Q(\mathbf{\beta })在何时取到最小值，再由回代法求出{\beta }_0^{\ast } \end{gathered}$$
由上述证明，下文仅仅讨论不含常数项${\beta }_0$的模型，即：
$${\mathbf{\beta }}^{\ast }=argmi{n}_{\mathbf{\beta }}\frac{1}{n}\sum _{i=1}^n{\left(y_i-{\mathbf{\beta }}^T{\mathbf{x}}_i\right)}^2$$
如果记$\mathbf{y}=[y_1,y_2,...,y_n{]}^T\in R^n,\mathbf{X}=[{\mathbf{x}}_1,{\mathbf{x}}_2,...,{\mathbf{x}}_n{]}^T\in R^{n\times d}$，由矩阵知识可知
$$\sum _{i=1}^n{\left(y_i-{\mathbf{\beta }}^T{\mathbf{x}}_i\right)}^2=\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2$$
其中$\Vert \mathbf{\alpha }{\Vert }_2$被称作向量的L2-范数，不明其意的读者可在下一章节看见其具体含义。
故线性规划模型可简化为：
$${\mathbf{\beta }}^{\ast }=argmi{n}_{\mathbf{\beta }}\frac{1}{n}\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2$$

二、Ridge Regression(岭回归)与Lasso

一般来说，回归问题是一个函数拟合的过程，那么我们希望模型不要太复杂，否则很容易发生过拟合现象，所以我们要加入正则化项，而不同的正则化项就产生了不同的回归方法，其中以Ridge Regression (岭回归）和Lasso最为经典，前者是加入了L2正则化项，后者加入的是L1正则化项。

2.1 向量的范数

本文先对向量的L1范数和L2范数进行介绍，以便后面的介绍。
$\Vert \mathbf{\alpha }{\Vert }_1$被称为向量的L1范数，$\Vert \mathbf{\alpha }{\Vert }_2$被称为向量的L2范数
$$\begin{gathered} \Vert \mathbf{\alpha }{\Vert }_1=\sum _{i=1}^d\left|a_i\right| \\ \Vert \mathbf{\alpha }{\Vert }_2=\sqrt{\sum _{i=1}^d{a}_i^2} \end{gathered}$$

2.2 Ridge Regression与Lasso的模型

Ridge Regression与Lasso都是普通线性回归本质上也是优化模型，而前文提到过，相比于普通线性回归，Ridge Regression加入了L2正则化项，Lasso加入了L1正则化项，故它们的目标函数如下：
$$\begin{gathered} Lasso:{\mathbf{\beta }}^{\ast }=argmi{n}_{\mathbf{\beta }}\frac{1}{n}\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2+\lambda \Vert \mathbf{\beta }{\Vert }_1 \\ RidgeRegression:{\mathbf{\beta }}^{\ast }=argmi{n}_{\mathbf{\beta }}\frac{1}{n}\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2+\lambda \Vert \mathbf{\beta }{\Vert }_2^2 \end{gathered}$$

2.3 Ridge Regression与Lasso的差异

Ridge Regression与Lasso作为普通线性回归的改进方法，二者的差异仅仅在于加入的正则化项不同，那么二者的差异在何处呢？它们分别能起到什么样的作用？

2.3.1 模型转化

岭回归模型有一个等价的优化问题，即考虑一个带有约束的线性回归问题，其中约束条件是将系数向量$\mathbf{\beta }$限制在一个半径为$t$的超球内部，即下面这个模型是与原问题等价的：
$$\begin{gathered} {\mathbf{\beta }}^{\ast }=argmi{n}_{\mathbf{\beta }}\frac{1}{n}\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2 \\ s.t.\Vert \mathbf{\beta }{\Vert }_2^2\le t \end{gathered}$$
$证明$
$$\begin{gathered} 对上面带约束的优化模型使用广义的拉格朗日乘数法，将其转化为无约束优化问题，便有 \\ L(\mathbf{\beta },\lambda )=\frac{1}{n}\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2+\lambda \Vert \mathbf{\beta }{\Vert }_2^2-\lambda t \\ 当\lambda 与t给定时，可以看出L(\mathbf{\beta },\lambda )是RidgeRegression的目标函数(\lambda t是一个常数) \\ 故我们证明了上面的优化模型可推得原问题，又由于所求优化问题是一个凸优化(线 \\ 性回归是一个典型的凸优化问题)，具有强对偶性，故两个模型等价 \end{gathered}$$
同理地，Lasso也有一个带约束形式的等价优化模型：
$$\begin{gathered} {\mathbf{\beta }}^{\ast }=argmi{n}_{\mathbf{\beta }}\frac{1}{n}\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2 \\ s.t.\Vert \mathbf{\beta }{\Vert }_1\le t \end{gathered}$$

2.3.2 二者差异比较

考虑$d=2$的情况，记$\mathbf{\beta }=[{\beta }_1,{\beta }_2{]}^T$，给出一组样本数据：

容易验证这组数据是满足均值为0条件的，即
$$\sum _{i=1}^n{\mathbf{x}}_i=0,\sum _{i=1}^n{y}_i=0$$
故我们研究不含常数项的回归模型
$${\mathbf{\beta }}^{\ast }=argmi{n}_{\mathbf{\beta }}Q(\mathbf{\beta })=\frac{1}{n}\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2$$
由于此时$\mathbf{\beta }=[{\beta }_1,{\beta }_2{]}^T$，则$Q$是${\beta }_1,{\beta }_2$的函数，不妨作出其函数图像

可以看出该函数有且仅有一个极小值，为更加直观的观察，我们绘制出其等值线图

等值线中，红色点是$Q$的极小值点，同一颜色的圈代表着$Q$在该圈上函数值相同，越靠近红色点则$Q$的值越小，很显然对于这个问题而已普通线性回归的解就是红色点。
下面我们来看看Ridge Regression与Lasso的解，从它们的带约束条件的优化模型可以看出，在$d=2$场合下，Ridge Regression将可行解约束在一个圆内，而Lasso将可行解约束在一个菱形内，具体情况可见下图：

从上图可以看见，Ridge Regression的最优解在某条等值线与约束圆的切点上取到，而Lasso的最优解通常在约束菱形的端点处取到（即坐标轴上），也就是说Lasso可能会使某个或某些回归系数为0，也就是说达到了“属性约简”或是“特征提取”的效果。而Ridge Regression会导致“权值衰减”。

三、求解方法

从目标函数的表达式就可以看出，普通线性回归与Ridge Regression的目标函数是处处可导的，而Lasso目标函数就没有那么好的性质了。实际上，我们可以给出前两个模型的解析解，而Lasso仅能给出迭代解法。

3.1 矩阵解法

3.1.1 矩阵的迹及其导数

矩阵的迹大家应该非常了解，即对角线元素相加（要求是方阵）：
$$tr(A)=\sum _{i=1}^n{a}_{ii}$$
下面介绍矩阵迹的一些性质，仅仅叙述，不加以证明
性质1：$tr(AB)=tr(BA)$
性质2：$\frac{\partial tr(AB)}{\partial (A)}=\frac{\partial tr(BA)}{\partial (A)}=B^T$
性质3：$\frac{\partial tr(A^{T}B)}{\partial (A)}=\frac{\partial tr(B{A}^T)}{\partial (A)}=B$
性质4：$\frac{\partial tr(A^{T}B)}{\partial (A)}=\frac{\partial tr(B{A}^T)}{\partial (A)}=B$
性质5：$\frac{\partial tr(A^{T}A)}{\partial (A)}=\frac{\partial tr(A{A}^T)}{\partial (A)}=2A$
性质6：$\frac{\partial tr(A^{T}BA)}{\partial (A)}=(B+B^T)A$
性质7：$\frac{\partial tr(AB{A}^T)}{\partial (A)}=A(B+B^T)$
性质8：$\frac{\partial tr(BA{A}^T{B}^T)}{\partial (A)}=(A^T{B}^{T}B{)}^T+B^{T}BA=2B^{T}BA$
向量的L2-范数与矩阵的迹有着密不可分的关系，容易证明下式：
$$\begin{gathered} 若\mathbf{\alpha }是d维列向量，则有 \\ \Vert \mathbf{\alpha }{\Vert }_2^2={\mathbf{\alpha }}^T\mathbf{\alpha }=tr(\mathbf{\alpha }{\mathbf{\alpha }}^T) \end{gathered}$$

3.1.2 最小二乘法的矩阵解

由向量L2-范数与矩阵迹的关系，可以将最小二乘法的目标函数改写为：
$$\begin{gathered} Q(\mathbf{\beta })=\frac{1}{n}\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2=\frac{1}{n}tr[(\mathbf{y}-\mathbf{X}\mathbf{\beta })(\mathbf{y}-\mathbf{X}\mathbf{\beta }{)}^T] \\ =\frac{1}{n}tr[(\mathbf{y}-\mathbf{X}\mathbf{\beta })({\mathbf{y}}^T-{\mathbf{\beta }}^T{\mathbf{X}}^T)] \\ =\frac{1}{n}tr[\mathbf{y}{\mathbf{y}}^T-\mathbf{y}{\mathbf{\beta }}^T{\mathbf{X}}^T-\mathbf{X}\mathbf{\beta }{\mathbf{y}}^T+\mathbf{X}\mathbf{\beta }{\mathbf{\beta }}^T{\mathbf{X}}^T] \end{gathered}$$
$$\begin{gathered} 再令\frac{\partial Q(\mathbf{\beta })}{\partial \mathbf{\beta }}=\frac{1}{n}\frac{\partial tr[\mathbf{y}{\mathbf{y}}^T-\mathbf{y}{\mathbf{\beta }}^T{\mathbf{X}}^T-\mathbf{X}\mathbf{\beta }{\mathbf{y}}^T+\mathbf{X}\mathbf{\beta }{\mathbf{\beta }}^T{\mathbf{X}}^T]}{\partial \mathbf{\beta }} \\ =\frac{1}{n}[0-{\mathbf{X}}^T\mathbf{y}-{\mathbf{X}}^T\mathbf{y}+2{\mathbf{X}}^T\mathbf{X}\mathbf{\beta }] \\ =\frac{2}{n}({\mathbf{X}}^T\mathbf{X}\mathbf{\beta }-{\mathbf{X}}^T\mathbf{y})=0 \end{gathered}$$
$$得{\mathbf{\beta }}^{OLS}=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}$$

3.1.3 Ridge Regression的矩阵解

从最小二乘法解的形式可以看出，要求${\mathbf{X}}^T\mathbf{X}$是一个非奇异阵，即要求影响变量之间不存在共线性，一旦变量之间出现了共线性，最小二乘法就难以实施，Ridge Regression就此而生。
同样地，将Ridge Regression的目标函数写为
$$\begin{gathered} R(\mathbf{\beta })=\frac{1}{n}\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2+\lambda \Vert \mathbf{\beta }{\Vert }_2^2 \\ =Q(\mathbf{\beta })+\lambda tr(\mathbf{\beta }{\mathbf{\beta }}^T) \end{gathered}$$
$$\begin{gathered} 令\frac{\partial R(\mathbf{\beta })}{\partial \mathbf{\beta }}=\frac{\partial Q(\mathbf{\beta })}{\partial \mathbf{\beta }}+\lambda \frac{\partial tr(\mathbf{\beta }{\mathbf{\beta }}^T)}{\partial \mathbf{\beta }} \\ =\frac{2}{n}({\mathbf{X}}^T\mathbf{X}\mathbf{\beta }-{\mathbf{X}}^T\mathbf{y})+2\lambda \mathbf{\beta }=0 \end{gathered}$$
$$得{\mathbf{\beta }}^{Ridge}=({\mathbf{X}}^T\mathbf{X}+\lambda n\mathbf{I}{)}^{-1}{\mathbf{X}}^T\mathbf{y}$$
可以看见，通过适当的选择$\lambda$的值，可以使得原本是奇异阵的${\mathbf{X}}^T\mathbf{X}$变为一个可逆矩阵${\mathbf{X}}^T\mathbf{X}+\lambda n\mathbf{I}$，而$\lambda$的值可以人为设定，称为岭回归的参数(也有时称$\alpha =\lambda n$为岭回归参数)，从模型中可以看出，$\lambda$越大则对正则化约束越为看重，当$\lambda$较小时，可以得到一个与线性回归几乎一样的解。

3.2 Lasso的求解方法

3.2.1 解析解法

Lasso的目标函数不是一个处处可导的函数，且对向量的L1-范数求导需要引入次梯度的概念，且还需要假设${\mathbf{X}}^T\mathbf{X}=\mathbf{I}$，也就是任意两个特征向量相互正交，才能导出Lasso的解析解，由于条件过于苛刻，在实际运用中几乎不可能用到解析解，在此仅仅给出结论，省去推导过程。
$${\beta }_i^{Lasso}=sign({\beta }_i^{OLS})(|{\beta }_i^{OLS}|-\frac{\lambda n}{2}{)}_{+}$$
注：其中$sign$表示符号函数，$(x{)}_{+}$表示$x$的正部

3.2.2 数值解法

坐标下降法是用来求解Lasso的一个经典方法，它适用于在1维场合容易解决但高维场合难以解决的问题，下面对算法进行叙述
$Step1:随机化确定{\mathbf{\beta }}^{(0)}=({\beta }_1^{(0)},{\beta }_2^{(0)},...,{\beta }_d^{(0)})$
$$\begin{gathered} Step2:假设现在有{\mathbf{\beta }}^{(k)}=({\beta }_1^{(k)},{\beta }_2^{(k)},...,{\beta }_d^{(k)})，则根据如下公式进行迭代： \\ {\beta }_i^{(k+1)}=argmi{n}_{{\beta }_i}L({\beta }_1^{(k)},{\beta }_2^{(k)},...,{\beta }_{i-1}^{(k)},{\beta }_i,{\beta }_{i+1}^{(k)},...,{\beta }_d^{(k)}),这里L是Lasso的目标函数 \\ (也就是说此时只有{\beta }_i是变量，其他都为常量，是一个一维搜索问题，可采用一维搜索方法求解或直接使用3.2.1中 \\ 提到的解析解) \end{gathered}$$
$Step3:对\forall i,若|{\beta }_i^{(k+1)}-{\beta }_i^{(k)}|<ϵ，则跳出循环，否则回到Step2$

四、程序实现

下面介绍一下如何用Matlab软件实现上述模型。先给出一个背景，我们来研究这样一个问题：

例1 ：财政收入与国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资等因素有关。下表列出了1952-1981年的原始数据，试构造预测模型。

4.1 线性回归

只要你的Matlab版本不低于2012，那么它就帮你集成了lasso这个函数(线性回归与Ridge Regression更是在更老的版本就有了)，先从最简单的线性回归说起。
线性回归在Matlab中的命令为regress，其基本用法为：
b = regress(y,X)
其中y为因变量组成的列向量，X为自变量组成的矩阵(还有更多用法，读者请自行探索)

%% 读取数据
data = [1952	598	349	461	57482	20729	44	184
    1953	586	455	475	58796	21364	89	216
    1954	707	520	491	60266	21832	97	248
    1955	737	558	529	61465	22328	98	254
    1956	825	715	556	62828	23018	150	268
    1957	837	798	575	64653	23711	139	286
    1958	1028	1235	598	65994	26600	256	357
    1959	1114	1681	509	67207	26173	338	444
    1960	1079	1870	444	66207	25880	380	506
    1961	757	1156	434	65859	25590	138	271
    1962	677	964	461	67295	25110	66	230
    1963	779	1046	514	69172	26640	85	266
    1964	943	1250	584	70499	27736	129	323
    1965	1152	1581	632	72538	28670	175	393
    1966	1322	1911	687	74542	29805	212	466
    1967	1249	1647	697	76368	30814	156	352
    1968	1187	1565	680	78534	31915	127	303
    1969	1372	2101	688	80671	33225	207	447
    1970	1638	2747	767	82992	34432	312	564
    1971	1780	3156	790	85229	35620	355	638
    1972	1833	3365	789	87177	35854	354	658
    1973	1978	3684	855	89211	36652	374	691
    1974	1993	3696	891	90859	37369	393	655
    1975	2121	4254	932	92421	38168	462	692
    1976	2052	4309	955	93717	38834	443	657
    1977	2189	4925	971	94974	39377	454	723
    1978	2475	5590	1058	96259	39856	550	922
    1979	2702	6065	1150	97542	40581	564	890
    1980	2791	6592	1194	98705	41896	568	826
    1981	2927	6862	1273	100072	73280	496	810];

%% 数据处理
x = data(1:end-5,2:end-1);     % 留5个用于预测
y = data(1:end-5,end);  % 留5个用于预测
X = [ones(size(y)) x];  % 加上常数项

b_OLS = regress(y,X) %  线性回归

运行结果为：

【output】
b_OLS =

  365.5804
    0.5610
    0.0137
   -0.5819
   -0.0021
   -0.0050
    0.0713

即${\beta }_0=365.5804,{\mathbf{\beta }}_1=(0.5610,0.0137,-0.5819,-0.0021,-0.0050,0.0713{)}^T$

4.2 Ridge Regression

b = ridge(y,X,k,scaled)
这是岭回归的命令，其中y表示因变量组成的向量，X表示自变量组成的向量(无需带有常数项)，k为岭回归参数，scaled表示是否输出标准化系数(0为否，1为是)

b_Ridge = ridge(y,x,1,0) % Ridge

首先选用岭回归参数=1进行岭回归，可以发现结果与线性回归相去甚远

【output】
b_Ridge =

   47.0837
    0.1266
    0.0305
   -0.0705
    0.0006
    0.0017
    0.5002

若选用岭回归参数=0.001，输出结果为：

b_Ridge =

  363.9392
    0.5577
    0.0140
   -0.5781
   -0.0021
   -0.0049
    0.0744

这个结果与线性回归的结果非常接近了！

最近时间比较忙, 待续…

---

1. 引用自“LASSO——百度百科” ↩︎

2. 引用自“wikipedia” ↩︎