源代码请参考：实验一 GitHub 仓库
运行效果请参考：主程序

哈尔滨工业大学计算学部

实验报告


《机器学习》
实验一：多项式拟合正弦函数

学号：1183710109 姓名：郭茁宁

一、实验目的

掌握最小二乘法求解（无惩罚项的损失函数），掌握增加惩罚项（2 范数）的损失函数优化，梯度下降法、共轭梯度法，理解过拟合、客服过拟合的方法（如增加惩罚项、增加样本）。

二、实验要求及实验环境

实验要求

1. 生成数据，加入噪声；
2. 用高阶多项式函数拟合曲线；
3. 用解析解求解两种 loss 的最优解（无正则项和有正则项）
4. 优化方法求解最优解（梯度下降，共轭梯度）；
5. 用你得到的实验数据，解释过拟合。
6. 用不同数据量，不同超参数，不同的多项式阶数，比较实验效果。
7. 语言不限，可以用 matlab，python。求解解析解时可以利用现成的矩阵求逆。梯度下降，共轭梯度要求自己求梯度，迭代优化自己写。不许用现成的平台，例如 pytorch，tensorflow 的自动微分工具。

实验环境

• OS: Win 10
• Python 3.6.8

三、算法原理和设计

1、生成数据算法

主要是利用 s i n ( 2 π x ) sin(2\pi x) sin(2πx)函数产生样本，其中 x x x均匀分布在 [ 0 , 1 ] [0, 1] [0,1]之间，对于每一个目标值 t = s i n ( 2 π x ) t=sin(2\pi x) t=sin(2πx)增加一个 0 0 0均值，方差为 0.25 0.25 0.25的高斯噪声。

def get_data(
    x_range: (float, float) = (0, 1),
    sample_num: int = 10,
    base_func=lambda x: np.sin(2 * np.pi * x),
    noise_scale=0.25,
) -> "pd.DataFrame":
    X = np.linspace(x_range[0], x_range[1], num=sample_num)
    Y = base_func(X) + np.random.normal(scale=noise_scale, size=X.shape)
    data = pd.DataFrame(data=np.dstack((X, Y))[0], columns=["X", "Y"])
    return data

2、利用高阶多项式函数拟合曲线(不带惩罚项)

利用训练集合，对于每个新的 x ^ \hat{x} x^，预测目标值 t ^ \hat{t} t^。采用多项式函数进行学习，即利用式 ( 1 ) (1) (1)来确定参数 w w w，假设阶数 m m m已知。

y ( x , w ) = w 0 + w 1 x + ⋯ + w m x m = ∑ i = 0 m w i x i (1) y(x, w) = w_0 + w_1x + \dots + w_mx^m = \sum_{i = 0}^{m}w_ix^i \tag{1} y(x,w)=w0​+w1​x+⋯+wm​xm=i=0∑m​wi​xi(1)

采用最小二乘法，即建立误差函数来测量每个样本点目标值 t t t与预测函数 y ( x , w ) y(x, w) y(x,w)之间的误差，误差函数即式 ( 2 ) (2) (2)

E ( w ) = 1 2 ∑ i = 1 N { y ( x i , w ) − t i } 2 (2) E(\bold{w}) = \frac{1}{2} \sum_{i = 1}^{N} \{y(x_i, \bold{w}) - t_i\}^2 \tag{2} E(w)=21​i=1∑N​{y(xi​,w)−ti​}2(2)

将上式写成矩阵形式如式 ( 3 ) (3) (3)

E ( w ) = 1 2 ( X w − T ) ′ ( X w − T ) (3) E(\bold{w}) = \frac{1}{2} (\bold{Xw} - \bold{T})'(\bold{Xw} - \bold{T})\tag{3} E(w)=21​(Xw−T)′(Xw−T)(3)

其中 X = [ 1 x 1 ⋯ x 1 m 1 x 2 ⋯ x 2 m ⋮ ⋮ ⋱ ⋮ 1 x N ⋯ x N m ] , w = [ w 0 w 1 ⋮ w m ] , T = [ t 1 t 2 ⋮ t N ] \bold{X} = \left[ \begin{matrix} 1 & x_1 & \cdots & x_1^m \\ 1 & x_2 & \cdots & x_2^m \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_N & \cdots & x_N^m \\ \end{matrix} \right], \bold{w} = \left[ \begin{matrix} w_0 \\ w_1 \\ \vdots \\ w_m \end{matrix} \right], \bold{T} = \left[ \begin{matrix} t_1 \\ t_2 \\ \vdots \\ t_N \end{matrix} \right] X=⎣⎢⎢⎢⎡​11⋮1​x1​x2​⋮xN​​⋯⋯⋱⋯​x1m​x2m​⋮xNm​​⎦⎥⎥⎥⎤​,w=⎣⎢⎢⎢⎡​w0​w1​⋮wm​​⎦⎥⎥⎥⎤​,T=⎣⎢⎢⎢⎡​t1​t2​⋮tN​​⎦⎥⎥⎥⎤​
通过将上式求导我们可以得到式 ( 4 ) (4) (4)

