最小二乘法
前提
假设数据集D = { ( x 1 , y 1 ) , ( x 2 , y 2 ) , ⋯ , ( x N , y N ) } \mathcal{D}=\{(x_1, y_1),(x_2, y_2),\cdots,(x_N, y_N)\} D = {( x 1 , y 1 ) , ( x 2 , y 2 ) , ⋯ , ( x N , y N )} x i ∈ R , y i ∈ R , i = 1 , 2 , ⋯ , N x_i\in\mathbb{R},y_i\in\mathbb{R},i=1,2,\cdots,N x i ∈ R , y i ∈ R , i = 1 , 2 , ⋯ , N
X = ( x 1 , x 2 , ⋯ , x N ) T = [ x 1 T x 2 T ⋮ x N T ] N × p = [ x 11 x 12 ⋯ x 1 p x 21 x 22 ⋯ x 2 p ⋮ ⋮ ⋱ ⋮ x N 1 x N 2 ⋯ x N p ] N × p \begin{align} X&=(x_1,x_2,\cdots,x_N)^T \\ &= \begin{bmatrix} x_{1}^T \\x_{2}^T \\\vdots \\x_{N}^T \\ \end{bmatrix}_{N\times p} \\ & = \begin{bmatrix}x_{11} & x_{12} & \cdots & x_{1p} \\x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\x_{N1} & x_{N2} & \cdots & x_{Np} \\ \end{bmatrix}_{N\times p} \\ \end{align} X = ( x 1 , x 2 , ⋯ , x N ) T = x 1 T x 2 T ⋮ x N T N × p = x 11 x 21 ⋮ x N 1 x 12 x 22 ⋮ x N 2 ⋯ ⋯ ⋱ ⋯ x 1 p x 2 p ⋮ x Np N × p Y = ( y 1 , y 2 , ⋯ , y N ) T = [ y 1 T y 2 T ⋮ y N T ] N × p = [ y 11 y 12 ⋯ y 1 p y 21 y 22 ⋯ y 2 p ⋮ ⋮ ⋱ ⋮ y N 1 y N 2 ⋯ y N p ] N × p \begin{align} Y&=(y_1,y_2,\cdots,y_N)^T \\ &= \begin{bmatrix} y_{1}^T \\ y_{2}^T \\ \vdots \\ y_{N}^T \\ \end{bmatrix}_{N\times p} \\ & = \begin{bmatrix} y_{11} & y_{12} & \cdots & y_{1p} \\ y_{21} & y_{22} & \cdots & y_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ y_{N1} & y_{N2} & \cdots & y_{Np} \\ \end{bmatrix}_{N\times p} \\ \end{align} Y = ( y 1 , y 2 , ⋯ , y N ) T = y 1 T y 2 T ⋮ y N T N × p = y 11 y 21 ⋮ y N 1 y 12 y 22 ⋮ y N 2 ⋯ ⋯ ⋱ ⋯ y 1 p y 2 p ⋮ y Np N × p 线性回归假设：
f ( w ) = w T x f(w)=w^Tx f ( w ) = w T x 构造数据
X = [ x T , 1 ] T X=[x^T,1]^T X = [ x T , 1 ] T y = 3 x + 2 y=3x+2 y = 3 x + 2 x x x w = [ 3 , 2 ] T w=[3,2]^T w = [ 3 , 2 ] T X . d o t ( w ) X.dot(w) X . d o t ( w )
复制 import numpy as np
"""
构造数据
"""
X = np.linspace(0,100,100) # 生成100个样本数据
X = np.c_[X,np.ones(100)]
w = np.asarray([3,2])
Y = X.dot(w) # 矩阵点乘
X = X.astype('float')
Y = Y.astype('float')
X[:,0] += np.random.normal(size=(X[:,0].shape))*3 #添加噪声
Y = Y.reshape(Y.size,1) x x x
X X X
绘制图形
复制 import matplotlib.pyplot as plt
# 绘制(x,y)
plt.scatter(X[:,0],Y)
plt.plot(np.arange(0,100).reshape((100,1)),Y,'y.')
plt.xlabel("X")
plt.ylabel("Y") 如图：
定义
直接求闭式解
采用二范数定义的平方误差来定义损失函数：
l o s s f u n c t i o n ： L ( w ) = ∑ i = 1 N ∣ ∣ w T x i − y i ∣ ∣ 2 2 = ∑ i = 1 N ( w T x i − y i ) 2 = ( w T x 1 − y 1 , ⋯ , w T x N − y N ) ⋅ ( w T x 1 − y 1 , ⋯ , w T x N − y N ) T = ( w T X T − Y T ) ⋅ ( X w − Y ) = w T X T X w − Y T X w − w T X T Y + Y T Y = w T X T X w − 2 w T X T Y + Y T Y \begin{align} loss function： L(w)&=\sum\limits_{i=1}^N||w^Tx_i-y_i||^2_2 \\ & = \sum\limits_{i=1}^N(w^Tx_i-y_i)^2 \\ &=(w^Tx_1-y_1,\cdots,w^Tx_N-y_N)\cdot (w^Tx_1-y_1,\cdots,w^Tx_N-y_N)^T\nonumber \\ &=(w^TX^T-Y^T)\cdot (Xw-Y) \\ &=w^TX^TXw-Y^TXw-w^TX^TY+Y^TY\nonumber \\ &=w^TX^TXw-2w^TX^TY+Y^TY \end{align} l oss f u n c t i o n ： L ( w ) = i = 1 ∑ N ∣∣ w T x i − y i ∣ ∣ 2 2 = i = 1 ∑ N ( w T x i − y i ) 2 = ( w T x 1 − y 1 , ⋯ , w T x N − y N ) ⋅ ( w T x 1 − y 1 , ⋯ , w T x N − y N ) T = ( w T X T − Y T ) ⋅ ( Xw − Y ) = w T X T Xw − Y T Xw − w T X T Y + Y T Y = w T X T Xw − 2 w T X T Y + Y T Y ( w T x i − y i ) 2 = ( Y − X w ) T ( Y − X w ) (w^Tx_i-y_i)^2=(Y−Xw)^T(Y−Xw) ( w T x i − y i ) 2 = ( Y − Xw ) T ( Y − Xw )
最小化值w ^ \hat{w} w ^
w ^ = a r g m i n w L ( w ) ⟶ ∂ ∂ w L ( w ) = 0 ⟶ 2 X T X w ^ − 2 X T Y = 0 ⟶ w ^ = ( X T X ) − 1 X T Y = X + Y \begin{align} \hat{w}=\mathop{argmin}\limits_wL(w)&\longrightarrow\frac{\partial}{\partial w}L(w)=0\nonumber\\ &\longrightarrow2X^TX\hat{w}-2X^TY=0\nonumber\\ &\longrightarrow \hat{w}=(X^TX)^{-1}X^TY=X^+Y \end{align} w ^ = w a r g min L ( w ) ⟶ ∂ w ∂ L ( w ) = 0 ⟶ 2 X T X w ^ − 2 X T Y = 0 ⟶ w ^ = ( X T X ) − 1 X T Y = X + Y ( X T X ) − 1 X T (X^TX)^{-1}X^T ( X T X ) − 1 X T X + X^+ X +
这个式子中 ( X T X ) − 1 X T (X^TX)^{-1}X^T ( X T X ) − 1 X T X X X X X X
X = U Σ V T X=U\Sigma V^T X = U Σ V T 于是：
X + = V Σ − 1 U T X^+=V\Sigma^{-1}U^T X + = V Σ − 1 U T 在几何上，最小二乘法相当于模型（这里就是直线）和试验值的距离的平方求和，假设我们的试验样本张成一个 $p$ 维空间（满秩的情况）：X = S p a n ( x 1 , ⋯ , x N ) X=Span(x_1,\cdots,x_N) X = Sp an ( x 1 , ⋯ , x N ) f ( w ) = w T x = X T β f(w)=w^Tx=X^T\beta f ( w ) = w T x = X T β x 1 , ⋯ , x N x_1,\cdots,x_N x 1 , ⋯ , x N Y Y Y
X T ⋅ ( Y − X β ) = 0 ⟶ β = ( X T X ) − 1 X T Y X^T\cdot(Y-X\beta)=0\longrightarrow\beta=(X^TX)^{-1}X^TY X T ⋅ ( Y − Xβ ) = 0 ⟶ β = ( X T X ) − 1 X T Y 提示：a ⃗ ⊥ b ⃗ = > a T ⋅ b = 0 \vec a \perp\vec b => a^T \cdot b = 0 a ⊥ b => a T ⋅ b = 0
通过伪逆求解到的结果有如下优点：
（1）当w w w w = X + Y w=X^+Y w = X + Y ∣ ∣ w ∣ ∣ 2 ||w||^2 ∣∣ w ∣ ∣ 2
（2）当w w w w w w X w Xw Xw Y Y Y ( w T x i − y i ) 2 (w^Tx_i-y_i)^2 ( w T x i − y i ) 2
复制 # np.linalg.inv(X) <==> np.matrix(X).I X X X
Y Y Y
求出参数w w w
复制 np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y) 求出Y 1 Y1 Y 1
复制 w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y)
Y1 = X.dot(w)
Y1[:15] 绘图：
复制 # 绘制(x,y)
plt.scatter(X[:,0],Y)
plt.plot(X[:,0],Y1,'y')
plt.xlabel("X")
plt.ylabel("Y") 梯度下降求解
随机梯度下降法
但对于数据量很大的情况，求闭式解的方式会让内存很吃力，我们可以通过随机梯度下降法（SGD）对w w w w w w w w w
w : = w − η d L d w ，其中 η 参数取值一般很小，比如 0.0001 w:=w-\eta \frac{dL}{dw}，其中\eta 参数取值一般很小，比如0.0001 w := w − η d w d L ，其中 η 参数取值一般很小，比如 0.0001 其中d L d w = 2 X T X w ^ − 2 X T Y = 2 X T ( X w ^ − Y ) （ 1 ） \frac{dL}{dw}=2X^TX\hat{w}-2X^TY=2X^T(X\hat{w}-Y) （1） d w d L = 2 X T X w ^ − 2 X T Y = 2 X T ( X w ^ − Y ) （ 1 ） ( w T x i − y i ) 2 = ( Y − X w ) T ( Y − X w ) （ 2 ） (w^Tx_i-y_i)^2=(Y−Xw)^T(Y−Xw) （2） ( w T x i − y i ) 2 = ( Y − Xw ) T ( Y − Xw ) （ 2 ）
复制 # 初始化
w = np.random.random(size=(2,1))
# 迭代次数
epoches = 100
# eta参数
eta = 0.0000001
#记录loss变化(损失)
losses = []
for e in range(epoches):
nw = 2*X.T.dot(X.dot(w)-Y) #（1）
w = w - eta*nw
# reshape(-1)：变成一维数组
losses.append((Y-X.dot(w)).T.dot(Y-X.dot(w)).reshape(-1)) # （2）
w 参数η \eta η
求出w w w
可视化：
loss变化：
复制 plt.plot(losses)
plt.xlabel("iterations")
plt.ylabel("loss") 小批量梯度下降法
问题
在上面的梯度下降的例子中存在一个问题，w 1 w1 w 1 w 2 w2 w 2 w 1 比 w 2 w1比w2 w 1 比 w 2 w = [ w 1 , w 2 ] T w=[w1,w2]^T w = [ w 1 , w 2 ] T f ( x ) = x ∗ w 1 + 1 ⋅ w 2 f(x)=x∗w1+1⋅w2 f ( x ) = x ∗ w 1 + 1 ⋅ w 2 x x x f ( x ) f(x) f ( x ) w 1 w1 w 1 w 2 w2 w 2
（1）对输入X X X w 1 , w 2 w1,w2 w 1 , w 2
代码实现
复制 from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(X[:,:-1],Y) #训练数据
predict = lr.predict(X[:,:-1])
#查看w系数,b常量
print('w:',lr.coef_,'b:',lr.intercept_)
#查看标准差
np.std(Y-predict) 假设函数为y = w x + ϵ y=wx+\epsilon y = w x + ϵ
lr.intercept_:为ϵ \epsilon ϵ
噪声为高斯分布的 MLE（极大似然估计）
对于一维的情况，记 y = w T x + ϵ , ϵ ∼ N ( 0 , σ 2 ) y=w^Tx+\epsilon,\epsilon\sim\mathcal{N}(0,\sigma^2) y = w T x + ϵ , ϵ ∼ N ( 0 , σ 2 ) y ∼ N ( w T x , σ 2 ) y\sim\mathcal{N}(w^Tx,\sigma^2) y ∼ N ( w T x , σ 2 )
其中正态分布函数的概率密度为p ( y ∣ x i w ) = 1 2 π σ e − ( y i − w T x i ) 2 2 σ 2 p(y|x_iw)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y_i-w^Tx_i)^2}{2\sigma^2}} p ( y ∣ x i w ) = 2 π σ 1 e − 2 σ 2 ( y i − w T x i ) 2
L ( w ) = log p ( Y ∣ X , w ) = log ∏ i = 1 N p ( y i ∣ x i , w ) = ∑ i = 1 N log ( 1 2 π σ e − ( y i − w T x i ) 2 2 σ 2 ) = ∑ i = 1 N ( log 1 2 π σ + log e − ( y i − w T x i ) 2 2 σ 2 ) = ∑ i = 1 N ( log 1 2 π σ − ( y i − w T x i ) 2 2 σ 2 ) a r g m a x w L ( w ) = a r g m a x w ∑ i = 1 N − 1 2 σ 2 ( y i − w T x i ) 2 = a r g m i n w ∑ i = 1 N ( y i − w T x i ) 2 ( 进而转化为求最小 ) \begin{align} L(w)=\log p(Y|X,w)&=\log\prod\limits_{i=1}^Np(y_i|x_i,w)\nonumber\\ &=\sum\limits_{i=1}^N\log(\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y_i-w^Tx_i)^2}{2\sigma^2}})\\ &=\sum\limits_{i=1}^N\big ( \log\frac{1}{\sqrt{2\pi}\sigma}+\log e^{-\frac{(y_i-w^Tx_i)^2}{2\sigma^2}} \big )\\ &=\sum\limits_{i=1}^N \big ( \log\frac{1}{\sqrt{2\pi}\sigma}- \frac{(y_i-w^Tx_i)^2}{2\sigma^2} \big )\\ \mathop{argmax}\limits_wL(w) &=\mathop{argmax}\limits_w\sum \limits_{i=1}^N -\frac 1{2\sigma^2} (y_i-w^Tx_i)^2\\ &=\mathop{argmin}\limits_w\sum\limits_{i=1}^N (y_i-w^Tx_i)^2 (进而转化为求最小) \end{align} L ( w ) = log p ( Y ∣ X , w ) w a r g ma x L ( w ) = log i = 1 ∏ N p ( y i ∣ x i , w ) = i = 1 ∑ N log ( 2 π σ 1 e − 2 σ 2 ( y i − w T x i ) 2 ) = i = 1 ∑ N ( log 2 π σ 1 + log e − 2 σ 2 ( y i − w T x i ) 2 ) = i = 1 ∑ N ( log 2 π σ 1 − 2 σ 2 ( y i − w T x i ) 2 ) = w a r g ma x i = 1 ∑ N − 2 σ 2 1 ( y i − w T x i ) 2 = w a r g min i = 1 ∑ N ( y i − w T x i ) 2 ( 进而转化为求最小 ) LSE（最小二乘估计）<==> MLE（极大似然估计），其中噪声为高斯分布
（正则化）Regularized LSE <==> MAP，其中噪声为高斯分布
权重先验为高斯分布的 MAP
取先验分布 w ∼ N ( 0 , σ 0 2 ) w\sim\mathcal{N}(0,\sigma_0^2) w ∼ N ( 0 , σ 0 2 )
概率密度为p ( w ) = 1 2 π σ e − w 2 2 σ 0 2 p(w)=\frac{1}{\sqrt{2\pi\sigma}}e^{-\frac{w^2}{2\sigma_0^2}} p ( w ) = 2 πσ 1 e − 2 σ 0 2 w 2 p ( y ∣ w ) = 1 2 π σ e − ( y − w T x ) 2 2 σ 2 p(y|w)=\frac{1}{\sqrt{2\pi\sigma}}e^{-\frac{(y-w^Tx)^2}{2\sigma^2}} p ( y ∣ w ) = 2 πσ 1 e − 2 σ 2 ( y − w T x ) 2
有p ( y ∣ w ) ⋅ p ( w ) = 1 2 π σ 0 ⋅ 1 2 π σ ⋅ e − ( y − w T x ) 2 2 σ 2 ⋅ e − ∣ ∣ w ∣ ∣ 2 2 σ 2 p(y|w) \cdot p(w) =\frac{1}{\sqrt{2\pi}\sigma_0}\cdot \frac{1}{\sqrt{2\pi}\sigma} \cdot e^{-\frac{(y-w^Tx)^2}{2\sigma^2}} \cdot e^{-\frac{||w||^2}{2\sigma^2}} p ( y ∣ w ) ⋅ p ( w ) = 2 π σ 0 1 ⋅ 2 π σ 1 ⋅ e − 2 σ 2 ( y − w T x ) 2 ⋅ e − 2 σ 2 ∣∣ w ∣ ∣ 2
w ^ = a r g m a x w p ( w ∣ y ) = a r g m a x w p ( y ∣ w ) p ( w ) = a r g m a x w log p ( y ∣ w ) p ( w ) = a r g m a x w ( log p ( y ∣ w ) + log p ( w ) ) = a r g m i n w [ ( y − w T x ) 2 2 σ 2 + ∣ ∣ w ∣ ∣ 2 2 σ 0 2 ] = a r g m i n w [ ( y − w T x ) 2 + σ 2 σ 0 2 w T w ] \begin{align} \hat{w}&=\mathop{argmax}\limits_wp(w|y) \\ &=\mathop{argmax}\limits_wp(y|w)p(w)\nonumber\\ &=\mathop{argmax}\limits_w\log p(y|w)p(w)\nonumber\\ &=\mathop{argmax}\limits_w(\log p(y|w)+\log p(w))\nonumber\\ &=\mathop{argmin}\limits_w[\frac {(y-w^Tx)^2} {2\sigma^2} +\frac {||w||^2} {2\sigma_0^2}] \\ &=\mathop{argmin}\limits_w[(y-w^Tx)^2+\frac{\sigma^2}{\sigma_0^2}w^Tw]\\ \end{align} w ^ = w a r g ma x p ( w ∣ y ) = w a r g ma x p ( y ∣ w ) p ( w ) = w a r g ma x log p ( y ∣ w ) p ( w ) = w a r g ma x ( log p ( y ∣ w ) + log p ( w )) = w a r g min [ 2 σ 2 ( y − w T x ) 2 + 2 σ 0 2 ∣∣ w ∣ ∣ 2 ] = w a r g min [( y − w T x ) 2 + σ 0 2 σ 2 w T w ] 这里省略了 $X$，$p(Y)$和 $w$ 没有关系，同时也利用了上面高斯分布的 MLE的结果。
将会看到，超参数 $\sigma_0$的存在和下面会介绍的 Ridge 正则项可以对应，同样的如果将先验分布取为 Laplace 分布，那么就会得到和 L1 正则类似的结果。
正则化
过拟合
演示
复制 import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
#造伪样本
X=np.linspace(0,100,100)
X=np.c_[X,np.ones(100)]
w=np.asarray([3,2])
Y=X.dot(w)
X=X.astype('float')
Y=Y.astype('float')
X[:,0]+=np.random.normal(size=(X[:,0].shape))*3#添加噪声
Y=Y.reshape(100,1) # 一维数组转变为二维矩阵
复制 #拟合数据并可视化
w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y)
Y1 = X.dot(w)
# 绘制(x,y)
plt.scatter(X[:,0],Y)
plt.plot(X[:,0],Y1,'y')
plt.xlabel("X")
plt.ylabel("Y")
复制 # 初始化
w = np.random.random(size=(2,1))
# 迭代次数
epoches = 100
# eta参数
eta = 0.0000001
#记录loss变化(损失)
losses = []
for e in range(epoches):
nw = 2*X.T.dot(X.dot(w)-Y) #（1）
w = w - eta*nw
# reshape(-1)：变成一维数组
losses.append((Y-X.dot(w)).T.dot(Y-X.dot(w)).reshape(-1)) # （2）
plt.plot(losses)
复制 # np.concatenate(x,y)将x和y两个矩阵垂直拼接在一起
X=np.concatenate([X,np.array([[100,1],[101,1],[102,1],[103,1],[104,1]])])
Y=np.concatenate([Y,np.array([[3000],[3300],[3600],[3800],[3900]])])
#拟合数据并可视化
w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y)
Y1 = X.dot(w)
# 绘制(x,y)
plt.scatter(X[:,0],Y)
plt.plot(X[:,0],Y1,'y')
plt.xlabel("X")
plt.ylabel("Y") 因此在实际应用时，如果样本容量不远远大于样本的特征维度，很可能造成过拟合，对这种情况，有下面三个解决方式：
特征选择/特征提取（降低特征维度）如 PCA 算法。
正则化一般是在损失函数（如上面介绍的最小二乘损失）上加入正则化项（表示模型的复杂度对模型的惩罚），下面我们介绍一般情况下的两种正则化框架。
( l a s s o r e g r e s s i o n ) L 1 : a r g m i n w L ( w ) + λ ∣ ∣ w ∣ ∣ 1 , λ > 0 ( r i d g e r e g r e s s i o n / 岭回归 / 权值衰减 ) L 2 : a r g m i n w L ( w ) + λ ∣ ∣ w ∣ ∣ 2 2 , λ > 0 \begin{align} (lasso \quad regression) \quad L1&:\mathop{argmin}\limits_wL(w)+\lambda||w||_1,\lambda\gt0\\ (ridge \quad regression/岭回归/权值衰减) \quad L2&:\mathop{argmin}\limits_wL(w)+\lambda||w||^2_2,\lambda \gt 0 \end{align} ( l a sso re g ress i o n ) L 1 ( r i d g e re g ress i o n / 岭回归 / 权值衰减 ) L 2 : w a r g min L ( w ) + λ ∣∣ w ∣ ∣ 1 , λ > 0 : w a r g min L ( w ) + λ ∣∣ w ∣ ∣ 2 2 , λ > 0 ∣ ∣ w ∣ ∣ 2 = w T ⋅ w ||w||^2 = w^T \cdot w ∣∣ w ∣ ∣ 2 = w T ⋅ w
$\text{当}||w|| = |w_1| + |w_2| + ··· + |w_n| \text{，称之为L1正则化}$。
当 ∣ ∣ w ∣ ∣ = w 1 2 + w 2 2 + ⋅ ⋅ ⋅ + w n 2 \text{当}||w|| = w_1^2 + w_2^2 + ··· + w_n^2 当 ∣∣ w ∣∣ = w 1 2 + w 2 2 + ⋅⋅⋅ + w n 2
这里∣ ∣ w ∣ ∣ ||w|| ∣∣ w ∣∣
L1 Lasso
L1正则化可以引起稀疏解。
从最小化损失的角度看，由于 L1 项求导在0附近的左右导数都不是0，因此更容易取到0解。
从另一个方面看，L1 正则化相当于：
a r g m i n w L ( w ) s . t . ∣ ∣ w ∣ ∣ 1 < C w ^ = ( X T X ) − 1 ( X T Y − λ I 2 ) （后面会用这个） \mathop{argmin}\limits_wL(w)\\ s.t. ||w||_1\lt C\\ \hat{w}=(X^TX)^{-1}(X^TY-\frac{\lambda I}2)（后面会用这个） w a r g min L ( w ) s . t .∣∣ w ∣ ∣ 1 < C w ^ = ( X T X ) − 1 ( X T Y − 2 λ I ) （后面会用这个） 我们已经看到平方误差损失函数在 $w$ 空间是一个椭球，因此上式求解就是椭球和 $||w||_1=C$的切点，因此更容易相切在坐标轴上。
下面分别对LASSO的cost function的两部分求解:
R S S ( w ) = ∑ i = 1 m ( w T x i − y i ) 2 = ∑ i = 1 m ( ∑ j = 1 n x i j w j − y i ) 2 \begin{align} RSS(w) &= \sum_{i=1}^m(w^Tx_i-y_i)^2 \\ &= \sum_{i=1}^m(\sum_{j=1}^n{x_{ij}w_j}-y_i)^2 \end{align} RSS ( w ) = i = 1 ∑ m ( w T x i − y i ) 2 = i = 1 ∑ m ( j = 1 ∑ n x ij w j − y i ) 2 求导：
( 对 w k 进行求导， w i ≠ k 为常数 ) (对w_k进行求导，w_{i≠k}为常数) ( 对 w k 进行求导， w i = k 为常数 )
∂ ∂ w k R S S ( w ) = 2 ∑ i = 1 m x i k ( ∑ j = 1 n x i j w j − y i ) = 2 ∑ i = 1 m ( x i k ∑ j = 1 n x i j w j − y i x i k ) = 2 ∑ i = 1 m ( x i k ∑ j = 1 , j ≠ k n x i j w j + x i k 2 w k − y i x i k ) = 2 ∑ i = 1 m x i k ( ∑ j = 1 , j ≠ k n x i j w j − y i ) + 2 w k ∑ i = 1 m x i k 2 \begin{align} \frac{\partial}{\partial w_k}RSS(w) &= 2\sum_{i=1}^m{x_{ik}}(\sum_{j=1}^nx_{ij}w_j-y_i) \\ &= 2\sum_{i=1}^m(x_{ik}\sum_{j=1}^nx_{ij}w_j-y_ix_{ik}) \\ &=2\sum_{i=1}^m(x_{ik}\sum_{j=1,j≠k}^nx_{ij}w_j+x_{ik}^2w_k-y_ix_{ik}) \\ &=2\sum_{i=1}^mx_{ik}(\sum_{j=1,j≠k}^nx_{ij}w_j-y_i)+2w_k\sum_{i=1}^m x_{ik}^2 \\ \end{align} ∂ w k ∂ RSS ( w ) = 2 i = 1 ∑ m x ik ( j = 1 ∑ n x ij w j − y i ) = 2 i = 1 ∑ m ( x ik j = 1 ∑ n x ij w j − y i x ik ) = 2 i = 1 ∑ m ( x ik j = 1 , j = k ∑ n x ij w j + x ik 2 w k − y i x ik ) = 2 i = 1 ∑ m x ik ( j = 1 , j = k ∑ n x ij w j − y i ) + 2 w k i = 1 ∑ m x ik 2 令 p k = ∑ i = 1 m x i k ( y i − ∑ j = 1 , j ≠ k n x i j w j ) , z k = ∑ i = 1 m x i k 2 令p_k=\sum_{i=1}^mx_{ik}(y_i-\sum_{j=1,j≠k}^nx_{ij}w_j),z_k=\sum_{i=1}^m x_{ik}^2 令 p k = i = 1 ∑ m x ik ( y i − j = 1 , j = k ∑ n x ij w j ) , z k = i = 1 ∑ m x ik 2 有，
∂ ∂ w k R S S ( w ) = − 2 p k + 2 z k w k \begin{align} \frac{\partial}{\partial w_k}RSS(w) &= -2p_k+2z_kw_k \end{align} ∂ w k ∂ RSS ( w ) = − 2 p k + 2 z k w k \begin{align} \lambda \frac{\partial \sum_{i=1}^n|w_i| }{\partial w_k} &= \left\{ \begin{array}{**lr**} -\lambda, & w_k<0\\ [-\lambda,\lambda], & w_k=0\\ \lambda, & w_k>0\\ \end{array} \right. \end{align}
进行整体的求导可得：
\begin{align} \frac{\partial}{\partial w_k}L(w) &=\frac{\partial}{\partial w_k}RSS(w) +\lambda \frac{\partial \sum_{i=1}^n|w_i| }{\partial w_k}\\ &= 2z_kw_k-2p_k+ \left\{ \begin{array}{**lr**} -\lambda, & w_k<0\\ [-\lambda,\lambda], & w_k=0\\ \lambda, & w_k>0\\ \end{array} \right. \\ &= \left\{ \begin{array}{**lr**} 2z_kw_k-2p_k-\lambda, & w_k<0\\ [-2p_k-\lambda,-2p_k+\lambda], & w_k=0\\ 2z_kw_k-2p_k+\lambda, & w_k>0\\ \end{array} \right. \\ & = 0 \end{align}
得到：
\begin{align} \hat{w_k} &= \left\{ \begin{array}{**lr**} (\lambda/2+p_k)/z_k, & p_k<-\lambda/2,\\ 0, & -\lambda/2 \leq p_k \leq \lambda/2\\ (p_k-\lambda/2)/z_k, & p_k>\lambda/2\\ \end{array} \right. \\ 其中，p_k&=\sum_{i=1}^mx_{ik}(y_i-\sum_{j=1,j≠k}^nx_{ij}w_j),z_k=\sum_{i=1}^m x_{ik}^2，\\& (w^Tx_i-y_i)^2=(Y−Xw)^T(Y−Xw) \end{align}
手动实现代码
复制 import numpy as np
import matplotlib.pyplot as plt
# 计算残差
def RSS(X,y,w):
return (Y-X.dot(w)).T.dot(Y-X.dot(w))
#造伪样本
X=np.linspace(0,100,100)
X=np.c_[X,np.ones(100)]
w=np.asarray([3,2])
print(w)
Y=X.dot(w)
X=X.astype('float')
Y=Y.astype('float')
X[:,0]+=np.random.normal(size=(X[:,0].shape))*3#添加噪声
Y=Y.reshape(Y.shape[0],1) # 一维数组转变为二维矩阵
# np.concatenate(x,y)将x和y两个矩阵垂直拼接在一起
X=np.concatenate([X,np.array([[100,1],[101,1],[102,1],[103,1],[104,1]])])
Y=np.concatenate([Y,np.array([[3000],[3300],[3600],[3800],[3900]])])
#拟合数据并可视化
lam = 1000000
i = np.eye(X.shape[1])
w = np.linalg.pinv(X.T.dot(X)).dot((X.T.dot(Y))-(lam*i)/2)
Y1=X.dot(w[:,0])
plt.scatter(X[:,:-1],Y)
plt.plot(X[:,:-1],Y1)
复制 import numpy as np
import matplotlib.pyplot as plt
# 计算残差
def RSS(X,y,w):
return (Y-X.dot(w)).T.dot(Y-X.dot(w))
#造伪样本
X=np.linspace(0,100,100)
X=np.c_[X,np.ones(100)]
w=np.asarray([3,2])
print(w)
Y=X.dot(w)
X=X.astype('float')
Y=Y.astype('float')
X[:,0]+=np.random.normal(size=(X[:,0].shape))*3#添加噪声
Y=Y.reshape(Y.shape[0],1) # 一维数组转变为二维矩阵
# np.concatenate(x,y)将x和y两个矩阵垂直拼接在一起
X=np.concatenate([X,np.array([[100,1],[101,1],[102,1],[103,1],[104,1]])])
Y=np.concatenate([Y,np.array([[3000],[3300],[3600],[3800],[3900]])])
#拟合数据并可视化
m,n = X.shape
lambda_params = 100
z = [0 for i in range(n)]
p = [0 for i in range(n)]
w = np.eye(n,1)
for k in range(n):
z[k] = np.sum(X[:,k]*X[:,k],axis=0)
p[k] = np.sum([X[i,k]* (Y[i]-np.sum([X[i,j]*w[j] for j in range(n)])) for i in range(m)],axis=0)
if p[k] > lambda_params/2:
w_k = (p[k]-lambda_params/2) / z[k]
elif p[k] < -lambda_params/2:
w_k = (lambda_params/2+p[k]) / z[k]
else:
w_k = 0
w[k] = w_k
y = X.dot(w)
plt.scatter(X[:,:-1],Y)
plt.plot(X[:,:-1],y)
w 使用sklearn
复制 from sklearn.linear_model import Lasso
from sklearn.preprocessing import StandardScaler
# 特征值和目标值是都必须进行标准化处理, 实例化两个标准化API
std_x = StandardScaler()
std_y = StandardScaler()
X = std_x.fit_transform(X)
Y = std_y.fit_transform(Y.reshape(-1,1))
lasso = Lasso()
lasso.fit(X[:,:-1],Y)
predict = lasso.predict(X[:,:-1])
plt.plot(X[:,:-1],predict)
plt.scatter(X[:,:-1],Y) L2 Ridge
s . t . ∣ ∣ w ∣ ∣ 2 2 < C s.t. ||w||_2^2\lt C\\ s . t .∣∣ w ∣ ∣ 2 2 < C w ^ = a r g m i n w L ( w ) + λ w T w ⟶ ∂ ∂ w L ( w ) + 2 λ w = 0 ⟶ 2 X T X w ^ − 2 X T Y + 2 λ w ^ = 0 ⟶ w ^ = ( X T X + λ I ) − 1 X T Y \begin{align} \hat{w}=\mathop{argmin}\limits_wL(w)+\lambda w^Tw&\longrightarrow\frac{\partial}{\partial w}L(w)+2\lambda w=0\nonumber\\ &\longrightarrow2X^TX\hat{w}-2X^TY+2\lambda \hat w=0\nonumber\\ &\longrightarrow \hat{w}=(X^TX+\lambda \mathbb{I})^{-1}X^TY \end{align} w ^ = w a r g min L ( w ) + λ w T w ⟶ ∂ w ∂ L ( w ) + 2 λ w = 0 ⟶ 2 X T X w ^ − 2 X T Y + 2 λ w ^ = 0 ⟶ w ^ = ( X T X + λ I ) − 1 X T Y 可以看到，这个正则化参数和前面的 MAP 结果不谋而合。利用2范数进行正则化不仅可以是模型选择 w w w X T X X^TX X T X
I 是一个单位矩阵 \mathbb{I}是一个单位矩阵 I 是一个单位矩阵
定义符号函数并进行向量化，用于对L1正则化项的梯度计算:
手动实现代码
复制 alpha = 100.0
w = np.linalg.pinv(
X.T.dot(X)+alpha*np.eye(X.T.dot(X).shape[0])
).dot(X.T.dot(Y))
y = X.dot(w)
plt.scatter(X[:,:-1],Y)
plt.plot(X[:,:-1],y) 使用sklearn
复制 from sklearn.linear_model import Ridge
from sklearn.preprocessing import StandardScaler
std_x = StandardScaler()
std_y = StandardScaler()
X2 = std_x.fit_transform(X)
Y2 = std_y.fit_transform(Y.reshape(-1,1))
ridge = Ridge(alpha=100)
ridge.fit(X2[:,:-1],Y2)
predict = ridge.predict(X2[:,:-1])
plt.plot(X2[:,:-1],predict)
plt.scatter(X2[:,:-1],Y2) L1与L2结合--ElasticNet回归
复制 from sklearn.linear_model import ElasticNet
std_x = StandardScaler()
std_y = StandardScaler()
X3 = std_x.fit_transform(X)
Y3 = std_y.fit_transform(Y.reshape(-1,1))
lasso = ElasticNet()
lasso.fit(X3[:,:-1],Y3)
predict = lasso.predict(X3[:,:-1])
plt.plot(X3[:,:-1],predict)
plt.scatter(X3[:,:-1],Y3) 当岭参数λ = 0 \lambda=0 λ = 0
当岭参数 λ \lambda λ w i w_i w i t t t
