# 要求使用共轭梯度法和QR分解法求解方程组
# 分析
## 共轭梯度法输入：Ax=b想法：构造迭代关系$ x^k = x^{k−1} + a^{k−1}p^{k−1} $其中$ \{p^0,....,p^k\} $都共轭于A，实际意义上代表每步迭代的方向，而a表示对应p上迭代的步长，原理是将每个方向先走到最终，所以我们总共只需要n步完成全部迭代。经过一系列运算得$ r_k=r_{k-1}-a_{k-1}Ap{k-1} $ $a_{k-1}=/$ $p^k=r^k-/*p^{k-1}$所以算法如下
## QR分解法
输入：Ax=b
想法：将A分解为QR，其中Q为正交矩阵，R为上三角矩阵，所以$$ Rx=Q^Tb $$在自底部向上求解x。其中分解为QR有2种算法HouseHolder法和Gram_Schmidt法。
# 实现
## 共轭梯度法
```python
#ConjugateGradient Solve
def ConjugateGradient(A,x,b):
r0=b-np.dot(A,x)
p0=r0
for k in range(0,A.shape[0]):
a0=np.dot(r0.T,r0)/np.dot(p0.T,np.dot(A,p0))
x=x+a0*p0
r0=r0-a0*np.dot(A,p0)
p0=r0-np.dot(r0.T,np.dot(A,p0))/np.dot(p0,np.dot(A,p0))*p0
rate=np.sqrt(np.sum((np.dot(A,x)-b)**2))/np.sqrt(np.sum(b**2))
if(rate<10**-6):
print k
return x
```
## QR法
```python
#Based on Gram-Schmidt method decomposite
def QR_Gram_Schmidt(A):
Q = np.zeros_like(A)
k = 0
for a in A.T:
u = np.copy(a)
for i in range(0,k):
u -= np.dot(np.dot(Q[:,i].T,a),Q[:,i])
e = u / np.linalg.norm(u)
Q[:,k] = e
k += 1
R =np.dot(Q.T,A)
return Q,R
#based on Householder method decomposite
def QR_Householder(A):
(r, c) = np.shape(A)
Q = np.identity(r)
R = np.copy(A)
for k in range(r - 1):
x = R[k:, k]
e = np.zeros_like(x)
e[0] = np.linalg.norm(x)
u = x - e
v = u / np.linalg.norm(u)
Q_cnt = np.identity(r)
Q_cnt[k:, k:] -= 2.0 * np.outer(v, v)
R = np.dot(Q_cnt, R) # R=H(n-1)*...*H(2)*H(1)*A
Q = np.dot(Q, Q_cnt) # Q=H(n-1)*...*H(2)*H(1) H为自逆矩阵
return Q,R
#calculate Rx=Q.T b
def CalcQR(A,b,i):
timea=time.time()
if i==0:
time0=time.time()
Q,R=QR_Gram_Schmidt(A)
print "Gram_Schmidt",time.time()-time0
else:
time0=time.time()
Q,R=QR_Householder(A)
print "Householder",time.time()-time0
x=np.dot(Q.T,b)
m=x.shape[0]
for k in range(m-1,-1,-1):
x[k]=(x[k]-np.dot(R[k,(k+1):m],x[(k+1):m]))/R[k,k]
print "QR",i,"time :",time.time()-timea
return x
```
## Main
```python
if __name__=='__main__':
A=generateA()
b=generateb()
# print np.dot(np.linalg.inv(A),b)
x=np.zeros(100).T
time0=time.time()
x_result=ConjugateGradient(A,x,b)
print "ConjugateGradient",time.time()-time0
x_result
print CalcQR(A,b,0)
print CalcQR(A,b,1)
```
# 结果分析
k=49
ConjugateGradient 0.00170302391052
Gram_Schmidt 0.0211570262909
QR 0 time : 0.0214040279388
Householder 0.0368831157684
QR 1 time : 0.0372591018677
所以共轭梯度法的速度远远高于QR法，且他们都获得的是真实解，而在QR分解上Gram_Schmidt优于Householder
网友评论