这个实验的要求写的还是挺清楚的(与上学期相比),本博客采用python实现,科学计算库采用numpy
,作图采用matplotlib.pyplot
,为了简便在文件开头import如下:
import numpy as npimport matplotlib.pyplot as plt
本实验用到的numpy函数
一般把numpy
简写为np
(import numpy as np)。下面简单介绍一下实验中用到的numpy函数。下面的代码均需要在最前面加上import numpy as np
。
np.array
该函数返回一个numpy.ndarray
对象,可以理解为一个多维数组(本实验中仅会用到一维(可以当作列向量)和二维(矩阵))。下面用小写的 x \pmb x xx表示列向量,大写的 A A A表示矩阵。A.T
表示 A A A的转置。对ndarray
的运算一般都是逐元素的。
>>> x = np.array([1,2,3])>>> xarray([1, 2, 3])>>> A = np.array([[2,3,4],[5,6,7]])>>> Aarray([[2, 3, 4], [5, 6, 7]])>>> A.T # 转置array([[2, 5], [3, 6], [4, 7]])>>> A + 1array([[3, 4, 5], [6, 7, 8]])>>> A * 2array([[ 4, 6, 8], [10, 12, 14]])
np.random
np.random
模块中包含几个生成随机数的函数。在本实验中用随机初始化参数(梯度下降法),给数据添加噪声。
>>> np.random.rand(3, 3) # 生成3 * 3 随机矩阵,每个元素服从[0,1)均匀分布array([[8.18713933e-01, 5.46592778e-01, 1.36380542e-01], [9.85514865e-01, 7.07323389e-01, 2.51858374e-04], [3.14683662e-01, 4.74980699e-02, 4.39658301e-01]]) >>> np.random.rand(1) # 生成单个随机数array([0.70944563])>>> np.random.rand(5) # 长为5的一维随机数组array([0.03911319, 0.67572368, 0.98884287, 0.12501456, 0.39870096])>>> np.random.randn(3, 3) # 同上,但每个元素服从N(0, 1)(标准正态)
数学函数
本实验中只用到了np.sin
。这些数学函数是对np.ndarray
逐元素操作的:
>>> x = np.array([0, 3.1415, 3.1415 / 2]) # 0, pi, pi / 2>>> np.round(np.sin(x)) # 先求sin再四舍五入: 0, 0, 1array([0., 0., 1.])
此外,还有np.log
、np.exp
等与python的math
库相似的函数(只不过是对多维数组进行逐元素运算)。
np.dot
返回两个矩阵的乘积。与线性代数中的矩阵乘法一致。要求第一个矩阵的列等于第二个矩阵的行数。特殊地,当其中一个为一维数组时,形状会自动适配为 n × 1 n\times1 n×1或 1 × n . 1\times n. 1×n.
>>> x = np.array([1,2,3]) # 一维数组>>> A = np.array([[1,1,1],[2,2,2],[3,3,3]]) # 3 * 3矩阵>>> np.dot(x,A)array([14, 14, 14])>>> np.dot(A,x)array([ 6, 12, 18])>>> x_2D = np.array([[1,2,3]]) # 这是一个二维数组(1 * 3矩阵)>>> np.dot(x_2D, A) # 可以运算array([[14, 14, 14]])>>> np.dot(A, x_2D) # 行列不匹配Traceback (most recent call last): File "", line 1, in <module> File "", line 5, in dotValueError: shapes (3,3) and (1,3) not aligned: 3 (dim 1) != 1 (dim 0)
np.eye
np.eye(n)
返回一个n阶单位阵。
>>> A = np.eye(3)>>> Aarray([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])
线性代数相关
np.linalg
是与线性代数有关的库。
>>> Aarray([[1, 0, 0], [0, 2, 0], [0, 0, 3]])>>> np.linalg.inv(A) # 求逆(本实验不考虑逆不存在)array([[1. , 0. , 0. ], [0. , 0.5 , 0. ], [0. , 0. , 0.33333333]])>>> x = np.array([1,2,3])>>> np.linalg.norm(x) # 返回向量x的模长(平方求和开根号)3.7416573867739413>>> np.linalg.eigvals(A) # A的特征值array([1., 2., 3.])
生成数据
生成数据要求加入噪声(误差)。上课讲的时候举的例子就是正弦函数,我们这里也采用标准的正弦函数 y = sin x . y=\sin x. y=sinx.(加入噪声后即为 y = sin x + ϵ , y=\sin x+\epsilon, y=sinx+ϵ,其中 ϵ ∼ N ( 0 , σ 2 ) \epsilon\sim N(0, \sigma^2) ϵ∼N(0,σ2),由于 sin x \sin x sinx的最大值为 1 1 1,我们把误差的方差设小一点,这里设成 1 25 \frac{1}{25} 251)。
'''返回数据集,形如[[x_1, y_1], [x_2, y_2], ..., [x_N, y_N]]保证 bound[0] <= x_i < bound[1].- N 数据集大小, 默认为 100- bound 产生数据横坐标的上下界, 应满足 bound[0] < bound[1], 默认为(0, 10)'''def get_dataset(N = 100, bound = (0, 10)): l, r = bound # np.random.rand 产生[0, 1)的均匀分布,再根据l, r缩放平移 # 这里sort是为了画图时不会乱,可以去掉sorted试一试 x = sorted(np.random.rand(N) * (r - l) + l)# np.random.randn 产生N(0,1),除以5会变为N(0, 1 / 25) y = np.sin(x) + np.random.randn(N) / 5 return np.array([x,y]).T
产生的数据集每行为一个平面上的点。产生的数据看起来像这样:
隐隐约约能看出来是个正弦函数的形状。产生上面图像的代码如下:
dataset = get_dataset(bound = (-3, 3))# 绘制数据集散点图for [x, y] in dataset: plt.scatter(x, y, color = 'red')plt.show()
最小二乘法拟合
下面我们分别用四种方法(最小二乘,正则项/岭回归,梯度下降法,共轭梯度法)以用多项式拟合上述干扰过的正弦曲线。
解析解推导
简单回忆一下最小二乘法的原理:现在我们想用一个 m m m次多项式
f ( x ) = w 0 + w 1 x + w 2 x 2 + . . . + w m x m f(x)=w_0+w_1x+w_2x^2+…+w_mx^m f(x)=w0+w1x+w2x2+…+wmxm
来近似真实函数 y = sin x . y=\sin x. y=sinx.我们的目标是最小化数据集 ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( x N , y N ) (x_1,y_1),(x_2,y_2),…,(x_N,y_N) (x1,y1),(x2,y2),…,(xN,yN)上的损失 L L L(loss),这里损失函数采用平方误差:
L = ∑ i = 1 N [ y i − f ( x i ) ] 2 L=\sum\limits_{i=1}^N[y_i-f(x_i)]^2 L=i=1∑N[yi−f(xi)]2
为了求得使均方误差最小(因此最贴合目标曲线)的参数 w 0 , w 1 , . . . , w m , w_0,w_1,…,w_m, w0,w1,…,wm,我们需要分别求损失 L L L关于 w 0 , w 1 , . . . , w m w_0,w_1,…,w_m w0,w1,…,wm的导数。为了方便,我们采用线性代数的记法:
X = ( 1 x 1 x 1 2 ⋯ x 1 m 1 x 2 x 2 2 ⋯ x 2 m ⋮ ⋮ 1 x N x N 2 ⋯ x N m ) N × ( m + 1 ) , Y = ( y 1 y 2 ⋮ y N ) N × 1 , W = ( w 0 w 1 ⋮ w m ) ( m + 1 ) × 1 . X=\begin{pmatrix}1 & x_1 & x_1^2 & \cdots & x_1^m\\ 1 & x_2 & x_2^2 & \cdots & x_2^m\\ \vdots & & & &\vdots\\ 1 & x_N & x_N^2 & \cdots & x_N^m\\\end{pmatrix}_{N\times(m+1)},Y=\begin{pmatrix}y_1 \\ y_2 \\ \vdots \\y_N\end{pmatrix}_{N\times1},W=\begin{pmatrix}w_0 \\ w_1 \\ \vdots \\w_m\end{pmatrix}_{(m+1)\times1}. X=⎝ ⎛11⋮1x1x2xNx12x22xN2⋯⋯⋯x1mx2m⋮xNm⎠ ⎞N×(m+1),Y=⎝ ⎛y1y2⋮yN⎠ ⎞N×1,W=⎝ ⎛w0w1⋮wm⎠ ⎞(m+1)×1.
在这种表示方法下,有
( f ( x 1 ) f ( x 2 ) ⋮ f ( x N ) ) = X W . \begin{pmatrix}f(x_1)\\ f(x_2) \\ \vdots \\ f(x_N)\end{pmatrix}= XW. ⎝ ⎛f(x1)f(x2)⋮f(xN)⎠ ⎞=XW.
如果有疑问可以自己拿矩阵乘法验证一下。继续,误差项之和可以表示为
( f ( x 1 ) − y 1 f ( x 2 ) − y 2 ⋮ f ( x N ) − y N ) = X W − Y . \begin{pmatrix}f(x_1)-y_1 \\ f(x_2)-y_2 \\ \vdots \\ f(x_N)-y_N\end{pmatrix}=XW-Y. ⎝ ⎛f(x1)−y1f(x2)−y2⋮f(xN)−yN⎠ ⎞=XW−Y.
因此,损失函数
L = ( X W − Y ) T ( X W − Y ) . L=(XW-Y)^T(XW-Y). L=(XW−Y)T(XW−Y).
(为了求得向量 x = ( x 1 , x 2 , . . . , x N ) T \pmb x=(x_1,x_2,…,x_N)^T xx=(x1,x2,…,xN)T各分量的平方和,可以对 x \pmb x xx作内积,即 x T x . \pmb x^T \pmb x. xxTxx.)
为了求得使 L L L最小的 W W W(这个 W W W是一个列向量),我们需要对 L L L求偏导数,并令其为 0 : 0: 0:
∂ L ∂ W = ∂ ∂ W [ ( X W − Y ) T ( X W − Y ) ] = ∂ ∂ W [ ( W T X T − Y T ) ( X W − Y ) ] = ∂ ∂ W ( W T X T X W − W T X T Y − Y T X W + Y T Y ) = ∂ ∂ W ( W T X T X W − 2 Y T X W + Y T Y ) ( 容易验证 , W T X T Y = Y T X W , 因而可以将其合并 ) = 2 X T X W − 2 X T Y \begin{aligned}\frac{\partial L}{\partial W}&=\frac{\partial}{\partial W}[(XW-Y)^T(XW-Y)]\\ &=\frac{\partial}{\partial W}[(W^TX^T-Y^T)(XW-Y)] \\ &=\frac{\partial}{\partial W}(W^TX^TXW-W^TX^TY-Y^TXW+Y^TY)\\ &=\frac{\partial}{\partial W}(W^TX^TXW-2Y^TXW+Y^TY)(容易验证,W^TX^TY=Y^TXW,因而可以将其合并)\\ &=2X^TXW-2X^TY\end{aligned} ∂W∂L=∂W∂[(XW−Y)T(XW−Y)]=∂W∂[(WTXT−YT)(XW−Y)]=∂W∂(WTXTXW−WTXTY−YTXW+YTY)=∂W∂(WTXTXW−2YTXW+YTY)(容易验证,WTXTY=YTXW,因而可以将其合并)=2XTXW−2XTY
说明:
(1)从第3行到第4行,由于 W T X T Y W^TX^TY WTXTY和 Y T X W Y^TXW YTXW都是数(或者说 1 × 1 1\times1 1×1矩阵),二者互为转置,因此值相同,可以合并成一项。
(2)从第4行到第5行的矩阵求导,第一项 ∂ ∂ W ( W T ( X T X ) W ) \frac{\partial}{\partial W}(W^T(X^TX)W) ∂W∂(WT(XTX)W)是一个关于 W W W的二次型,其导数就是 2 X T X W . 2X^TXW. 2XTXW.
(3)对于一次项 − 2 Y T X W -2Y^TXW −2YTXW的求导,如果按照实数域的求导应该得到 − 2 Y T X . -2Y^TX. −2YTX.但检查一下发现矩阵的型对不上,需要做一下转置,变为 − 2 X T Y . -2X^TY. −2XTY.
矩阵求导线性代数课上也没有系统教过,只对这里出现的做一下说明。(多了我也不会 )
令偏导数为0,得到
X T X W = Y T X , X^TXW=Y^TX, XTXW=YTX,
左乘 ( X T X ) − 1 (X^TX)^{-1} (XTX)−1( X T X X^TX XTX的可逆性见下方的补充说明),得到
W = ( X T X ) − 1 X T Y . W=(X^TX)^{-1}X^TY. W=(XTX)−1XTY.
这就是我们想求的 W W W的解析解,我们只需要调用函数算出这个值即可。
'''最小二乘求出解析解, m 为多项式次数最小二乘误差为 (XW - Y)^T*(XW - Y)- dataset 数据集- m 多项式次数, 默认为 5'''def fit(dataset, m = 5): X = np.array([dataset[:, 0] ** i for i in range(m + 1)]).T Y = dataset[:, 1] return np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), Y)
稍微解释一下代码:第一行即生成上面约定的 X X X矩阵,dataset[:,0]
即数据集第0列 ( x 1 , x 2 , . . . , x N ) T (x_1,x_2,…,x_N)^T (x1,x2,…,xN)T;第二行即 Y Y Y矩阵;第三行返回上面的解析解。(如果不熟悉python语法或者numpy
库还是挺不友好的)
简单地验证一下我们已经完成的函数的结果:为此,我们先写一个draw
函数,用于把求得的 W W W对应的多项式 f ( x ) f(x) f(x)画到pyplot
库的图像上去:
'''绘制给定系数W的, 在数据集上的多项式函数图像- dataset 数据集- w 通过上面四种方法求得的系数- color 绘制颜色, 默认为 red- label 图像的标签'''def draw(dataset, w, color = 'red', label = ''): X = np.array([dataset[:, 0] ** i for i in range(len(w))]).T Y = np.dot(X, w) plt.plot(dataset[:, 0], Y, c = color, label = label)
然后是主函数:
if __name__ == '__main__': dataset = get_dataset(bound = (-3, 3)) # 绘制数据集散点图 for [x, y] in dataset: plt.scatter(x, y, color = 'red') # 最小二乘 coef1 = fit(dataset) draw(dataset, coef1, color = 'black', label = 'OLS') # 绘制图像 plt.legend() plt.show()
可以看到5次多项式拟合的效果还是比较不错的(数据集每次随机生成,所以跟第一幅图不一样)。
截至这部分全部的代码,后面同名函数不再给出说明:
import numpy as npimport matplotlib.pyplot as plt'''返回数据集,形如[[x_1, y_1], [x_2, y_2], ..., [x_N, y_N]]保证 bound[0] <= x_i < bound[1].- N 数据集大小, 默认为 100- bound 产生数据横坐标的上下界, 应满足 bound[0] < bound[1]'''def get_dataset(N = 100, bound = (0, 10)): l, r = bound x = sorted(np.random.rand(N) * (r - l) + l) y = np.sin(x) + np.random.randn(N) / 5 return np.array([x,y]).T'''最小二乘求出解析解, m 为多项式次数最小二乘误差为 (XW - Y)^T*(XW - Y)- dataset 数据集- m 多项式次数, 默认为 5'''def fit(dataset, m = 5): X = np.array([dataset[:, 0] ** i for i in range(m + 1)]).T Y = dataset[:, 1] return np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), Y)'''绘制给定系数W的, 在数据集上的多项式函数图像- dataset 数据集- w 通过上面四种方法求得的系数- color 绘制颜色, 默认为 red- label 图像的标签'''def draw(dataset, w, color = 'red', label = ''): X = np.array([dataset[:, 0] ** i for i in range(len(w))]).T Y = np.dot(X, w) plt.plot(dataset[:, 0], Y, c = color, label = label)if __name__ == '__main__': dataset = get_dataset(bound = (-3, 3)) # 绘制数据集散点图 for [x, y] in dataset: plt.scatter(x, y, color = 'red') coef1 = fit(dataset) draw(dataset, coef1, color = 'black', label = 'OLS') plt.legend() plt.show()
补充说明
上面有一块不太严谨:对于一个矩阵 X X X而言, X T X X^TX XTX不一定可逆。然而在本实验中,可以证明其为可逆矩阵。由于这门课不是线性代数课,我们就不费太多篇幅介绍这个了,仅作简单提示:
(1) X X X是一个 N × ( m + 1 ) N\times(m+1) N×(m+1)的矩阵。其中数据数 N N N远大于多项式次数 m m m,有 N > m + 1 ; N>m+1; N>m+1;
(2)为了说明 X T X X^TX XTX可逆,需要说明 ( X T X ) ( m + 1 ) × ( m + 1 ) (X^TX)_{(m+1)\times(m+1)} (XTX)(m+1)×(m+1)满秩,即 R ( X T X ) = m + 1 ; R(X^TX)=m+1; R(XTX)=m+1;
(3)在线性代数中,我们证明过 R ( X ) = R ( X T ) = R ( X T X ) = R ( X X T ) ; R(X)=R(X^T)=R(X^TX)=R(XX^T); R(X)=R(XT)=R(XTX)=R(XXT);
(4) X X X是一个范德蒙矩阵,由其性质可知其秩等于 m i n { N , m + 1 } = m + 1. min\{N,m+1\}=m+1. min{N,m+1}=m+1.
添加正则项(岭回归)
最小二乘法容易造成过拟合。为了说明这种缺陷,我们用所生成数据集的前50个点进行训练(这样抽样不够均匀,这里只是为了说明过拟合),得出参数,再画出整个函数图像,查看拟合效果:
if __name__ == '__main__': dataset = get_dataset(bound = (-3, 3)) # 绘制数据集散点图 for [x, y] in dataset: plt.scatter(x, y, color = 'red') # 取前50个点进行训练 coef1 = fit(dataset[:50], m = 3) # 再画出整个数据集上的图像 draw(dataset, coef1, color = 'black', label = 'OLS')
过拟合在 m m m较大时尤为严重(上面图像为 m = 3 m=3 m=3时)。当多项式次数升高时,为了尽可能贴近所给数据集,计算出来的系数的数量级将会越来越大,在未见样本上的表现也就越差。如上图,可以看到拟合在前50个点(大约在横坐标 [ − 3 , 0 ] [-3,0] [−3,0]处)表现很好;而在测试集上表现就很差( [ 0 , 3 ] [0,3] [0,3]处)。为了防止过拟合,可以引入正则化项。此时损失函数 L L L变为
L = ( X W − Y ) T ( X W − Y ) + λ ∣ ∣ W ∣ ∣ 2 2 L=(XW-Y)^T(XW-Y)+\lambda||W||_2^2 L=(XW−Y)T(XW−Y)+λ∣∣W∣∣22
其中 ∣ ∣ ⋅ ∣ ∣ 2 2 ||\cdot||_2^2 ∣∣⋅∣∣22表示 L 2 L_2 L2范数的平方,在这里即 W T W ; λ W^TW;\lambda WTW;λ为正则化系数。该式子也称岭回归(Ridge Regression)。它的思想是兼顾损失函数与所得参数 W W W的模长(在 L 2 L_2 L2范数时),防止 W W W内的参数过大。
举个例子(数是随便编的):当正则化系数为 1 1 1,若方案1在数据集上的平方误差为 0.5 , 0.5, 0.5,此时 W = ( 100 , − 200 , 300 , 150 ) T W=(100,-200,300,150)^T W=(100,−200,300,150)T;方案2在数据集上的平方误差为 10 , 10, 10,此时 W = ( 1 , − 3 , 2 , 1 ) W=(1,-3,2,1) W=(1,−3,2,1),那我们选择方案2的 W . W. W.正则化系数 λ \lambda λ刻画了这种对于 W W W模长的重视程度: λ \lambda λ越大,说明 W W W的模长升高带来的惩罚也就越大。当 λ = 0 , \lambda=0, λ=0,岭回归即变为普通的最小二乘法。与岭回归相似的还有LASSO,就是将正则化项换为 L 1 L_1 L1范数。
重复上面的推导,我们可以得出解析解为
W = ( X T X + λ E m + 1 ) − 1 X T Y . W=(X^TX+\lambda E_{m+1})^{-1}X^TY. W=(XTX+λEm+1)−1XTY.
其中 E m + 1 E_{m+1} Em+1为 m + 1 m+1 m+1阶单位阵。容易得到 ( X T X + λ E m + 1 ) (X^TX+\lambda E_{m+1}) (XTX+λEm+1)也是可逆的。
该部分代码如下。
'''岭回归求解析解, m 为多项式次数, l 为 lambda 即正则项系数岭回归误差为 (XW - Y)^T*(XW - Y) + λ(W^T)*W- dataset 数据集- m 多项式次数, 默认为 5- l 正则化参数 lambda, 默认为 0.5'''def ridge_regression(dataset, m = 5, l = 0.5): X = np.array([dataset[:, 0] ** i for i in range(m + 1)]).T Y = dataset[:, 1] return np.dot(np.dot(np.linalg.inv(np.dot(X.T, X) + l * np.eye(m + 1)), X.T), Y)
两种方法的对比如下:
对比可以看出,岭回归显著减轻了过拟合(此时为 m = 3 , λ = 0.3 m=3,\lambda=0.3 m=3,λ=0.3)。
梯度下降法
梯度下降法并不是求解该问题的最好方法,很容易就无法收敛。先简单介绍梯度下降法的基本思想:若我们想求取复杂函数 f ( x ) f(x) f(x)的最小值(最值点)(这个 x x x可能是向量等),即
x m i n = arg min x f ( x ) x_{min}=\argmin_{x}f(x) xmin=xargminf(x)
梯度下降法重复如下操作:
(0)(随机)初始化 x 0 ( t = 0 ) x_0(t=0) x0(t=0);
(1)设 f ( x ) f(x) f(x)在 x t x_t xt处的梯度(当 x x x为一维时,即导数) ∇ f ( x t ) \nabla f(x_t) ∇f(xt);
(2) x t + 1 = x t − η ∇ f ( x t ) x_{t+1}=x_t-\eta\nabla f(x_t) xt+1=xt−η∇f(xt)
(3)若 x t + 1 x_{t+1} xt+1与 x t x_t xt相差不大(达到预先设定的范围)或迭代次数达到预设上限,停止算法;否则重复(1)(2).
其中 η \eta η为学习率,它决定了梯度下降的步长。
下面是一个用梯度下降法求取 y = x 2 y=x^2 y=x2的最小值点的示例程序:
import numpy as npimport matplotlib.pyplot as pltdef f(x): return x ** 2def draw(): x = np.linspace(-3, 3) y = f(x) plt.plot(x, y, c = 'red')cnt = 0# 初始化 xx = np.random.rand(1) * 3learning_rate = 0.05while True: grad = 2 * x # -----------作图用,非算法部分----------- plt.scatter(x, f(x), c = 'black') plt.text(x + 0.3, f(x) + 0.3, str(cnt)) # ------------------------------------- new_x = x - grad * learning_rate # 判断收敛 if abs(new_x - x) < 1e-3: break x = new_x cnt += 1draw()plt.show()
上图标明了 x x x随着迭代的演进,可以看到 x x x不断沿着正半轴向零点靠近。需要注意的是,学习率不能过大(虽然在上面的程序中,学习率设置得有点小了),需要手动进行尝试调整,否则容易想象, x x x在正负半轴来回震荡,难以收敛。
在最小二乘法中,我们需要优化的函数是损失函数
L = ( X W − Y ) T ( X W − Y ) . L=(XW-Y)^T(XW-Y). L=(XW−Y)T(XW−Y).
下面我们用梯度下降法求解该问题。在上面的推导中,
∂ L ∂ W = 2 X T X W − 2 X T Y , \begin{aligned}\frac{\partial L}{\partial W}=2X^TXW-2X^TY\end{aligned}, ∂W∂L=2XTXW−2XTY,
于是我们每次在迭代中对 W W W减去该梯度,直到参数 W W W收敛。不过经过实验,平方误差会使得梯度过大,过程无法收敛,因此采用均方误差(MSE)替换之,就是给原来的式子除以 N N N:
'''梯度下降法(Gradient Descent, GD)求优化解, m 为多项式次数, max_iteration 为最大迭代次数, lr 为学习率注: 此时拟合次数不宜太高(m <= 3), 且数据集的数据范围不能太大(这里设置为(-3, 3)), 否则很难收敛- dataset 数据集- m 多项式次数, 默认为 3(太高会溢出, 无法收敛)- max_iteration 最大迭代次数, 默认为 1000- lr 梯度下降的学习率, 默认为 0.01'''def GD(dataset, m = 3, max_iteration = 1000, lr = 0.01): # 初始化参数 w = np.random.rand(m + 1) N = len(dataset) X = np.array([dataset[:, 0] ** i for i in range(len(w))]).T Y = dataset[:, 1] try: for i in range(max_iteration): pred_Y = np.dot(X, w) # 均方误差(省略系数2) grad = np.dot(X.T, pred_Y - Y) / N w -= lr * grad ''' 为了能捕获这个溢出的 Warning,需要import warnings并在主程序中加上: warnings.simplefilter('error') ''' except RuntimeWarning: print('梯度下降法溢出, 无法收敛') return w
这时如果 m m m设置得稍微大一点(比如4),在迭代过程中梯度就会溢出,使参数无法收敛。在收敛时,拟合效果还算可以:
共轭梯度法
共轭梯度法(Conjugate Gradients)可以用来求解形如 A x = b A\pmb x=\pmb b Axx=bb的方程组,或最小化二次型 f ( x ) = 1 2 x T A x − b T x + c . f(\pmb x)=\frac12\pmb x^TA\pmb x-\pmb b^T \pmb x+c. f(xx)=21xxTAxx−bbTxx+c.(可以证明对于正定的 A A A,二者等价)其中 A A A为正定矩阵。在本问题中,我们要求解 X T X W = Y T X , X^TXW=Y^TX, XTXW=YTX,
就有 A ( m + 1 ) × ( m + 1 ) = X T X , b = Y T . A_{(m+1)\times(m+1)}=X^TX,\pmb b=Y^T. A(m+1)×(m+1)=XTX,bb=YT.若我们想加一个正则项,就变成求解
( X T X + λ E ) W = Y T X . (X^TX+\lambda E)W=Y^TX. (XTX+λE)W=YTX.
首先说明一点: X T X X^TX XTX不一定是正定的但一定是半正定的(证明见此)。但是在实验中我们基本不用担心这个问题,因为 X T X X^TX XTX有极大可能是正定的,我们只在代码中加一个断言(assert),不多关注这个条件。
共轭梯度法的思想来龙去脉和证明过程比较长,可以参考这个系列,这里只给出算法步骤(在上面链接的第三篇开头):
(0)初始化 x ( 0 ) ; x_{(0)}; x(0);
(1)初始化 d ( 0 ) = r ( 0 ) = b − A x ( 0 ) ; d_{(0)}=r_{(0)}=b-Ax_{(0)}; d(0)=r(0)=b−Ax(0);
(2)令 α ( i ) = r ( i ) T r ( i ) d ( i ) T A d ( i ) ; \alpha_{(i)}=\frac{r_{(i)}^Tr_{(i)}}{d_{(i)}^TAd_{(i)}}; α(i)=d(i)TAd(i)r(i)Tr(i);
(3)迭代 x ( i + 1 ) = x ( i ) + α ( i ) d ( i ) ; x_{(i+1)}=x_{(i)}+\alpha_{(i)}d_{(i)}; x(i+1)=x(i)+α(i)d(i);
(4)令 r ( i + 1 ) = r ( i ) − α ( i ) A d ( i ) ; r_{(i+1)}=r_{(i)}-\alpha_{(i)}Ad_{(i)}; r(i+1)=r(i)−α(i)Ad(i);
(5)令 β ( i + 1 ) = r ( i + 1 ) T r ( i + 1 ) r ( i ) T r ( i ) , d ( i + 1 ) = r ( i + 1 ) + β ( i + 1 ) d ( i ) . \beta_{(i+1)}=\frac{r_{(i+1)}^Tr_{(i+1)}}{r_{(i)}^Tr_{(i)}},d_{(i+1)}=r_{(i+1)}+\beta_{(i+1)}d_{(i)}. β(i+1)=r(i)Tr(i)r(i+1)Tr(i+1),d(i+1)=r(i+1)+β(i+1)d(i).
(6)当 ∣ ∣ r ( i ) ∣ ∣ ∣ ∣ r ( 0 ) ∣ ∣ < ϵ \frac{||r_{(i)}||}{||r_{(0)}||}<\epsilon ∣∣r(0)∣∣∣∣r(i)∣∣<ϵ时,停止算法;否则继续从(2)开始迭代。 ϵ \epsilon ϵ为预先设定好的很小的值,我这里取的是 1 0 − 5 . 10^{-5}. 10−5.
下面我们按照这个过程实现代码:
'''共轭梯度法(Conjugate Gradients, CG)求优化解, m 为多项式次数- dataset 数据集- m 多项式次数, 默认为 5- regularize 正则化参数, 若为 0 则不进行正则化'''def CG(dataset, m = 5, regularize = 0): X = np.array([dataset[:, 0] ** i for i in range(m + 1)]).T A = np.dot(X.T, X) + regularize * np.eye(m + 1) assert np.all(np.linalg.eigvals(A) > 0), '矩阵不满足正定!' b = np.dot(X.T, dataset[:, 1]) w = np.random.rand(m + 1) epsilon = 1e-5 # 初始化参数 d = r = b - np.dot(A, w) r0 = r while True: alpha = np.dot(r.T, r) / np.dot(np.dot(d, A), d) w += alpha * d new_r = r - alpha * np.dot(A, d) beta = np.dot(new_r.T, new_r) / np.dot(r.T, r) d = beta * d + new_r r = new_r # 基本收敛,停止迭代 if np.linalg.norm(r) / np.linalg.norm(r0) < epsilon: break return w
相比于朴素的梯度下降法,共轭梯度法收敛迅速且稳定。不过在多项式次数增加时拟合效果会变差:在 m = 7 m=7 m=7时,其与最小二乘法对比如下:
此时,仍然可以通过正则项部分缓解(图为 m = 7 , λ = 1 m=7,\lambda=1 m=7,λ=1):
最后附上四种方法的拟合图像(基本都一样)和主函数,可以根据实验要求调整参数:
if __name__ == '__main__': warnings.simplefilter('error') dataset = get_dataset(bound = (-3, 3)) # 绘制数据集散点图 for [x, y] in dataset: plt.scatter(x, y, color = 'red') # 最小二乘法 coef1 = fit(dataset) # 岭回归 coef2 = ridge_regression(dataset) # 梯度下降法 coef3 = GD(dataset, m = 3) # 共轭梯度法 coef4 = CG(dataset) # 绘制出四种方法的曲线 draw(dataset, coef1, color = 'red', label = 'OLS') draw(dataset, coef2, color = 'black', label = 'Ridge') draw(dataset, coef3, color = 'purple', label = 'GD') draw(dataset, coef4, color = 'green', label = 'CG(lambda:0)') # 绘制标签, 显示图像 plt.legend() plt.show()