本题中我使用到了最速下降法来拟合函数，一开始我使用了高斯牛顿法计算了第一道题目，但是在代入最小二乘那个矩阵的公式时发现矩阵不可逆，因此直接否定了这个方法不适用于这道题目，后来决定使用最速下降法，因为该方法适用场合很广且十分简便，但是收敛速度不佳。算法思路十分简单，直接把上面的问题看成是无约束优化问题，寻找残差方程求和的最小值，首先求出该方程的梯度，然后随机选个初始点，一直沿着梯度方向迭代即可。
通过程序计算，最后优化执行的结果的精确度较高，运行速度很慢，我使用的一维线搜索方法采用了精确线搜索黄金分割法，初始边界取了0到3。设置的精度是1e-5，大约执行了8万次迭代耗时5分钟才得到最后结果，求出的结果，以及点带入残差方程值如下：
收敛性简单分析：首先通过表面运行时间来看发现这个算法的收敛速度是非常慢的，并且随的精度的提高，运行的时间越来越长，比如说我设置不同的epsion值，1e-4时迭代了20000次左右，1e-5迭代了80000次左右。
在一般情况下，当用最速下降法寻找极小点时，其搜索路径呈直角锯齿状，如下图，
这也就意味着开头几步下降的快，但是越到后面，下降速度越来越慢，通过证明可以的知梯度下降法的收敛速度是线性的。
在一维线搜索中，我采用了精确性搜索，这存在一些问题，首先是速度慢，而且在这里只是追求了该点下降最陡的方向，但是并不意味着该点就是下降最多的点，而且算法会因此变得不灵活因为，程序每次都要求出fai(x)=f(xk+akdk)的最小值，如果遇到更加复杂的函数，会有很大的麻烦，有可能找不到这个最小值。
import math
import numpy as np
e=math.e
def emi(a,b):
return a*(e**b)
def t(k):
return 0.1*k
def y(k):
return emi(1,-t(k))-emi(5,-10*t(k))+emi(3,-4*t(k))
def function(x): #返回函数即残差方程的求和，传入的x是array类型或matrix 6维向量
if (type(x) == np.matrix):
q = np.array(x)
x1,x2,x3,x4,x5,x6 = q[0][0],q[0][1],q[0][2],q[0][3],q[0][4],q[0][5]
else:
x1,x2,x3,x4,x5,x6 = x[0],x[1],x[2],x[3],x[4],x[5]
sum=0
for k in range(1,8):
sum+= (emi(x3,-t(k)*x1)-emi(x4,-t(k)*x2)+emi(x6,-t(k)*x5)-y(k))**2
return sum
def r(x,k): #返回点代入残差方程的值
if (type(x) == np.matrix):
q = np.array(x)
x1, x2, x3, x4, x5, x6 = q[0][0], q[0][1], q[0][2], q[0][3], q[0][4], q[0][5]
else:
x1, x2, x3, x4, x5, x6 = x[0], x[1], x[2], x[3], x[4], x[5]
return (emi(x3,-t(k)*x1)-emi(x4,-t(k)*x2)+emi(x6,-t(k)*x5)-y(k))**2
def grade(x): #返回梯度值
if (type(x) == np.matrix):
q = np.array(x)
x1, x2, x3, x4, x5, x6 = q[0][0], q[0][1], q[0][2], q[0][3], q[0][4], q[0][5]
else:
x1, x2, x3, x4, x5, x6 = x[0], x[1], x[2], x[3], x[4], x[5]
a1,a2,a3,a4,a5,a6=0,0,0,0,0,0
for k in range(1,8):
a1 += 2 * (emi(x3, -t(k) * x1) - emi(x4, -t(k) * x2) + emi(x6, -t(k) * x5) - y(k)) * (emi(-t(k) * x3, -t(k) * x1))
a2 += 2*(emi(x3, -t(k) * x1) - emi(x4, -t(k) * x2) + emi(x6, -t(k) * x5) - y(k)) * (emi(t(k) * x4, -t(k)*x2))
a3 += 2*(emi(x3, -t(k) * x1) - emi(x4, -t(k) * x2) + emi(x6, -t(k) * x5) - y(k)) * (emi(1, -t(k) * x1))
a4 += 2*(emi(x3, -t(k) * x1) - emi(x4, -t(k) * x2) + emi(x6, -t(k) * x5) - y(k)) * (emi(-1, -t(k) * x2))
a5 += 2 * (emi(x3, -t(k) * x1) - emi(x4, -t(k) * x2) + emi(x6, -t(k) * x5) - y(k)) * (emi(-t(k) * x6, -t(k) * x5))
a6 += 2 * (emi(x3, -t(k) * x1) - emi(x4, -t(k) * x2) + emi(x6, -t(k) * x5) - y(k)) * (emi(1, -t(k) * x5))
return np.matrix([a1,a2,a3,a4,a5,a6])
def fai(p, x, d): #用于一维搜索的函数
q = x + p * d
return function(q)
def gold(x, d): #黄金分割法
a, b, k, e = [0], [3], 0, 1e-5
l, r = [a[0] + 0.382 * (b[0] - a[0])], [a[0] + 0.618 * (b[0] - a[0])]
fl, fr = [fai(l[0], x, d)], [fai(r[0], x, d)]
while (1):
if (fl[k] > fr[k]):
if (b[k] - l[k] <= e):
return r[k]
a.append(l[k])
b.append(b[k])
l.append(r[k])
fl.append(fr[k])
r.append(a[k + 1] + 0.618 * (b[k + 1] - a[k + 1]))
fr.append(fai(r[k + 1], x, d))
else:
if (r[k] - a[k] <= e):
return l[k]
a.append(a[k])
b.append(r[k])
r.append(l[k])
fr.append(fl[k])
l.append(a[k + 1] + 0.382 * (b[k + 1] - a[k + 1]))
fl.append(fai(l[k + 1], x, d))
k += 1
eps = 1e-5
x = [np.array([1, 1, 1, 1,1,1])]
g = []
d = []
alpha=[]
k = 0
while (1):
print("第%d次迭代："%k)
print(x[k])
g.append(grade(x[k]))
if(np.linalg.norm(g[k])<=eps):
break
d.append(-g[k])
alpha.append(gold(x[k],d[k]))
x.append(x[k]+alpha[k]*d[k])
k += 1
for i in range (1,8):
print("第%d个点代入拟合函数的偏差值为：%lf"%(i,r(x[k],i)))