∂ E ∂ w = X ′ X w − X ′ T (4) \cfrac{\partial E}{\partial \bold{w}} = \bold{X'Xw} - \bold{X'T} \tag{4} ∂w∂E​=X′Xw−X′T(4)

令 ∂ E ∂ w = 0 \cfrac{\partial E}{\partial \bold{w}}=0 ∂w∂E​=0 我们有式 ( 5 ) (5) (5)即为 w ∗ \bold{w^*} w∗

w ∗ = ( X ′ X ) − 1 X ′ T (5) \bold{w^*} = (\bold{X'X})^{-1}\bold{X'T} \tag{5} w∗=(X′X)−1X′T(5)

由 ( 5 ) (5) (5)实现

def get_params(x_matrix, t_vec) -> [float]:
    return np.linalg.pinv(x_matrix.T @ x_matrix) @ x_matrix.T @ t_vec

3、带惩罚项的多项式函数拟合曲线

在不带惩罚项的多项式拟合曲线时，在参数多时 w ∗ \bold{w}^* w∗具有较大的绝对值，本质就是发生了过拟合。对于这种过拟合，我们可以通过在优化目标函数式 ( 3 ) (3) (3)中增加 w \bold{w} w的惩罚项，因此我们得到了式 ( 6 ) (6) (6)。

E ~ ( w ) = 1 2 ∑ i = 1 N { y ( x i , w ) − t i } 2 + λ 2 ∣ ∣ w ∣ ∣ 2 (6) \widetilde{E}(\bold{w}) = \frac{1}{2} \sum_{i=1}^{N} \{y(x_i, \bold{w}) - t_i\}^2 + \cfrac{\lambda}{2}|| \bold{w} || ^ 2 \tag{6} E (w)=21​i=1∑N​{y(xi​,w)−ti​}2+2λ​∣∣w∣∣2(6)

同样我们可以将式 ( 6 ) (6) (6)写成矩阵形式， 我们得到式 ( 7 ) (7) (7)

E ~ ( w ) = 1 2 [ ( X w − T ) ′ ( X w − T ) + λ w ′ w ] (7) \widetilde{E}(\bold{w}) = \frac{1}{2}[(\bold{Xw} - \bold{T})'(\bold{Xw} - \bold{T}) + \lambda \bold{w'w}]\tag{7} E (w)=21​[(Xw−T)′(Xw−T)+λw′w](7)

对式 ( 7 ) (7) (7)求导我们得到式 ( 8 ) (8) (8)

∂ E ~ ∂ w = X ′ X w − X ′ T + λ w (8) \cfrac{\partial \widetilde{E}}{\partial \bold{w}} = \bold{X'Xw} - \bold{X'T} + \lambda\bold{w} \tag{8} ∂w∂E ​=X′Xw−X′T+λw(8)

令 ∂ E ~ ∂ w = 0 \cfrac{\partial \widetilde{E}}{\partial \bold{w}} = 0 ∂w∂E ​=0 我们得到 w ∗ \bold{w^*} w∗即式 ( 9 ) (9) (9)，其中 I \bold{I} I为单位阵。

