缩减系数

一、岭回归
1、原理
如果对代价函数$J(\theta )=(X\theta -y{)}^T(X\theta -y)$不加任何限制，那么对$\theta$求导，可以得到$\theta$矩阵形式的解$\theta =(X^{T}X{)}^{-1}{X}^{T}y$。这里会出现一个问题：如果得到的样本个数小于样本的特征数，也就是说训练集矩阵$X$的行维度$m$小于列维度$d$,那么$X^{T}X$可能不可逆，也就是无法解出$\theta$。为了避免这个问题，我们给$X^{T}X$加上一个$\lambda I$，最终的解也变为$\theta =(X^{T}X+\lambda I{)}^{-1}{X}^{T}y$。
实质上，岭回归是一种缩减系数的方法，也就是剔除或者减弱不重要的参数：我们知道$L_2$正则化下的代价函数为：$J(\theta )=\frac{1}{2m}{\sum }_{i=1}^m(h_{\theta }(x_i)-y_i)+\frac{\lambda }{2m}{\sum }_{i=1}^m{\theta }_i^2$,如果写成矩阵形式则为：$J(\theta )=\frac{1}{2m}(X\theta -y{)}^T(X\theta -y)+\frac{\lambda }{2m}{\theta }^T\theta$，对$J(\theta )$关于$\theta$求导并令其为0：
$$\frac{\partial J(\theta )}{\partial \theta }=\frac{1}{m}{X}^T(X\theta -y)+\frac{\lambda }{m}\theta =0$$$$\to$$$$\theta =(X^{T}X+\lambda I)X^{T}y$$
也就是说岭回归等同于就是$L_2$正则化，而对于$L_2$正则化的问题$\min J(\theta )\to \min ({\sum }_{i=1}^m(h_{\theta }(x_i{)}^2-y_i)+\lambda {\sum }_{i=1}^m{\theta }_i^2)$就是约束问题：
$$\min \sum _{i=1}^m(h_{\theta }(x_i)-y_i{)}^2$$$$s.t.$$$$\sum _{i=1}^m{\theta }_i^2\le t$$
也就是说在满足条件${\sum }_{i=1}^m{\theta }_i^2\le t$的前提下最小化平方误差。而约束项${\sum }_{i=1}^m{\theta }_i^2\le t$会导致系数${\theta }_i$被压缩，甚至一些系数趋于0，这就达成了缩减系数的目的。
实际上，缩减系数和正则化也并不矛盾：正则化的目的为防止过拟合，而过拟合的原因一部分是由于样本太少，而特征太多导致模型会把一些不重要的数据特征当做所有数据的共性。而缩减系数正是排除或减弱了不重要特征。
2、实际应用
这里我们我们把岭回归应用于一组数据，看看其缩减系数的效果
①查看并处理数据
首先查看数据：

file = open('F:/MachineLearning/data/abalone.txt')
print(file.read())

得到：

F,0.67,0.51,0.155,1.278,0.5605,0.3045,0.358,11
M,0.67,0.54,0.195,1.217,0.532,0.2735,0.3315,11
F,0.67,0.54,0.2,1.46,0.6435,0.328,0.4165,9
F,0.675,0.535,0.185,1.5575,0.7035,0.402,0.4,11
M,0.675,0.51,0.17,1.527,0.809,0.318,0.341,11
F,0.675,0.53,0.195,1.4985,0.62,0.375,0.425,9
M,0.685,0.55,0.19,1.885,0.89,0.41,0.4895,10
M,0.685,0.535,0.175,1.432,0.637,0.247,0.46,11
M,0.705,0.55,0.21,1.4385,0.655,0.3255,0.462,11
F,0.705,0.53,0.17,1.564,0.612,0.394,0.44,10
M,0.71,0.555,0.175,2.14,1.2455,0.3725,0.434,11
F,0.725,0.56,0.185,1.792,0.873,0.367,0.435,11
M,0.78,0.6,0.21,2.548,1.1945,0.5745,0.6745,11
I,0.235,0.13,0.075,0.1585,0.0685,0.037,0.0465,5
I,0.35,0.25,0.1,0.4015,0.1725,0.063,0.1255,7
I,0.36,0.25,0.115,0.465,0.21,0.1055,0.128,7
I,0.38,0.28,0.095,0.2885,0.165,0.0435,0.067,7
F,0.38,0.32,0.115,0.6475,0.323,0.1325,0.164,7
M,0.43,0.31,0.13,0.6485,0.2735,0.163,0.184,9
I,0.465,0.36,0.105,0.452,0.22,0.159,0.1035,9
I,0.47,0.355,0.12,0.4915,0.1765,0.1125,0.1325,9
F,0.485,0.365,0.15,0.9145,0.4145,0.199,0.273,7
M,0.495,0.375,0.155,0.976,0.45,0.2285,0.2475,9
I,0.5,0.395,0.145,0.7865,0.332,0.1815,0.2455,8
M,0.505,0.4,0.15,0.775,0.3445,0.157,0.185,7
I,0.51,0.375,0.15,0.8415,0.3845,0.156,0.255,10
M,0.51,0.38,0.135,0.681,0.3435,0.142,0.17,9
M,0.515,0.37,0.115,0.6145,0.3415,0.155,0.146,9
F,0.55,0.415,0.18,1.1655,0.502,0.301,0.311,9
F,0.575,0.42,0.19,1.764,0.914,0.377,0.4095,10
M,0.605,0.455,0.16,1.1215,0.533,0.273,0.271,10
M,0.615,0.505,0.165,1.167,0.4895,0.2955,0.345,10
M,0.615,0.475,0.15,1.0375,0.476,0.2325,0.283,9
...

