哈工大2022机器学习实验一：曲线拟合

这个实验的要求写的还是挺清楚的（与上学期相比），本博客采用python实现，科学计算库采用numpy，作图采用matplotlib.pyplot，为了简便在文件开头import如下：

import numpy as np
import matplotlib.pyplot as plt
1
2

本实验用到的numpy函数

一般把numpy简写为np（import numpy as np）。下面简单介绍一下实验中用到的numpy函数。下面的代码均需要在最前面加上import numpy as np。

np.array

该函数返回一个numpy.ndarray对象，可以理解为一个多维数组（本实验中仅会用到一维（可以当作列向量）和二维（矩阵））。下面用小写的$xx$表示列向量，大写的$A$表示矩阵。A.T表示$A$的转置。对ndarray的运算一般都是逐元素的。

>>> x = np.array([1,2,3])
>>> x
array([1, 2, 3])
>>> A = np.array([[2,3,4],[5,6,7]])
>>> A
array([[2, 3, 4],
       [5, 6, 7]])
>>> A.T # 转置
array([[2, 5],
       [3, 6],
       [4, 7]])
>>> A + 1
array([[3, 4, 5],
       [6, 7, 8]])
>>> A * 2
array([[ 4,  6,  8],
       [10, 12, 14]])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

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)（标准正态）
1
2
3
4
5
6
7
8
9
10

数学函数

本实验中只用到了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, 1
array([0., 0., 1.])
1
2
3

此外，还有np.log、np.exp等与python的math库相似的函数（只不过是对多维数组进行逐元素运算）。

np.dot

返回两个矩阵的乘积。与线性代数中的矩阵乘法一致。要求第一个矩阵的列等于第二个矩阵的行数。特殊地，当其中一个为一维数组时，形状会自动适配为$n\times 1$或$1\times 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 "<__array_function__ internals>", line 5, in dot
ValueError: shapes (3,3) and (1,3) not aligned: 3 (dim 1) != 1 (dim 0)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

np.eye

np.eye(n)返回一个n阶单位阵。

>>> A = np.eye(3)
>>> A
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])
1
2
3
4
5

线性代数相关

np.linalg是与线性代数有关的库。

