牛顿迭代法
原标题:牛顿迭代法
原文来自:CSDN 原文链接:https://blog.csdn.net/qq_42002616/article/details/103794196
牛顿法是牛顿在17世纪提出的一种求解方程f(x)=0.多数方程不存在求根公式,从而求精确根非常困难,甚至不可能,从而寻找方程的近似根就显得特别重要。
牛顿迭代法是取x0之后,在这个基础上,找到比x0更接近的方程的跟,一步一步迭代,从而找到更接近方程根的近似跟。方法使用函数f(x)的泰勒级数的前面几项来寻找方程f(x) = 0的根。牛顿迭代法是求方程根的重要方法之一,其最大优点是在方程f(x) = 0的单根附近具有平方收敛,而且该法还可以用来求方程的重根、复根。另外该方法广泛用于计算机编程中。
设r是f(x)=0的根,选取x0作为r初始近似值,过点(x0,f(x0))做曲线y=f(x)的切线L,L的方程为y=f(x0)+f’(x0)(x-x0),求出L与x轴交点的横坐标 x1=x0-f(x0)/f’(x0),称x1为r的一次近似值,过点(x1,f(x1))做曲线y=f(x)的切线,并求该切线与x轴的横坐标 x2=x1-f(x1)/f’(x1)称x2为r的二次近似值,重复以上过程,得r的近似值序列{Xn},其中Xn+1=Xn-f(Xn)/f’(Xn),称为r的n+1次近似值。上式称为牛顿迭代公式。
''' f(x)=0的x值 x0处的切线方程L:y=f(x0)+f'(x0)(x-x0) 令y=f(x0)+f'(x0)(x-x0)=0求出L与x轴的交点x1 x1=x0-f(x0)/f'(x0) 之后一直迭代求到xn,满足精度即可''' '''求方程x^3+2x^2+3x+4=0在1附近的实数根''' '''写成函数形式''' def f(x): f=x**3+2*x**2+3*x+4 # f的方程 return f def f_first_order(x): f1=3*x**2+4*x+3 # f的一阶导数 return f1 def get_root(x0, max_iter=60, tol=1e-7): # 将初始值浮点化 p0 = x0 * 1.0 for i in range(max_iter): # f的一阶导数不能为0 p = p0 - f(p0) / f_first_order(p0) # 如果小于精度值则退出迭代 if abs(p - p0) < tol: # tol是判断迭代更新的阈值 return u'经过%s次的迭代,我们估计的参数值是%s' % (i + 1, p) p0 = p print(u'已达到最大迭代次数, 但是仍然无法收敛') if __name__ == '__main__': print(get_root(1)) # 由于牛顿迭代方法是局部最优解,不同的初始值有不同的结果。初始值分别取2、0、-2
'''由迭代次数'''x0=x1=1for i in range(60): f = x0 ** 3 + 2 * x0 ** 2 + 3 * x0 + 4 f1 = 3 * x0 ** 2 + 4 * x0 + 3 x1=x0-f/f1 x0=x1print(x1)
'''由迭代精度''' x0 = 1 #在1附近的实数根 x1=0.2 #第一次迭代的值(这个必须要有,且不能太接近于x0) while abs(x1-x0)>=1e-9: x0 = x1 f = x0 ** 3 + 2 * x0 ** 2 + 3 * x0 + 4 f1 = 3 * x0 ** 2 + 4 * x0 + 3 x1=x0-f/f1 print(x1) #-1.6506291914393882
'''由迭代次数和由迭代精度''' x0=x1=1 for i in range(60): f = x0 ** 3 + 2 * x0 ** 2 + 3 * x0 + 4 f1 = 3 * x0 ** 2 + 4 * x0 + 3 x1=x0-f/f1 if abs(x1-x0)<=1e-9: print(u'经过%s次的迭代,我们估计的参数值是%s' %(i + 1,x1)) break elif i>=60: print(u'已达到最大迭代次数, 但是仍然无法收敛') break x0=x1 print(x1)
import numpy as np from math import pow #导入训练数据 def load_data(file_path): '''导入数据 input:file_path(string):训练数据 output:feature(mat):特征 label(mat)标签 ''' f=open(file_path) feature=[] label=[] for line in f.readlines(): feature_tmp=[] lines=line.strip().split('t') #以tab键隔开 feature_tmp.append(1) #x0 for i in range(len(lines)-1): feature_tmp.append(float(lines[i])) feature.append(feature_tmp) label.append(float(lines[-1])) f.close() return np.mat(feature),np.mat(label).T def least_square(feature,label): '''最小二乘法 input: feature(mat):特征 label(mat):标签 output: w(mat):回归系数 ''' w = (feature.T * feature).I * feature.T * label # W=(X^TX)^(-1)*(X^TY) return w #一阶导数 def first_derivativ(feature,label,w): '''计算一阶导函数的值 input: feature(mat):特征 label(mat):标签 output: g(mat):一阶导数值 ''' m,n=np.shape(feature) g=np.mat(np.zeros((n,1))) for i in range(m): err=label[i,0]-feature[i,]*w for j in range(n): g[j,]-=err*feature[i,j] return g #二阶导数 def second_derivarive(feature): '''计算二阶导函数的值 input: feature(mat):特征 output: G(mat):二阶导数值 ''' m, n = np.shape(feature) G=np.mat(np.zeros((n,n))) for i in range(m): x_left = feature[i, ].T x_right = feature[i, ] G += x_left * x_right return G #损失函数的计算 def get_error(feature,label,w): '''计算误差 input: feature(mat):特征 label(mat):标签 w(mat):线性回归模型的参数 output: 损失函数值 ''' return (label-feature*w).T*(label-feature*w)/2 #全局牛顿法 def newton(feature,label,iterMax,sigma,delta): '''牛顿法 input: feature(mat):特征 label(mat):标签 iterMax(int):最大迭代次数 sigma(float), delta(float):牛顿法中的参数 output: w(mat):回归系数 ''' n=np.shape(feature)[1] w=np.mat(np.zeros((n,1))) #n行,1列 it=0 while it<=iterMax: print(it) g=first_derivativ(feature,label,w) #调用一阶导数 G=second_derivarive(feature) #调用二阶导数 d=-G.I*g m = get_min_m(feature, label, sigma, delta, d, w, g) # 得到最小的m w=w+pow(sigma,m)*d if it % 10==0: print("t---- itration: ", it, " , error: ", get_error(feature, label, w)[0, 0]) it+=1 return w def get_min_m(feature, label, sigma, delta, d, w, g): '''计算步长中最小的值m input: feature(mat):特征 label(mat):标签 sigma(float),delta(float):全局牛顿法的参数 d(mat):负的一阶导数除以二阶导数值 g(mat):一阶导数值 output: m(int):最小m值 ''' m = 0 while True: w_new = w + pow(sigma, m) * d left = get_error(feature, label, w_new) right = get_error(feature, label, w) + delta * pow(sigma, m) * g.T * d if left <= right: break else: m += 1 return m #保存模型的 save_model 函数 def save_model(file_name, w): '''保存最终的模型 input: file_name(string):要保存的文件的名称 w(mat):训练好的线性回归模型 ''' f_result = open(file_name, "w") m, n = np.shape(w) for i in range(m): w_tmp = [] for j in range(n): w_tmp.append(str(w[i, j])) f_result.write("t".join(w_tmp) + "n") f_result.close() import matplotlib.pyplot as plt # show your trained logistic regression model only available with 2-D data def showLogRegres(weights, feature, label): # notice: train_x and train_y is mat datatype import matplotlib.pyplot as plt numSamples, numFeatures = np.shape(feature) #此处shape应加上np # if numFeatures != 3: # print("Sorry! I can not draw because the dimension of your data is not 2!") # return 1 # draw all samples for i in range(numSamples): plt.plot(feature[i, 1], label[i], 'or') # y1=weights*feature[i] # y_yuce.append(y1) # plt.plot(feature[:1],y_yuce) # plt.savefig('nihe1,jpg') # plt.show() # draw the classify line min_x = min(feature[:, 1])[0, 0] max_x = max(feature[:, 1])[0, 0] weights = weights.getA() # convert mat to array y_min_x = float(- weights[0]+weights[1] * min_x) y_max_x = float(- weights[0]+weights[1] * max_x) plt.plot([min_x, max_x], [y_min_x, y_max_x], '-g') plt.xlabel('X1') plt.ylabel('X2') plt.savefig('nihe_newton.jpg')#可以保存下图像 plt.show() #------------ print(weights) # [[0.00310499] # [0.99450247]] if __name__ == "__main__": # 1、导入数据集 print("----------- 1.load data ----------") feature, label = load_data("data.txt") # 2.1、最小二乘求解 print("----------- 2.training ----------") # print "t ---------- least_square ----------" # w_ls = least_square(feature, label) # 2.2、牛顿法 print("t ---------- newton ----------") w_newton = newton(feature, label, 50, 0.1, 0.5) # 3、保存最终的结果 print("----------- 3.save result ----------") save_model("weights ", w_newton) showLogRegres(w_newton,feature, label)
运行结果如图所示,红色点为数据散点,绿色的线为拟合曲线
使用数据如下所示,命名为data.txt
免责声明:本文来自互联网新闻客户端自媒体,不代表本网的观点和立场。
合作及投稿邮箱:E-mail:editor@tusaishared.com
热门资源
Python 爬虫(二)...
所谓爬虫就是模拟客户端发送网络请求,获取网络响...
TensorFlow从1到2...
原文第四篇中,我们介绍了官方的入门案例MNIST,功...
TensorFlow从1到2...
“回归”这个词,既是Regression算法的名称,也代表...
TensorFlow2.0(10...
前面的博客中我们说过,在加载数据和预处理数据时...
机器学习中的熵、...
熵 (entropy) 这一词最初来源于热力学。1948年,克...
智能在线
400-630-6780
聆听.建议反馈
E-mail: support@tusaishared.com