这里的数据是一组鲍鱼一系列长度数据，如体长、厚度等；我们知道体长厚度一定是和其年龄有一定关系，即年龄越大，体长重量等数据越大，即成正相关，所以这里我们可以采用线性回归粗略的估计。这组数据需要注意的是：1>第一个特征为鲍鱼的性别，M代表母，F代表公，I代表幼儿，我们这里分别用0，1，-1代替；2>最后一个特征为鲍鱼年龄即我们需要预测的数据，所以应该单独提出；3>这里的数据中不同样本的一个属性的数据差值略大，所以应该进行归一化处理。所以处理数据的函数如下：

import numpy as np

def prepared_data(path):
    file=open(path)
    fr_list = file.readlines()
    data_list=[]
    for data in fr_list:
        data = data.strip()
        data = data.split(',')
        if data[0] == 'F':
            data[0] = 1
        elif data[0] == 'I':
            data[0] = 0
        else:
            data[0] =-1
        data = [float(i) for i in data]
        data_list.append(data)
        data_mat = np.array(data_list)
        X = data_mat[:,:-1]
        y = data_mat[:,-1]
        X = (X-np.mean(X,axis=0)) / np.std(X,axis=0)
        y = y - np.mean(y)
    return X,y

处理并查看数据：

X,y = prepared_data('F:/MachineLearning/data/abalone.txt')
print(X)
print(y)

得到：

[[-1.15198011 -0.57455813 -0.43214879 ... -0.60768536 -0.72621157
  -0.63821689]
 [-1.15198011 -1.44898585 -1.439929   ... -1.17090984 -1.20522124
  -1.21298732]
 [ 1.28068972  0.05003309  0.12213032 ... -0.4634999  -0.35668983
  -0.20713907]
 ...
 [-1.15198011  0.6329849   0.67640943 ...  0.74855917  0.97541324
   0.49695471]
 [ 1.28068972  0.84118198  0.77718745 ...  0.77334105  0.73362741
   0.41073914]
 [-1.15198011  1.54905203  1.48263359 ...  2.64099341  1.78744868
   1.84048058]]
   [ 5.06631554 -2.93368446 -0.93368446 ... -0.93368446  0.06631554
  2.06631554]

③岭回归代码
岭回归代码较为简单，这里注意的是$X^{T}X+\lambda I$不一定可逆，所以应判断：

def ridge_regression(X,y,lambd):
    m,d = X.shape
    X_X = X.T @ X + lambd * np.eye(d)
    if np.linalg.det(X_X) == 0:
        print('The mat can not be inverse')
        return 
    theta = np.linalg.inv(X_X) @ X.T @ y
    return theta

file = open('F:/MachineLearning/data/abalone.txt')
print(file.read())

④运行并查看缩减结果

theta = ridge_regression(X,y,100)
print(theta)

theta_inx=np.argsort(-np.abs(theta))
print(theta_inx)
print(theta[theta_inx])
fig1,ax1 = plt.subplots(figsize=(9,6))
ax1.plot(np.abs(theta[theta_inx]))
plt.show()

这里，我们设置$\lambda$为100，并查看相同$\lambda$下不同属性对应参数的绝对值：

可以看到最后一个属性对应的参数几乎趋近于0。如果设置大一些，如1000：