w ∗ = ( X ′ X + λ I ) − 1 X ′ T (9) \bold{w^*} = (\bold{X'X} + \lambda\bold{I})^{-1}\bold{X'T}\tag{9} w∗=(X′X+λI)−1X′T(9)

由 ( 9 ) (9) (9)实现

def get_params_with_penalty(x_matrix, t_vec, lambda_penalty):
    return (
        np.linalg.pinv(
            x_matrix.T @ x_matrix + lambda_penalty * np.identity(x_matrix.shape[1])
        )
        @ x_matrix.T
        @ t_vec
    )

4、梯度下降法求解最优解

对于 f ( x ) f(\bold{x}) f(x)如果在 x i \bold{x_i} xi​点可微且有定义，我们知道顺着梯度 ∇ f ( x i ) \nabla f(\bold{x_i}) ∇f(xi​)为增长最快的方向，因此梯度的反方向 − ∇ f ( x i ) -\nabla f(\bold{x_i}) −∇f(xi​) 即为下降最快的方向。因而如果有式 ( 10 ) (10) (10)对于 α > 0 \alpha > 0 α>0 成立,

x i + 1 = x i − α ∇ f ( x i ) (10) \bold{x_{i+1}}= \bold{x_i} - \alpha \nabla f(\bold{x_i}) \tag{10} xi+1​=xi​−α∇f(xi​)(10)

那么对于序列 x 0 , x 1 , … \bold{x_0}, \bold{x_1}, \dots x0​,x1​,… 我们有 f ( x 0 ) ≥ f ( x 1 ) ≥ … f(\bold{x_0}) \ge f(\bold{x_1}) \ge \dots f(x0​)≥f(x1​)≥…

因此，如果顺利我们可以得到一个 f ( x n ) f(\bold{x_n}) f(xn​) 收敛到期望的最小值，对于此次实验很大可能性可以收敛到最小值。

进而我们实现算法如下，其中 δ \delta δ为精度要求，通常可以设置为 δ = 1 × 1 0 − 6 \delta = 1\times10^{-6} δ=1×10−6：

def gradient_descent_fit(x_matrix, t_vec, lambda_penalty, w_vec_0, learning_rate=0.1, delta=1e-6):
    loss_0 = calc_loss(x_matrix, t_vec, lambda_penalty, w_vec_0)
    k = 0
    w = w_vec_0
    while True:
        w_ = w - learning_rate * calc_derivative(x_matrix, t_vec, lambda_penalty, w)
        loss = calc_loss(x_matrix, t_vec, lambda_penalty, w_)
        if np.abs(loss - loss_0) < delta:
            break
        else:
            k += 1
            if loss > loss_0:
                learning_rate *= 0.5
            loss_0 = loss
            w = w_
    return k, w

5、共轭梯度法求解最优解

共轭梯度法解决的主要是形如 A x = b \bold{Ax} = \bold{b} Ax=b的线性方程组解的问题，其中 A \bold{A} A必须是对称的、正定的。大概来说，共轭梯度下降就是在解空间的每一个维度分别取求解最优解的，每一维单独去做的时候不会影响到其他维，这与梯度下降方法，每次都选泽梯度的反方向去迭代，梯度下降不能保障每次在每个维度上都是靠近最优解的，这就是共轭梯度优于梯度下降的原因。

对于第 k 步的残差 r k = b − A x k \bold{r_k = b - Ax_k} rk​=b−Axk​，我们根据残差去构造下一步的搜索方向 p k \bold{p_k} pk​，初始时我们令 p 0 = r 0 \bold{p_0 = r_0} p0​=r0​。然后利用 G r a m − S c h m i d t Gram-Schmidt Gram−Schmidt方法依次构造互相共轭的搜素方向 p k \bold{p_k} pk​，具体构造的时候需要先得到第 k+1 步的残差，即 r k + 1 = r k − α k A p k \bold{r_{k+1} = r_k - \alpha_kAp_k} rk+1​=rk​−αk​Apk​，其中 α k \alpha_k αk​如后面的式 ( 11 ) (11) (11)。

根据第 k+1 步的残差构造下一步的搜索方向 p k + 1 = r k + 1 + β k + 1 p k \bold{p_{k+1} = r_{k+1} + \beta_{k+1}p_k} pk+1​=rk+1​+βk+1​pk​，其中 β k + 1 = r k + 1 T r k + 1 r k T r k \beta_{k+1} = \bold{\cfrac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k}} βk+1​=rkT​rk​rk+1T​rk+1​​。

然后可以得到 x k + 1 = x k + α k p k \bold{x_{k+1} = x_k + \alpha_kp_k} xk+1​=xk​+αk​pk​，其中 α k = p k T ( b − A x k ) p k T A p k = p k T r k p k T A p k \alpha_k = \cfrac{\bold{p_k}^T(\bold{b} - \bold{Ax_k})}{\bold{p_k}^T\bold{Ap_k}} = \cfrac{\bold{p_k}^T\bold{r_k}}{{\bold{p_k}^T\bold{Ap_k}}} αk​=pk​TApk​pk​T(b−Axk​)​=pk​TApk​pk​Trk​​

