基于梯度下降法的傅里叶级数拟合

#运算过程中主要运用矩阵（并行）运算，以提高运算速度

#梯度下降算法
#傅里叶拟合
import numpy as np
import math
import matplotlib.pyplot as plt

x=np.linspace(-math.pi,math.pi,200)
y=np.sin(x)#原函数
n=10
x1=np.linspace(1,n,n)
a=np.random.rand(2*n+1)#生成随机系数矩阵
grad_a=np.zeros(2*n+1)#用于保存预测值
learning_rate=1e-4

f1=np.array(x**0)#第一个系数
x1.shape=(n,1)#行向量变列向量
x.shape=(1,200)
f2=np.cos(x*x1)#广播乘法（并不是矩阵乘法）
f3=np.sin(x*x1)
f4=np.vstack([f1,f2,f3])#矩阵合并
f4=f4.transpose()#矩阵转置
#开始迭代
for i in range(3000):
    f4=f4.transpose()
    y_pred=np.dot(a,f4)#矩阵乘法，计算预测值
    #计算损失函数
    loss=np.square(y_pred-y).sum()#square获取矩阵元素的平方
    #每隔100项输出一次
    if i%100==99:
        print(i,loss)
        
    grad_y_pred=2.0*(y_pred-y) 
    f4=f4.transpose()
    grad_a=np.dot(grad_y_pred,f4)#矩阵乘法
    a-=grad_a*learning_rate
    
print(f'常数项：{a[0]}')
for i in range(n):
    print(f'cos({i+1}x)的系数为:{a[i+1]}')
for i in range(n):
    print(f'sin({i+1}x)的系数为:{a[i+n+1]}')
    
#使绘图可以显示中文
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
#原函数
x.shape=(200,)
plt.plot(x,y)
plt.title('原函数')
plt.show()
plt.plot(x,y_pred,'y')
plt.title('傅里叶级数拟合图像')
plt.show()