发现几乎所有参数绝对值都小于1，一定会于发生偏差，产生欠拟合。我们再把参数设置小一些，如0.01：

也起到了一定效果，但是效果并没有前面明显。或者说这是应为模型计算出的几个参数本身就很小。
从上面可以看出$\lambda$的选择对于岭回归模型即为重要，具体选择还需要选取一定数量的$\lambda$在测试集中验证，选择在测试集中得到最小错误率的$\lambda$。

二、LASSO
上面我们看到了利用$L_2$正则化变换来的岭回归的缩减系数的方法，$L_2$正则化为：
$$\min \sum _{i=1}^m(h_{\theta }(x_i)-y_i)+\sum _{i=1}^m{\theta }_i^2$$
如果我们把$L_2$正则化换为$L_1$正则化，那么问题又会怎样呢？这是问题就变为:
$$\min \sum _{i=1}^m(h_{\theta }(x_i)-y_i)+\sum _{i=1}^m|{\theta }_i|$$
(这里省去了与极小化无关的$\frac{1}{2m}$)
如果我们把两种情况用图来表示，则会有如下的图：

这里为了更直观的说明LASSO和岭回归的区别，我们假设模型只含有两个参数${\theta }_1$和${\theta }_2$:我们先画出平方误差项等值线，即在$({\theta }_1,{\theta }_2)$空间中平方误差项取值相同的点的连线，再画出$L_1$、$L_2$正则化的等值线。最小值一定出现在平方误差等值线和正则化项的等值线的相交处。而如果采用$L_2$正则化，相交处大多位于非坐标轴上；而采用$L_1$正则化，相交处大多位于坐标上。
所以总结一下，就是说岭回归得到的参数虽然被压缩，但是几乎不可能完全把参数压缩为0；而LASSO易于得到稀疏解，也就是$\theta$存在非0分量。由于LASSO这个特性，也可用来特征选择。
而如何得到LASSO的解呢？由于绝对值函数$|\theta |$的不可导，所以希望得到像岭回归那样的显式解是很困难的，所以采用近端梯度下降的方法。对于优化目标$\min f(\theta )+\lambda |\theta |$。如果$f(x)$ 可导，且$▽f(x)$满足L-Lipschitz条件，则有： 存在常数$L>0$:
$$|▽f(x_1)-▽f(x_2)|\le L|x_1-x_2|$$$$\to$$$$\frac{|▽f(x_1)-▽f(x_2)|}{|x_1-x_2|}\le L$$$$\to$$$$\underset{x_1\to x_2}{\lim }\frac{|▽f(x_1)-▽f(x_2)|}{|x_1-x_2|}={▽}^{2}f(x)\le L$$
如果我们把$f(x)$在$x_k$处泰勒展开:
$$f(x)=f(x_k)+▽f(x_k)(x-x_k)+\frac{{▽}^{2}f(x_k)}{2}(x-x_k{)}^2$$
$$\le f(x_k)+▽f(x_k)(x-x_k)+\frac{L}{2}(x-x_k{)}^2$$$$=\frac{L}{2}{x}^2-L{x}_{k}x+\frac{L}{2}{x}_k^2+▽f(x_k)x-▽f(x_k)x_k+f(x_k)$$
$$=\frac{L}{2}(x-(x_k-\frac{1}{L}▽f(x_k)){)}^2+const$$
所以有$$f(x)+\lambda |x|\le \frac{L}{2}(x-(x_k-\frac{1}{L}▽f(x_k)){)}^2+\lambda |x|+const$$
于是我们只要最小化小于等于号右边就可以把符号左边的式子控制在一个相对较小的范围，所以我们可以认为最小化右边近似最小化等式左边。所以问题转化为：
$$x_{k+1}={\arg }_x\min \frac{L}{2}(x-(x_k-\frac{1}{L}▽f(x_k)){)}^2+\lambda |x|$$
我们令$z=x_k-\frac{1}{L}▽f(x_k)$,所以
$$x_{k+1}={\arg }_x\min \frac{L}{2}(x-z{)}^2+\lambda |x|$$
右边函数可以转化为：
当$x\ge 0$,$g_1(x)=\frac{L}{2}(x-z{)}^2+\lambda x=\frac{L}{2}(x-\frac{Lz-\lambda }{L})+cons{t}_1$;
当$x<0$,$g_2(x)=\frac{L}{2}(x-z{)}^2-\lambda x=\frac{L}{2}(x-\frac{Lz+\lambda }{L})+cons{t}_2$
由于函数不存在$x^i{x}^j$项，所以$x$各分量互不影响。
这显然是一个分段函数，且两端函数的交点一定相交于$y$轴。下面为$g_1$、$g_2$函数图像：

