问题: 梯度下降法之使用各种函数(如多项式函数傅里叶级数)拟合正弦曲线

?1.三次函数拟合正弦余弦函数
?2.多项式函数拟合正弦余弦函数
?3.傅里叶级数拟合正弦余弦函数
?4.管他什么函数拟合什么函数
??理解了通通拿下好嘛 ! 怎么修改拟合函数见文末 q(≧▽≦q)
??ps：要是感觉博文乱没有思绪，拿鼠标对照左侧目录食用
?
??此处以三次函数为例，其他的函数拟合同理

?
?
?
?

1.梯度下降法原理

??梯度下降相关公式
?
$拟合的函数:h(x)=\sum _{i=0}^n{\theta }_i{x}_i$
$损失函数:J(\theta )=\frac{1}{2m}\sum _{i=1}^m(y^{(i)}?h_{\theta }(x^{(i)}){)}^2$
$求偏导:\frac{?J(\theta )}{?{\theta }_j}=\frac{1}{m}\sum _{i=1}^m(y^{(i)}?h_{\theta }(x^{(i)}))x_j^{(i)}$

$迭代公式:{\theta }_j={\theta }_j?\alpha \frac{?J(\theta )}{?{\theta }_j}$

??原理不再赘述, 请参考文档: https://blog.csdn.net/lilyth_lilyth/article/details/8973972
?
?
?
?
?

2.分析问题, 解题思路

??把公式写成矩阵形式(嫌麻烦，这里只打出了部分) :

$$\Theta ={?\begin{array}{c}{\theta }_1\\ {\theta }_2\\ {\theta }_3\\ {\theta }_4\end{array}?}_{4\times 1}X={?\begin{array}{cccc}{x}_0^{(1)}& x_1^{(1)}& x_2^{(1)}& x_3^{(1)}\\ x_0^{(2)}& x_1^{(2)}& x_2^{(2)}& x_3^{(2)}\\ ?& ?& ?& ?\\ x_0^{(m)}& x_1^{(m)}& x_2^{(m)}& x_3^{(m)}\end{array}?}_{m\times 4}Y={?\begin{array}{c}{y}_1\\ y_2\\ ?\\ y_m\end{array}?}_{m\times 1}$$

$$H(X)=X_{m\times 4}{\Theta }_{4\times 1}={?\begin{array}{c}h(x^{(1)})\\ h(x^{(2)})\\ ?\\ h(x^{(m)})\end{array}?}_{m\times 1}D=H(X)?Y={?\begin{array}{c}h(x^{(1)})?y_1\\ h(x^{(2)})?y_2\\ ?\\ h(x^{(m)})?y_m\end{array}?}_{m\times 1}$$

???注：这个$x_2^{(i)}$表示第 i 个样本值的 第2个分量， 大写通常指矩阵
?
?

??代码需要用到的矩阵 :
?
?????$拟合值:H(X)=X\Theta 残差:D=H(X)?Y$

?????$损失函数:J(\Theta )=\frac{1}{2m}{D}^{T}D梯度:\frac{?J(\Theta )}{?\Theta }=\frac{1}{m}{X}^{T}D$

?????$迭代公式:\Theta =\Theta ?\alpha \frac{?J(\Theta )}{?\Theta }$

???注：待求向量 $\Theta$，样本矩阵 $X$，实际值向量 $Y$ ，拟合值向量 $H(X)$，残差 $D$，学习率 $\alpha$，样本点个数 $m$
把这些还有上面的名字记好咯，别后面认不得了。
?
?
?
?

详细的解题步骤及思路 :
?
??先写一段废话, 都用Python了嘛, 处理数据就少用for循环, 多用用numpy, 要不然CPU跑起来可以烤萝卜吃了。
所以说就有了上面把梯度相关公式写成矩阵形式。
然后嘞， 我们用 ${\Theta }_{4\times 1}$ 来对应储存 a, b, c, d 。
$x\_vector=[x_0,x_1,x_2,x_3]$ 来对应储存 $x^0，x，x^2，x^3$，
然后我们一乘 $(x\_vector)\Theta$, 这不就是 $h(x)=a+bx+c{x}^2+d{x}^3$ 嘛。 (感谢我神凯利)
废话就不多说了，下面正经走起。。。。
?
?
?
?
?

1.取样本点

??先从区间(-π, π) 上取一些离散的样本点 $x^{(i)}$，用以带入 $h(x)=a+bx+c{x}^2+d{x}^3$ 来逼近 $y=sinx$

domain = np.arange(-np.pi, np.pi, 0.1)  # 设置定义域

?

2.先创建用于输出数组的函数

??输入一个离散点$x^{(i)}$，对应运算后, 得到关于这个$x^{(i)}$ 的四个分量，
并储存在数组 x_vector = $[x_0^{(i)},x_1^{(i)},x_2^{(i)},x_3^{(i)}]$ 中， 该题中即为向量 $[(x^{(i)}{)}^0,x^{(i)},(x^{(i)}{)}^2,(x^{(i)}{)}^3]$
??注 : 这一步很重要, 也是程序能扩展的核心代码
比如后面的用傅里叶级数拟合 $sinx$