对于第 k 步的残差 r k = b − A x \bold{r_k} = \bold{b} - \bold{Ax} rk​=b−Ax， r k \bold{r_k} rk​为 x = x k \bold{x} = \bold{x_k} x=xk​时的梯度反方向。由于我们仍然需要保证 p k \bold{p_k} pk​ 彼此共轭。因此我们通过当前的残差和之前所有的搜索方向来构建 p k \bold{p_k} pk​，得到式 ( 11 ) (11) (11)

p k = r k − ∑ i < k p i T A r k p i T A p i p i (11) \bold{p_k} = \bold{r_k} - \sum_{i < k}\cfrac{\bold{p_i}^T\bold{Ar_k}}{\bold{p_i}^T\bold{Ap_i}}\bold{p_i}\tag{11} pk​=rk​−i<k∑​pi​TApi​pi​TArk​​pi​(11)

进而通过当前的搜索方向 p k \bold{p_k} pk​得到下一步优化解 x k + 1 = x k + α k p k \bold{x_{k+1}} = \bold{x_k} + \alpha_k\bold{p_k} xk+1​=xk​+αk​pk​，其中 α k = p k T ( b − A x k ) p k T A p k = p k T r k p k T A p k \alpha_k = \cfrac{\bold{p_k}^T(\bold{b} - \bold{Ax_k})}{\bold{p_k}^T\bold{Ap_k}} = \cfrac{\bold{p_k}^T\bold{r_k}}{{\bold{p_k}^T\bold{Ap_k}}} αk​=pk​TApk​pk​T(b−Axk​)​=pk​TApk​pk​Trk​​

求解的方法就是我们先猜一个解 x 0 \bold{x_0} x0​，然后取梯度的反方向 p 0 = b − A x \bold{p_0} = \bold{b} - \bold{Ax} p0​=b−Ax，在 n 维空间的基中 p 0 \bold{p_0} p0​要与其与的基共轭并且为初始残差。

对于共轭梯度下降，算法实现如下，初始时取 w 0 = [ 0 0 ⋮ 0 ] \bold{w_0} = \left[ \begin{matrix} 0 \\ 0 \\ \vdots \\ 0 \end{matrix} \right] w0​=⎣⎢⎢⎢⎡​00⋮0​⎦⎥⎥⎥⎤​

由上文 ( 8 ) ( 9 ) (8)(9) (8)(9)得