红色曲线为$g_1$,蓝色为$g_2$。因为$\frac{Lz-\lambda }{L}＜\frac{Lz+\lambda }{L}$。所以如果$\frac{Lz-\lambda }{L}＞0$，则两个函数的最小值对应的$x$都大于0，去除$g_1$小于0的部分，$g_2$大于0的部分，所以只能取到$g_1$的最小值。也就是，当$\frac{Lz-\lambda }{L}＞0\to z>\frac{\lambda }{L}$，$x=z-\frac{\lambda }{L}$;同理，$z＜\frac{\lambda }{L}$，$x=z+\frac{\lambda }{L}$；当$\frac{Lz-\lambda }{L}\le 0$且$\frac{Lz+\lambda }{L}\ge 0$，即$|z|\le \frac{\lambda }{L}$时，这时两个函数的最小值都取不到，这时最小值在$x$=0取得。我们把得到的结论写成分量形式：

从上式可以看出，如果$\lambda /L＜z^i$,最终会减去$\lambda /L$，经过迭代，最终一定会使$z^k\le \lambda /L$，最终一定会有$x^k$为0；$z^i＜-\lambda /L$时情况一定相同。最终可以得到稀疏解。

三、前向逐步回归
前向逐步算法开始将所有权重设置为1，然后每一步所做的决策是对某个权重增加或减少一个很小的值。具体算法的伪代码如下：

算法的代码部分如下：

def stageWise_regression(X,y,NumOfIt,episilon):
    X = (X-np.mean(X,axis=0)) / np.std(X,axis=0)
    y = y - np.mean(y)
    m,d = X.shape
    theta = np.zeros(d)
    thetaMat = np.zeros((NumOfIt,d))
    for i in range(NumOfIt):
        lowestError = 100000
        for j in range(d):
            for sign in [-1,1]:
                thetaTest=theta.copy()
                thetaTest[j] += sign * episilon
                Error = cost(thetaTest,X,y)
                if Error < lowestError:
                    lowestError=Error
                    thetaBest = thetaTest
        theta = thetaBest
        thetaMat[i,:]=theta
    return thetaMat
                
            



def cost(theta,X,y):
    m,d=X.shape
    cost=(X @ theta - y).T @ (X @ theta - y)/(2*m)
    return cost

这实际上是一种贪心算法，因为实际上每轮迭代中我们都会重新设置一个lowestError，也就是说，我们不能保证这一轮迭代的误差比上一轮迭代跟小，但是可以保证在这一轮迭代中我们的误差最小。我们再次使用鲍鱼数据及其处理数据函数进行试验，但是这里我们已经在stageWise_regression函数中已经将数据归一化，所以可以去掉

def prepared_data(path):
    file=open(path)
    fr_list = file.readlines()
    data_list=[]
    for data in fr_list:
        data = data.strip()
        data = data.split(',')
        if data[0] == 'F':
            data[0] = 1
        elif data[0] == 'I':
            data[0] = 0
        else:
            data[0] =-1
        data = [float(i) for i in data]
        data_list.append(data)
        data_mat = np.array(data_list)
        X = data_mat[:,:-1]
        y = data_mat[:,-1]
  return X,y
  
  thetaMat=stageWise_regression(X,y,5000,0.001)
print(thetaMat)

得到：

[[ 0.000e+00  0.000e+00  0.000e+00 ...  0.000e+00  0.000e+00  1.000e-03]
 [ 0.000e+00  0.000e+00  0.000e+00 ...  0.000e+00  0.000e+00  2.000e-03]
 [ 0.000e+00  0.000e+00  0.000e+00 ...  0.000e+00  0.000e+00  3.000e-03]
 ...
 [ 0.000e+00  0.000e+00  6.440e-01 ... -1.552e+00  0.000e+00  2.351e+00]
 [ 0.000e+00  0.000e+00  6.440e-01 ... -1.553e+00  0.000e+00  2.351e+00]
 [ 0.000e+00  0.000e+00  6.450e-01 ... -1.553e+00  0.000e+00  2.351e+00]]

（这里数据归一化函数为$x_{norm}=\frac{x-x_{mean}}{\sigma }$,这里的$\sigma$为标准差，而很多地方这里的分子为方差${\sigma }^2$,我觉得应该是错误的）