x_vector = lambda x: np.array([ x**0, x, x**2, x**3 ])  # 为简洁, 使用lambda表达式

?

3.创建列向量$\Theta$

??用于对应储存待定常数a, b, c, d
直接用np.zeros()初始化这个数组 , 然后再用.reshape()给它转化成4x1的数组
这里的参数4是因为样本点有4个分量, 即 $[x_0^{(i)},x_1^{(i)},x_2^{(i)},x_3^{(i)}]$
??注 : 此步重要程度及原因同步骤2

theta = np.zeros(4).reshape(4, 1)

?

4.创建样本矩阵$X$

??把你第1步取好的离散点数组 domain 带入函数 x_vector ，你会得到一个 4×m 的数组，
你用.transpose()给它转一下，样本矩阵 x 就是m×4的了。

m = len(domain)
alpha = 0.01                            # 学习率
error = 1                               # 初始化拟合误差
x = hx(domain).transpose()              # x的样本矩阵 
y = np.sin(domain).reshape(m, 1)        # 被拟合的函数

?

5. 定义损失函数和梯度函数

??不用多说了吧，参照上面给的损失函数矩阵，对应乘就行了。
但是我还是要说一下，fun_error用来后面主程序中确定拟合精度用的，fun_error值越小，h(x)拟合sinx的效果越好

def fun_error(theta, x, y):             # 定义损失函数
    diff = np.dot(x, theta) - y			# 残差 diff
    return (1/2m)* np.dot(diff.transpose(), diff)

def fun_gradient(theta, x, y):          # 定义梯度函数
    diff = np.dot(x, theta) - y			# 残差 diff
    return (1/m) * np.dot(x.transpose(), diff)

?

6. 主程序：迭代$\Theta$了

??根据上面给的迭代公式算就行了。

while abs(error) >= 0.01:               # 主程序: 梯度下降
    error = fun_error(theta, x, y)
    gradient = fun_gradient(theta, x, y)
    theta = theta - alpha * gradient

?

7. 最后呢，就画个图乐一下吧

??毕竟对于老司机来说，图像更有视觉冲击力

print('所求参数向量θ为:\n', theta)
h = np.dot(x, theta)
plt.plot(domain, y, "*")
plt.plot(domain, h)
plt.show()

?
?
?
?
?

3. 完整代码

import numpy as np
import matplotlib.pyplot as plt

###############此段代码根据拟合题目修改######################
x_vector = lambda x: np.array([ x**0, x, x**2, x**3 ])
theta = np.zeros(4).reshape(4, 1)       # 待求参数Θ 此处4为变量个数
domain = np.arange(-np.pi, np.pi, 0.1)  # 设置定义域
#########################################################

m = len(domain)
alpha = 0.01                            # 学习率
error = 1                               # 初始化拟合误差
x = x_vector(domain).transpose()              # x的样本矩阵 
y = np.sin(domain).reshape(m, 1)        # 被拟合的函数

def fun_error(theta, x, y):             # 定义损失函数
    diff = np.dot(x, theta) - y			# 残差 diff
    return (1/2m)* np.dot(diff.transpose(), diff)

def fun_gradient(theta, x, y):          # 定义梯度函数
    diff = np.dot(x, theta) - y			# 残差 diff
    return (1/m) * np.dot(x.transpose(), diff) 
                             
while abs(error) >= 0.01:               # 主程序: 梯度下降
    error = fun_error(theta, x, y)
    gradient = fun_gradient(theta, x, y)
    theta = theta - alpha * gradient

print('所求参数向量θ为:\n', theta)
h = np.dot(x, theta)
plt.plot(domain, y, "*")
plt.plot(domain, h)
plt.show()

?
?
?
?
?

4.关于题目的扩展及总结

??什么用 $h(x)=a+bx+c{x}^2+d{x}^3+sinx$ 啊，傅里叶拟合正弦函数或者拟合其他函数啊，改一改 $x\_vector$，$\Theta$ 等参数(见下列代码块)就行了。

x_vector = lambda x: np.array([ x**0, x, x**2, x**3, np.sin(x) ])
theta = np.zeros(5).reshape(5, 1)       # 待求参数Θ 

x_vector = lambda x: np.array([ x**0, np.sin(x), np.cos(x), np.sin(2*x), np.cos(2*x) ])
theta = np.zeros(5).reshape(5, 1)       # 待求参数Θ 

??用数组或者矩阵来处理梯度下降这个问题简洁且快速。
?
?
?
?
?

写在最后的话

??第一次发CSDN，敲打了一下午公式，排版了一晚上，话有点多有点啰嗦，看完一遍感觉有思路点乱可以结合目录再看一遍，要还是感觉乱还请多多包涵。观众老爷给个三连吧。
?
?
本篇也参考了其它一些文章，特此感谢！
参考文档：
https://blog.csdn.net/yhao2014/article/details/51554910
https://blog.csdn.net/lilyth_lilyth/article/details/8973972
https://blog.csdn.net/asahinokawa/article/details/80846439
?
?
?