数值计算与优化（共轭梯度法和QR）

# -*- coding: utf-8 -*-
import numpy as np
import time #calculate time of diffient method
np.random.seed(2)   #set seed to make x0 unchanged
def Create_Tridiagonal_Matrices(a,b,c,n):#创建一种特殊的三对角矩阵，主对角线元素b，比主对角线低一行的对角线上a，比主对角线高一行的对角线上cA=np.zeros((n,n))if(b<=c or b<=a or b<a+c or n<2): # 判断是否符合三对角矩阵的基本定义,以及n的个数，1行的三对角毫无意义print "this is not a Create_Tridiagonal_Matrices,please check a,b,c "A[0][0]=b;A[0][1]=cA[n-1][n-1]=b;A[n-1][n-2]=afor j in range(1,n-1):A[j][j-1]=aA[j][j]=bA[j][j+1]=creturn A
def timecmp(A,b,fun1,fun2):#通过运行100次的平均时间较各个函数的效率time1=np.zeros((100,1))time2 = np.zeros((100, 1))for i in range(100):t1=time.clock()fun1(A, b)time1[i][0]=time.clock()-t1t2 = time.clock()fun2(A, b)time2[i][0] = time.clock() - t2return np.mean(time1),np.mean(time2)def conjugate_gradient_method(A, b):m = A.shape[0]x = np.zeros((n, 1))#r0 = b - np.dot(A, x)r = b - np.dot(A, x)p = rwhile (np.sqrt(np.sum(np.power(np.dot(A, x) - b, 2))) / np.sqrt(np.sum(np.power(b, 2))) > 1e-6):#终止条件A_dot_p = np.dot(A, p)  # 后面多处会用到，保存结果减少时间# a=np.dot(r0.T,p)/np.dot(p.T,A_dot_p) 这个地方关于系数a的更新有3种，分别利用不同方式推导得到，都可以用，=。# a=np.dot(r.T,p)/np.dot(p.T,A_dot_p)a = np.dot(r.T, r) / np.dot(p.T, A_dot_p)x = x + a * p #更新xr = r - a * A_dot_p#更新rb0 = - np.dot(r.T, A_dot_p) / np.dot(p.T, A_dot_p)#更新b0p = r + b0 * p #更新p# print np.sqrt(np.sum(np.power(np.dot(A, x) - b, 2)))return xdef Q_R(A, b):n = A.shape[0]#施密特正交基的求解Q = np.zeros((n, n))Q[:, 0] = A[:, 0] / np.sqrt(np.sum(np.power(A[:, 0], 2)))#初始化第一个向量，并且进行了单位化for j in range(1, n): #循环迭代求得所有单位正交基tmp = np.zeros((n, 1))for i in range(j):tmp = tmp + np.dot(A[:, j], Q[:, i].T) * Q[:, i].reshape(n, 1)Q[:, j] = (A[:, j].reshape(n, 1) - tmp).reshape(n, ) / np.sqrt((np.sum(np.power(A[:, j].reshape(n, 1) - tmp, 2))))#求解r的过程，利用rij=<vi,qj>r = np.zeros((n, n))r[0][0] = np.sum(np.power(A[:, 0], 2))for i in range(n):for j in range(i + 1):r[j][i] = np.dot(A[:, i], Q[:, j].T)#计算出来了QRx=b，Q.TQRx=Q.Tb， Rx=Q.tb,R是一个上三角矩阵，可以直接进行求解这个方程组b = np.dot(Q.T, b)newx = np.zeros((n, 1))newx[n - 1] = (b[n - 1] / r[n - 1][n - 1])#迭代法求解xi = n - 2while (i >= 0):sum1 = 0for j in range(i + 1, n):sum1 += r[i][j] * newx[j]newx[i] = (b[i] - sum1) / r[i][i]i -= 1return newxn=100       #可以改变参数，方便调试
A=Create_Tridiagonal_Matrices(-1,2,-1,n)
print"生成的系数矩阵是：", A
b=np.ones((n,1))
print"\n生成的B是", bx=conjugate_gradient_method(A.copy(),b.copy())
print "\nconjugate_gradient_method 的结果是:",x
x=Q_R(A.copy(),b.copy())
print  "\nQ_R_Method 的结果是:",x
print "\n函数计算出来的正确解是",np.linalg.solve(A,b)
t1,t2=timecmp(A,b,conjugate_gradient_method,Q_R)
print "conjugate_gradient_method，Q_R_Method运行100次的平均时间分别为",t1,t2

数值计算与优化（共轭梯度法和QR）相关推荐

机器学习基础（五十九）—— 高级优化算法（梯度下降、L-BFGS、共轭梯度）
优化算法两大核心,一曰:方向,比如由负梯度方向给出:二曰:步长. 在机器学习领域,不管是基础的梯度下降,还是更为高级的 L-BFGS.共轭梯度,都对应的是参数学习,用来学习相关模型的参数. 1. 梯度 ...
【优化理论】共轭梯度下降算法实现
实验目的实验步骤此次实验要求采用共轭梯度下降来解决二次优化问题,在解决此问题前,首先分析下什么是共轭梯度法,共轭梯度法能解决什么问题. 我们来看一个线性方程组Ax=b,求解此方程组的过程可 ...
python实现共轭梯度算法
python实现共轭梯度优化算法一.共轭梯度算法简介二.实现共轭梯度方法的两块重要积木 1.共轭方向的确定 2.方向优化步长的确定 note 三.共轭梯度算法优化过程四.python实现共轭梯度 ...
机器学习：共轭梯度算法（PCG）
今天介绍数值计算和优化方法中非常有效的一种数值解法,共轭梯度法.我们知道,在解大型线性方程组的时候,很少会有一步到位的精确解析解,一般都需要通过迭代来进行逼近,而 PCG 就是这样一种迭代逼近算法. ...
几种常用的优化方法梯度下降法、牛顿法、）
几种常用的优化方法 1. 前言熟悉机器学习的童鞋都知道,优化方法 ...
CG共轭梯度下降法【学习笔记、例题与代码】
资料参考视频: [详细推导][本视频还证明了收敛性]https://www.bilibili.com/video/BV16a4y1t76z?from=search&seid=35153938 ...
matlab ncg,2步非线性共轭梯度(NCG)法
2步非线性共轭梯度(NCG)法提出一种求解非线性方程组的2步非线性共轭梯度法(2步NCG法),并在一定条件下证明 (本文共9页) 阅读全文>> 作者提出的自适应共轭梯度方法是对共轭梯度方 ...
共轭梯度算法求最小值-scipy
1 # coding=utf-8 2 3 #共轭梯度算法求最小值 4 import numpy as np 5 6 from scipy import optimize 7 8 9 10 11 def ...
入门 | 一文简述深度学习优化方法——梯度下降
http://www.sohu.com/a/241298990_129720 本文是一篇关于深度学习优化方法--梯度下降的介绍性文章.作者通过长长的博文,简单介绍了梯度下降的概念.优势以及两大挑战.文 ...

数值计算与优化（共轭梯度法和QR）

数值计算与优化（共轭梯度法和QR）相关推荐

最新文章

热门文章