哈工大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对象,可以理解为一个多维数组(本实验中仅会用到一维(可以当作列向量)和二维(矩阵))。下面用小写的

x

\pmb x

xx表示列向量,大写的

A

A

A表示矩阵。A.T表示

A

A

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逐元素操作的:

APK名称修改时间版本

>>> 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.lognp.exp等与python的math库相似的函数(只不过是对多维数组进行逐元素运算)。

标准模板库

np.dot

返回两个矩阵的乘积。与线性代数中的矩阵乘法一致。要求第一个矩阵的列等于第二个矩阵的行数。特殊地,当其中一个为一维数组时,形状会自动适配为

n

×

1

n\times1

n×1

1

×

n

.

1\times n.

1×n.

Multi-Head Atte

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

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=1N[yif(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=

111x1x2xNx12x22xN2x1mx2mxNm

N×(m+1),Y=

y1y2yN

N×1,W=

w0w1wm

(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)y2f(xN)yN

=XWY.
因此,损失函数

L

=

(

X

W

Y

)

T

(

X

W

Y

)

.

L=(XW-Y)^T(XW-Y).

L=(XWY)T(XWY).
(为了求得向量

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}

WL=W[(XWY)T(XWY)]=W[(WTXTYT)(XWY)]=W(WTXTXWWTXTYYTXW+YTY)=W(WTXTXW2YTXW+YTY)(容易验证,WTXTY=YTXW,因而可以将其合并)=2XTXW2XTY
说明:
(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.

stylegan

矩阵求导线性代数课上也没有系统教过,只对这里出现的做一下说明。(多了我也不会
令偏导数为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的解析解,我们只需要调用函数算出这个值即可。

limit

'''
最小二乘求出解析解, 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库还是挺不友好的)

3D

简单地验证一下我们已经完成的函数的结果:为此,我们先写一个draw函数,用于把求得的

W

W

W对应的多项式

f

(

x

)

f(x)

f(x)画到pyplot库的图像上去:

gateway

'''
绘制给定系数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次多项式拟合的效果还是比较不错的(数据集每次随机生成,所以跟第一幅图不一样)。

低代码

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

BCELoss

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()

补充说明

上面有一块不太严谨:对于一个矩阵

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.

Istio

添加正则项(岭回归)

最小二乘法容易造成过拟合。为了说明这种缺陷,我们用所生成数据集的前50个点进行训练(这样抽样不够均匀,这里只是为了说明过拟合),得出参数,再画出整个函数图像,查看拟合效果:

MODBUS

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=(XWY)T(XWY)+λ∣∣W22
其中

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内的参数过大。

iOS 启动图适配

举个例子(数是随便编的):当正则化系数为

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范数。

SPL

重复上面的推导,我们可以得出解析解为

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)也是可逆的。

电赛

该部分代码如下。

Mask R-CNN

'''
岭回归求解析解, 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 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()

在这里插入图片描述
上图标明了

x

x

x随着迭代的演进,可以看到

x

x

x不断沿着正半轴向零点靠近。需要注意的是,学习率不能过大(虽然在上面的程序中,学习率设置得有点小了),需要手动进行尝试调整,否则容易想象,

x

x

x在正负半轴来回震荡,难以收敛。

在最小二乘法中,我们需要优化的函数是损失函数

L

=

(

X

W

Y

)

T

(

X

W

Y

)

.

L=(XW-Y)^T(XW-Y).

L=(XWY)T(XWY).
下面我们用梯度下降法求解该问题。在上面的推导中,

L

W

=

2

X

T

X

W

2

X

T

Y

,

\begin{aligned}\frac{\partial L}{\partial W}=2X^TXW-2X^TY\end{aligned},

WL=2XTXW2XTY,
于是我们每次在迭代中对

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)=21xxTAxxbbTxx+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)=bAx(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}.

105.
下面我们按照这个过程实现代码:

'''
共轭梯度法(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):
blog.csdnimg.cn/49f5b3380f1d45e48033c94208ed2b2c.png)
最后附上四种方法的拟合图像(基本都一样)和主函数,可以根据实验要求调整参数:
在这里插入图片描述

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()

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注