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

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

import numpy as np
import matplotlib.pyplot as plt

本实验用到的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]])

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, 1
array([0., 0., 1.])

此外，还有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 "<stdin>", 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)

np.eye

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

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

线性代数相关

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.])

生成数据

生成数据要求加入噪声（误差）。上课讲的时候举的例子就是正弦函数，我们这里也采用标准的正弦函数$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

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

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

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

最小二乘法拟合

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

解析解推导

简单回忆一下最小二乘法的原理：现在我们想用一个$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=\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⋮1​x