( X ′ X + λ I ) w ∗ = X ′ T (12) (\bold{X'X} + \lambda\bold{I})\bold{w^*} = \bold{X'T}\tag{12} (X′X+λI)w∗=X′T(12)

令

{ A = X ′ X + λ I b = X ′ T (13) \left\{ \begin{array}{lr} \bold{A} = \bold{X'X} + \lambda\bold{I} & \\ \bold{b} = \bold{X'T} \end{array} \right.\tag{13} {A=X′X+λIb=X′T​(13)

def switch_deri_func_for_conjugate_gradient(x_matrix, t_vec, lambda_penalty, w_vec):
    A = x_matrix.T @ x_matrix - lambda_penalty * np.identity(len(x_matrix.T))
    x = w_vec
    b = x_matrix.T @ t_vec
    return A, x, b

通过共轭梯度下降法求解 A x = b Ax=b Ax=b

def conjugate_gradient_fit(A, x_0, b, delta=1e-6):
    x = x_0
    r_0 = b - A @ x
    p = r_0
    k = 0
    while True:
        alpha = (r_0.T @ r_0) / (p.T @ A @ p)
        x = x + alpha * p
        r = r_0 - alpha * A @ p
        if r_0.T @ r_0 < delta:
            break
        beta = (r.T @ r) / (r_0.T @ r_0)
        p = r + beta * p
        r_0 = r
        k += 1
    return k, x

解得 x x x，并记录迭代次数 k k k

四、实验结果分析

1、不带惩罚项的解析解

固定训练样本的大小为 15，分别使用不同多项式阶数，测试：

我们可以看到在固定训练样本的大小之后，在多项式阶数为 3 时的拟合效果已经很好。继续提高多项式的阶数，尤其在阶数为 14 的时候曲线“完美的”经过了所有的节点，这种剧烈的震荡并没有很好的拟合真实的背后的函数 s i n ( 2 π x ) sin(2\pi x) sin(2πx)，反而将所有噪声均很好的拟合，即表现出来一种过拟合的情况。

其出现过拟合的本质原因是，在阶数过大的情况下，模型的复杂度和拟合的能力都增强，因此可以通过过大或者过小的系数来实现震荡以拟合所有的数据点，以至于甚至拟合了所有的噪声。在这里由于我们的数据样本大小只有 15，所以在阶数为 14 的时候，其对应的系数向量 w \bold{w} w恰好有唯一解，因此可以穿过所有的样本点。

对于过拟合我们可以通过增加样本的数据或者通过增加惩罚项的方式来解决。增加数据集样本数量，使其超过参数向量的大小，就会在一定程度上解决过拟合问题。

固定多项式阶数，使用不同数量的样本数据，测试

可以看到在固定多项式阶数为 7 的情况下，随着样本数量逐渐增加，过拟合的现象有所解决。特别是对比左上图与右下图的差距，可以看到样本数量对于过拟合问题是有影响的。

2、带惩罚项的解析解

首先根据式 ( 6 ) (6) (6)我们需要确定最佳的超参数 λ \lambda λ，因此我们通过根均方(RMS)误差来确定

观察上面的四张图， 我们可以发现对于超参数 λ \lambda λ的选择，在 ( e − 50 , e − 30 ) (e^{-50}, e^{-30}) (e−50,e−30)左右保持一种相对稳定的错误率；但是在 ( e − 30 , e − 5 ) (e^{-30}, e^{-5}) (e−30,e−5)错误率有一个明显的下降，所以下面在下面的完整 100 次实验中我们可以看到最佳参数的分布区间也大都在这个范围内；在大于 e − 5 e^{-5} e−5的区间内，错误率有一个急剧的升高。

对于不同训练集，超参数 λ \lambda λ的选择不同。
因此每生成新训练集，则通过带惩罚项的解析解计算出当前的 b e s t _ l a m b d a best\_lambda best_lambda，用于带惩罚项的解析解、梯度下降和共轭梯度下降。

一般来说，最佳的超参数范围在 ( e − 10 , e − 6 ) (e^{-10}, e^{-6}) (e−10,e−6)之间。

比较是否带惩罚项的拟合曲线

3、优化解

实验利用梯度下降(Gradient descent)和共轭梯度法(Conjugate gradient)两种方法来求优化解。由于该问题是有解析解存在，并且在之前的实验报告部分已经列出，所以在此处直接利用上面的解析解进行对比。此处使用的学习率固定为 l e a r n i n g _ r a t e = 0.01 learning\_rate = 0.01 learning_rate=0.01，停止的精度要求为 1 × 1 0 − 6 1 \times 10 ^ {-6} 1×10−6。

阶数	训练集规模	梯度下降迭代次数	共轭梯度下降迭代次数
3	10	110939	4
3	20	92551	4
3	50	56037	4
5	10	38038	4
5	20	22867	5
5	50	118788	5
8	10	81387	5
8	20	71743	7
8	50	47445	8

首先在固定多项式阶数的情况下，随着训练样本的增加，梯度下降的迭代次数均有所下降，但是对于共轭梯度迭代次数变化不大。

其次在固定训练样本的情况下，梯度下降迭代次数的变化，对于 3 阶的情况下多于 8 阶的情况对于共轭梯度的而言，迭代次数仍然较少。

总的来说，对于梯度下降法，迭代次数在 10000 次以上；而对于共轭梯度下降，则需要的迭代次数均不超过 10 次(即小于解空间的维度 M)。

下面是不同阶数和不同训练集规模上，两种方法的比较

梯度下降和共轭梯度下降两种方法拟合结果相似，但共轭梯度下降收敛速度却远优于梯度下降、梯度下降稳定程度略优于共轭梯度下降。

4、四种拟合对比

3 阶，训练集规模 10

w w w	Analytical_Without_Penalty	Analytical_With_Penalty	Gradient_Descent	Conjugate_Gradient
0	0.086879	0.086879	0.221466	0.086879
1	10.030641	10.030641	8.037745	10.030641
2	-29.339163	-29.339163	-24.321837	-29.339163
3	18.911480	18.911480	15.658977	18.911480

5 阶，训练集规模 15

w w w	Analytical_Without_Penalty	Analytical_With_Penalty	Gradient_Descent	Conjugate_Gradient
0	0.179098	0.179098	0.252970	0.127280
1	6.186679	6.186679	6.120921	8.802575
2	-4.197192	-4.197192	-14.547026	-24.093028
3	-47.562082	-47.562082	-0.720464	6.293567
4	73.581708	73.581708	6.179084	13.416870
5	-28.074476	-28.074476	2.945509	-4.398728

8 阶，训练集规模 15

w w w	Analytical_Without_Penalty	Analytical_With_Penalty	Gradient_Descent	Conjugate_Gradient
0	0.106875	-0.090665	-0.015363	-0.394036
1	-3.384695	9.291901	7.889612	18.583773
2	37.377871	-19.656246	-16.016235	-65.870523
3	451.762357	-6.583293	-3.968438	54.941302
4	-3775.172514	16.087370	5.885670	25.695923
5	10453.330011	10.996164	7.857581	-17.630151
6	-13828.893913	-3.553833	4.750929	-30.267690
7	8931.020178	-9.096546	-0.458928	-12.072342
8	-2266.131841	2.553330	-5.970533	27.074041

15 阶，训练集规模 20

w w w	Analytical_Without_Penalty	Analytical_With_Penalty	Gradient_Descent	Conjugate_Gradient
0	0.333171	0.056738	0.132317	-0.336080
1	-43.880559	7.449353	6.433572	14.722545
2	1033.334292	-13.398546	-12.067511	-36.803428
3	-8505.117304	-8.201093	-5.256489	3.989214
4	34811.021660	3.048536	2.141618	17.244891
5	-75466.068594	8.612557	5.268958	13.535997
6	75267.404153	8.529704	5.150508	5.089805
7	786.304140	5.220844	3.410453	-2.234687
8	-53324.724219	0.884139	1.230193	-6.603244
9	-2109.290031	-2.981579	-0.689982	-8.038983
10	42045.345941	-5.509223	-1.994826	-7.225030
11	14623.791181	-6.291021	-2.546300	-4.952217
12	-33191.453129	-5.219866	-2.332749	-1.904370
13	-22428.414650	-2.369747	-1.410270	1.399483
14	38870.202234	2.086673	0.133219	4.598064
15	-12368.829171	7.931668	2.199092	7.459104

五、结论

• 增加训练样本的数据可以有效的解决过拟合的问题。
• 对于训练样本限制较多的问题，通过增加惩罚项仍然可以有效解决过拟合问题。
• 对于梯度下降法和共轭梯度法而言，梯度下降收敛速度较慢，共轭梯度法的收敛速度快；且二者相对于解析解而言，共轭梯度法的拟合效果解析解的效果更好。
• 相比之下，共轭梯度下降能够有效的解决梯度下降法迭代次数多，和复杂度高的有效方法。

六、参考文献

• Pattern Recognition and Machine Learning.
• Gradient descent wiki
• Conjugate gradient method wiki
• Shewchuk J R. An introduction to the conjugate gradient method without the agonizing pain[J]. 1994.

七、附录（代码）

源代码：实验一 GitHub 仓库

main.ipynb

Jupyter Notebook 运行效果：主程序

import numpy as np
import pandas as pd
import seaborn as sns
from matplotlib import pyplot as plt

from DataGenerator import *
from AnalyticalSolution import *
from Switcher import *
from Calculator import *
from GradientDescent import *
from ConjugateGradient import *

%matplotlib inline
sns.set(palette="Set2")

# 超参数
TRAIN_NUM = 10 # 训练集规模
VALIDATION_NUM = 100 # 验证集规模
TEST_NUM = 1000 # 测试集规模
ORDER = 7 # 阶数
X_LEFT = 0 # 左界限
X_RIGHT = 1 # 右界限
NOISE_SCALE = 0.25 # 噪音标准差
LEARNING_RATE = 0.01 # 梯度下降学习率
DELTA = 1e-6 # 优化残差界限
W_SOLUTION = pd.DataFrame() # 不同方法的w解
LAMBDA = np.power(np.e, -7)

def base_func(x):
    """原始函数
    """
    return np.sin(2 * np.pi * x)

# 训练集
train_data = get_data(
    x_range=(X_LEFT, X_RIGHT),
    sample_num=TRAIN_NUM,
    base_func=base_func,
    noise_scale=NOISE_SCALE,
)
X_TRAIN = train_data["X"]
Y_TRAIN = train_data["Y"]

# 验证集
X_VALIDATION = np.linspace(X_LEFT, X_RIGHT, VALIDATION_NUM)
Y_VALIDATION = base_func(X_VALIDATION)

# 测试集
X_TEST = np.linspace(X_LEFT, X_RIGHT, TEST_NUM)
Y_TEST = base_func(X_TEST)

train_data

	X	Y
0	0.000000	0.087444
1	0.111111	0.667822
2	0.222222	0.826025
3	0.333333	0.563550
4	0.444444	0.548170
5	0.555556	-0.388687
6	0.666667	-0.803972
7	0.777778	-1.049368
8	0.888889	-0.648516
9	1.000000	0.234352

可视化拟合结果

def show_train_data(x_train, y_train):
    plt.scatter(x=x_train, y=y_train, color="white", edgecolors="darkblue")

def show_test_data(x_test, y_test):
    plt.plot(x_test, y_test, linewidth=2, color="gray", linestyle="-.")

def show_order_and_train_num():
    plt.title("ORDER = " + str(ORDER) + ", TRAIN_NUM = " + str(TRAIN_NUM))

def show_comparation(x, y1, label1, y2, label2):
    """可视化两个函数
    """
    d = pd.DataFrame(
        np.concatenate(
            (
                np.transpose([X_TEST, y1, [label1] * TEST_NUM]),
                np.transpose([X_TEST, y2, [label2] * TEST_NUM]),
            )
        ),
        columns=["X", "Y", "type"],
    )
    d[["X", "Y"]] = d[["X", "Y"]].astype("float")
    sns.lineplot(x="X", y="Y", data=d, hue="type")
    show_order_and_train_num()

解析解

def get_y_pred_by_analytical(x_train, y_train, x, with_penalty=False, lambda_penalty=None):
    assert not with_penalty or lambda_penalty is not None  # 有惩罚项时，lambda不为空
    x_vec, t_vec = x_train, y_train
    w_vec = (
        get_params_with_penalty(
            x_matrix=get_x_matrix(x_vec, order=ORDER),
            t_vec=t_vec,
            lambda_penalty=lambda_penalty,
        )
        if with_penalty
        else get_params(x_matrix=get_x_matrix(x_vec, order=ORDER), t_vec=t_vec)
    )
    return np.dot(get_x_matrix(x, order=ORDER), w_vec), w_vec

不带惩罚项解析解

y_pred_without_penalty, w_wec_without_penalty = get_y_pred_by_analytical(
    x_train=X_TRAIN, y_train=Y_TRAIN, x=X_TEST, with_penalty=False
)
W_SOLUTION["Analytical_Without_Penalty"] = w_wec_without_penalty

show_comparation(x=X_TEST, y1=Y_TEST, label1="True", y2=y_pred_without_penalty, label2="Pred")
show_train_data(x_train=X_TRAIN, y_train=Y_TRAIN)

带惩罚项解析解

def show_lambda_error(error_ln_lambda):
    data = pd.DataFrame(error_ln_lambda, columns=["$ln{\lambda}$", "$E_{rms}$", "type"])
    sns.lineplot(x="$ln{\lambda}$", y="$E_{rms}$", data=data, hue="type")
    idxmin = error_ln_lambda[data[data["type"]=="VALIDATION"]["$E_{rms}$"].idxmin()]
    plt.title(label=("Min: $e^{" + str(int(idxmin[0])) + "}, " + "{:.3f}".format(idxmin[1]) + "$"))
    return np.power(np.e, idxmin[0]), idxmin[1]

error_ln_lambda = []
for i in range(-50, 0):
    y_pred, w_vec = get_y_pred_by_analytical(
        x_train=X_TRAIN,
        y_train=Y_TRAIN,
        x=X_TRAIN,  # 训练集
        with_penalty=True,
        lambda_penalty=np.exp(i),
    )
    error_ln_lambda.append(
        [i, calc_e_rms(y_pred=y_pred, y_true=Y_TRAIN), "TRAIN"]
    )  # 训练集上的根均方误差
    y_pred, w_vec = get_y_pred_by_analytical(
        x_train=X_TRAIN,
        y_train=Y_TRAIN,
        x=X_VALIDATION,  # 验证集
        with_penalty=True,
        lambda_penalty=np.exp(i),
    )
    error_ln_lambda.append(
        [i, calc_e_rms(y_pred=y_pred, y_true=Y_VALIDATION), "VALIDATION"]
    )  # 测试集上的根均方误差

best_lambda, least_loss = show_lambda_error(error_ln_lambda)
LAMBDA = best_lambda

y_pred_with_penalty, w_wec_with_penalty = get_y_pred_by_analytical(
    x_train=X_TRAIN, y_train=Y_TRAIN, x=X_TEST, with_penalty=True, lambda_penalty=LAMBDA,
)
W_SOLUTION["Analytical_With_Penalty"] = w_wec_with_penalty
show_comparation(x=X_TEST, y1=Y_TEST, label1="True", y2=y_pred_with_penalty, label2="Pred")
show_train_data(x_train=X_TRAIN, y_train=Y_TRAIN)

对比是否带惩罚项的拟合结果

show_comparation(
    x=X_TEST,
    y1=y_pred_without_penalty,
    label1="Analytical_Without_Penalty",
    y2=y_pred_with_penalty,
    label2="Analytical_With_Penalty",
)
plt.legend(["Analytical_Without_Penalty", "Analytical_With_Penalty"])

<matplotlib.legend.Legend at 0x24269fcaa90>

梯度下降法

def get_y_pred_by_gradient_descent(x_train, y_train, x, lambda_penalty):
    x_vec, t_vec = x_train, y_train
    k, w_vec = gradient_descent_fit(
        x_matrix=get_x_matrix(x_vec, order=ORDER),
        t_vec=t_vec,
        lambda_penalty=lambda_penalty,
        w_vec_0=np.zeros(ORDER + 1),
        learning_rate=LEARNING_RATE,
        delta=DELTA,
    )
    return k, np.dot(get_x_matrix(x, order=ORDER), w_vec), w_vec

k_gradient_descent, y_pred_gradient_descent, w_wec_gradient_descent = get_y_pred_by_gradient_descent(
    x_train=X_TRAIN, y_train=Y_TRAIN, x=X_TEST, lambda_penalty=LAMBDA
)
W_SOLUTION["Gradient_Descent"] = w_wec_gradient_descent
k_gradient_descent

52893

show_comparation(x=X_TEST, y1=Y_TEST, label1="True", y2=y_pred_gradient_descent, label2="Pred")
show_train_data(x_train=X_TRAIN, y_train=Y_TRAIN)

共轭梯度法

def get_y_pred_by_conjugate_gradient(x_train, y_train, x, lambda_penalty):
    x_vec, t_vec = x_train, y_train
    A, x_0, b = switch_deri_func_for_conjugate_gradient(
        x_matrix=get_x_matrix(x_vec, order=ORDER),
        t_vec=t_vec,
        lambda_penalty=lambda_penalty,
        w_vec=np.zeros(ORDER + 1),
    )
    k, w_vec = conjugate_gradient_fit(A=A, x_0=x_0, b=b, delta=DELTA)
    return k, np.dot(get_x_matrix(x, order=ORDER), w_vec), w_vec

k_conjugate_gradient, y_pred_conjugate_gradient, w_wec_conjugate_gradient = get_y_pred_by_conjugate_gradient(
    x_train=X_TRAIN, y_train=Y_TRAIN, x=X_TEST, lambda_penalty=LAMBDA,
)
W_SOLUTION["Conjugate_Gradient"] = w_wec_conjugate_gradient
k_conjugate_gradient

5

show_comparation(x=X_TEST, y1=Y_TEST, label1="True", y2=y_pred_conjugate_gradient, label2="Pred")
show_train_data(x_train=X_TRAIN, y_train=Y_TRAIN)

对比梯度下降和共轭梯度的拟合结果

show_comparation(
    x=X_TEST,
    y1=y_pred_gradient_descent,   label1="Gradient_Descent",
    y2=y_pred_conjugate_gradient, label2="Conjugate_Gradient",
)
show_test_data(x_test=X_TEST, y_test=Y_TEST)
plt.legend(["Gradient_Descent", "Conjugate_Gradient"])

<matplotlib.legend.Legend at 0x24269d1cc88>

四种拟合方法汇总

W_SOLUTION

	Analytical_Without_Penalty	Analytical_With_Penalty	Gradient_Descent	Conjugate_Gradient
0	0.084599	0.103922	0.233123	0.110288
1	13.779036	5.694382	3.978839	5.693982
2	-125.344309	-8.782488	-7.604623	-8.090070
3	601.300442	-9.708832	-3.345961	-12.785307
4	-1540.525403	-1.287225	0.884203	2.779658
5	2029.729045	14.740948	2.743226	14.930732
6	-1308.321291	13.958697	2.525709	10.353532
7	329.531802	-14.476569	0.907662	-12.711493

plt.figure(figsize=(10, 6))
show_train_data(x_train=X_TRAIN, y_train=Y_TRAIN)
show_test_data(x_test=X_TEST, y_test=Y_TEST)
for column in W_SOLUTION.columns:
    sns.lineplot(
        x=X_TEST,
        y=[W_SOLUTION[column] @ get_x_series(x=x, order=ORDER) for x in X_TEST],
    )
plt.legend(["$\sin (2 \pi x)$"] + list(W_SOLUTION.columns))
show_order_and_train_num()