>>> A
array([[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.])
1
2
3
4
5
6
7
8
9
10
11
12
13

生成数据

生成数据要求加入噪声（误差）。上课讲的时候举的例子就是正弦函数，我们这里也采用标准的正弦函数$y=\sin x.$（加入噪声后即为$y=\sin x+ϵ,$其中$ϵ\sim N(0,{\sigma }^2)$,由于$\sin x$的最大值为$1$,我们把误差的方差设小一点，这里设成$\frac{1}{25}$）。

'''
返回数据集，形如[[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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

产生的数据集每行为一个平面上的点。产生的数据看起来像这样：

隐隐约约能看出来是个正弦函数的形状。产生上面图像的代码如下：

dataset = get_dataset(bound = (-3, 3))
# 绘制数据集散点图
for [x, y] in dataset:
    plt.scatter(x, y, color = 'red')
plt.show()
1
2
3
4
5

最小二乘法拟合

下面我们分别用四种方法（最小二乘，正则项/岭回归，梯度下降法，共轭梯度法）以用多项式拟合上述干扰过的正弦曲线。

解析解推导

简单回忆一下最小二乘法的原理：现在我们想用一个$m$次多项式
$$f(x)=w_0+w_{1}x+w_2{x}^2+...+w_m{x}^m$$
来近似真实函数$y=\sin x.$我们的目标是最小化数据集$(x_1,y_1),(x_2,y_2),...,(x_N,y_N)$上的损失$L$（loss），这里损失函数采用平方误差：
$$L=\sum _{i=1}^N[y_i-f(x_i){]}^2$$
为了求得使均方误差最小（因此最贴合目标曲线）的参数$w_0,w_1,...,w_m,$我们需要分别求损失$L$关于$w_0,w_1,...,w_m$的导数。为了方便，我们采用线性代数的记法：
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=

_{N\times(m+1)},Y=_{N\times1},W=_{(m+1)\times1}. X=⎝ ⎛​11⋮1​x1​x2​xN​​x12​x22​xN2​​⋯⋯⋯​x1m​x2m​⋮xNm​​⎠ ⎞​N×(m+1)​,Y=⎝ ⎛​y1​y2​⋮yN​​⎠ ⎞​N×1​,W=⎝ ⎛​w0​w1​⋮wm​​⎠ ⎞​(m+1)×1​.
在这种表示方法下，有
( f ( x 1 ) f ( x 2 ) ⋮ f ( x N ) ) = X W .

(f(x1)f(x2)⋮f(xN))" role="presentation" style="position: relative;">⎛⎝⎜⎜⎜⎜f(x1)f(x2)⋮f(xN)⎞⎠⎟⎟⎟⎟$$\left(\begin{array}{c}f(x_1)\\ f(x_2)\\ ⋮\\ f(x_N)\end{array}\right)$$

= 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 .

(f(x1)−y1f(x2)−y2⋮f(xN)−yN)" role="presentation" style="position: relative;">⎛⎝⎜⎜⎜⎜f(x1)−y1f(x2)−y2⋮f(xN)−yN⎞⎠⎟⎟⎟⎟$$\left(\begin{array}{c}f(x_1)-y_1\\ f(x_2)-y_2\\ ⋮\\ f(x_N)-y_N\end{array}\right)$$

=XW-Y. ⎝ ⎛​f(x1​)−y1​f(x2​)−y2​⋮f(xN​)−yN​​⎠ ⎞​=XW−Y.
因此，损失函数
$$L=(XW-Y{)}^T(XW-Y).$$
（为了求得向量$xx=(x_1,x_2,...,x_N{)}^T$各分量的平方和，可以对$xx$作内积，即${xx}^{T}xx.$）
为了求得使$L$最小的$W$(这个$W$是一个列向量)，我们需要对$L$求偏导数，并令其为$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

∂L∂W=∂∂W[(XW−Y)T(XW−Y)]=∂∂W[(WTXT−YT)(XW−Y)]=∂∂W(WTXTXW−WTXTY−YTXW+YTY)=∂∂W(WTXTXW−2YTXW+YTY)(容易验证,WTXTY=YTXW,因而可以将其合并)=2XTXW−2XTY" role="presentation">∂L∂W=∂∂W[(XW−Y)T(XW−Y)]=∂∂W[(WTXT−YT)(XW−Y)]=∂∂W(WTXTXW−WTXTY−YTXW+YTY)=∂∂W(WTXTXW−2YTXW+YTY)(容易验证,WTXTY=YTXW,因而可以将其合并)=2XTXW−2XTY$$\begin{array}{rl}\frac{\mathrm{\partial }L}{\mathrm{\partial }W}& =\frac{\mathrm{\partial }}{\mathrm{\partial }W}[(XW-Y{)}^T(XW-Y)]\\ & =\frac{\mathrm{\partial }}{\mathrm{\partial }W}[(W^T{X}^T-Y^T)(XW-Y)]\\ & =\frac{\mathrm{\partial }}{\mathrm{\partial }W}(W^T{X}^{T}XW-W^T{X}^{T}Y-Y^{T}XW+Y^{T}Y)\\ & =\frac{\mathrm{\partial }}{\mathrm{\partial }W}(W^T{X}^{T}XW-2Y^{T}XW+Y^{T}Y)(容易验证,W^T{X}^{T}Y=Y^{T}XW,因而可以将其合并)\\ & =2X^{T}XW-2X^{T}Y\end{array}$$

∂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$和$Y^{T}XW$都是数（或者说$1\times 1$矩阵），二者互为转置，因此值相同，可以合并成一项。
（2）从第4行到第5行的矩阵求导，第一项$\frac{\partial }{\partial W}(W^T(X^{T}X)W)$是一个关于$W$的二次型，其导数就是$2X^{T}XW.$
（3）对于一次项$-2Y^{T}XW$的求导，如果按照实数域的求导应该得到$-2Y^{T}X.$但检查一下发现矩阵的型对不上，需要做一下转置，变为$-2X^{T}Y.$

矩阵求导线性代数课上也没有系统教过，只对这里出现的做一下说明。（多了我也不会 ）
令偏导数为0，得到
$$X^{T}XW=Y^{T}X,$$
左乘$(X^{T}X{)}^{-1}$（$X^{T}X$的可逆性见下方的补充说明），得到
$$W=(X^{T}X{)}^{-1}{X}^{T}Y.$$
这就是我们想求的$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)
1
2
3
4
5
6
7
8
9
10

稍微解释一下代码：第一行即生成上面约定的$X$矩阵，dataset[:,0]即数据集第0列$(x_1,x_2,...,x_N{)}^T$；第二行即$Y$矩阵；第三行返回上面的解析解。（如果不熟悉python语法或者numpy库还是挺不友好的）

简单地验证一下我们已经完成的函数的结果：为此，我们先写一个draw函数，用于把求得的$W$对应的多项式$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)
1
2
3
4
5
6
7
8
9
10
11
12

然后是主函数：

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()
1
2
3
4
5
6
7
8
9
10
11
12

可以看到5次多项式拟合的效果还是比较不错的（数据集每次随机生成，所以跟第一幅图不一样）。

截至这部分全部的代码，后面同名函数不再给出说明：

import numpy as np
import 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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

补充说明

上面有一块不太严谨：对于一个矩阵$X$而言，$X^{T}X$不一定可逆。然而在本实验中，可以证明其为可逆矩阵。由于这门课不是线性代数课，我们就不费太多篇幅介绍这个了，仅作简单提示：
（1）$X$是一个$N\times (m+1)$的矩阵。其中数据数$N$远大于多项式次数$m$，有$N>m+1;$
（2）为了说明$X^{T}X$可逆，需要说明$(X^{T}X{)}_{(m+1)\times (m+1)}$满秩，即$R(X^{T}X)=m+1;$
（3）在线性代数中，我们证明过$R(X)=R(X^T)=R(X^{T}X)=R(X{X}^T);$
（4）$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')
1
2
3
4
5
6
7
8
9

过拟合在$m$较大时尤为严重（上面图像为$m=3$时）。当多项式次数升高时，为了尽可能贴近所给数据集，计算出来的系数的数量级将会越来越大，在未见样本上的表现也就越差。如上图，可以看到拟合在前50个点（大约在横坐标$[-3,0]$处）表现很好；而在测试集上表现就很差（$[0,3]$处）。为了防止过拟合，可以引入正则化项。此时损失函数$L$变为
$$L=(XW-Y{)}^T(XW-Y)+\lambda ||W|{|}_2^2$$
其中$||\cdot |{|}_2^2$表示$L_2$范数的平方，在这里即$W^{T}W;\lambda$为正则化系数。该式子也称岭回归（Ridge Regression）。它的思想是兼顾损失函数与所得参数$W$的模长（在$L_2$范数时），防止$W$内的参数过大。

举个例子(数是随便编的)：当正则化系数为$1$，若方案1在数据集上的平方误差为$0.5,$此时$W=(100,-200,300,150{)}^T$；方案2在数据集上的平方误差为$10,$此时$W=(1,-3,2,1)$,那我们选择方案2的$W.$正则化系数$\lambda$刻画了这种对于$W$模长的重视程度：$\lambda$越大，说明$W$的模长升高带来的惩罚也就越大。当$\lambda =0,$岭回归即变为普通的最小二乘法。与岭回归相似的还有LASSO，就是将正则化项换为$L_1$范数。

重复上面的推导，我们可以得出解析解为
$$W=(X^{T}X+\lambda E_{m+1}{)}^{-1}{X}^{T}Y.$$
其中$E_{m+1}$为$m+1$阶单位阵。容易得到$(X^{T}X+\lambda E_{m+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)
1
2
3
4
5
6
7
8
9
10
11

两种方法的对比如下：

对比可以看出，岭回归显著减轻了过拟合（此时为$m=3,\lambda =0.3$）。

梯度下降法

梯度下降法并不是求解该问题的最好方法，很容易就无法收敛。先简单介绍梯度下降法的基本思想：若我们想求取复杂函数$f(x)$的最小值（最值点）（这个$x$可能是向量等），即
$$x_{min}=\underset{x}{\mathrm{argmin}}f(x)$$
梯度下降法重复如下操作：
（0）（随机）初始化$x_0(t=0)$；
（1）设$f(x)$在$x_t$处的梯度（当$x$为一维时，即导数）$\nabla f(x_t)$；
（2）$x_{t+1}=x_t-\eta \nabla f(x_t)$
（3）若$x_{t+1}$与$x_t$相差不大（达到预先设定的范围）或迭代次数达到预设上限，停止算法；否则重复（1）（2）.

其中$\eta$为学习率，它决定了梯度下降的步长。
下面是一个用梯度下降法求取$y=x^2$的最小值点的示例程序：

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return x ** 2

def draw():
    x = np.linspace(-3, 3)
    y = f(x)
    plt.plot(x, y, c = 'red')

cnt = 0
# 初始化 x
x = np.random.rand(1) * 3
learning_rate = 0.05

while 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 += 1

draw()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

上图标明了$x$随着迭代的演进，可以看到$x$不断沿着正半轴向零点靠近。需要注意的是，学习率不能过大（虽然在上面的程序中，学习率设置得有点小了），需要手动进行尝试调整，否则容易想象，$x$在正负半轴来回震荡，难以收敛。

在最小二乘法中，我们需要优化的函数是损失函数
$$L=(XW-Y{)}^T(XW-Y).$$
下面我们用梯度下降法求解该问题。在上面的推导中，
∂ L ∂ W = 2 X T X W − 2 X T Y ,

, ∂W∂L​=2XTXW−2XTY​,
于是我们每次在迭代中对$W$减去该梯度，直到参数$W$收敛。不过经过实验，平方误差会使得梯度过大，过程无法收敛，因此采用均方误差(MSE)替换之，就是给原来的式子除以$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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

这时如果$m$设置得稍微大一点（比如4），在迭代过程中梯度就会溢出，使参数无法收敛。在收敛时，拟合效果还算可以：

共轭梯度法

共轭梯度法（Conjugate Gradients）可以用来求解形如$Axx=bb$的方程组，或最小化二次型$f(xx)=\frac{1}{2}{xx}^{T}Axx-{bb}^{T}xx+c.$（可以证明对于正定的$A$,二者等价）其中$A$为正定矩阵。在本问题中，我们要求解$$X^{T}XW=Y^{T}X,$$
就有$A_{(m+1)\times (m+1)}=X^{T}X,bb=Y^{T}X.$若我们想加一个正则项，就变成求解
$$(X^{T}X+\lambda E)W=Y^{T}X.$$
首先说明一点：$X^{T}X$不一定是正定的但一定是半正定的（证明见此）。但是在实验中我们基本不用担心这个问题，因为$X^{T}X$有极大可能是正定的，我们只在代码中加一个断言（assert），不多关注这个条件。
共轭梯度法的思想来龙去脉和证明过程比较长，可以参考这个系列，这里只给出算法步骤（在上面链接的第三篇开头）：

（0）初始化$x_{(0)};$
（1）初始化$d_{(0)}=r_{(0)}=b-A{x}_{(0)};$
（2）令$${\alpha }_{(i)}=\frac{r_{(i)}^T{r}_{(i)}}{d_{(i)}^{T}A{d}_{(i)}};$$
（3）迭代$x_{(i+1)}=x_{(i)}+{\alpha }_{(i)}{d}_{(i)};$
（4）令$r_{(i+1)}=r_{(i)}-{\alpha }_{(i)}A{d}_{(i)};$
（5）令$${\beta }_{(i+1)}=\frac{r_{(i+1)}^T{r}_{(i+1)}}{r_{(i)}^T{r}_{(i)}},d_{(i+1)}=r_{(i+1)}+{\beta }_{(i+1)}{d}_{(i)}.$$
（6）当$\frac{||r_{(i)}||}{||r_{(0)}||}<ϵ$时，停止算法；否则继续从（2）开始迭代。$ϵ$为预先设定好的很小的值，我这里取的是$10^{-8}.$
下面我们按照这个过程实现代码：

'''
共轭梯度法(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-8

    # 初始化参数
    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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

相比于朴素的梯度下降法，共轭梯度法收敛迅速且稳定。

最后附上四种方法的拟合图像（基本都一样）和主函数，可以根据实验要求调整参数：

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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27