【Julia1.5.2】求解对称正定线性方程组的共轭梯度法

流星雨

续watermelon主题：https://www.0xaa55.com/thread-26183-1-1.html中的雅各比迭代法，整理另一种“半迭代半直接”的求解对称正定线性方程组的共轭梯度法。

参考数目：
【1】《数值分析》(李乃成、梅立泉编著，科学出版社，ISBN：978-7-03-032192-3)
采用Julia语言(1.5.2)编制脚本Conjugate Gradient，主要功能均为求解方程组Ax=b。
脚本中的函数Conjugate_Gradient为共轭梯度求解函数，其输入量为矩阵A、尾向量b和迭代偏差ε，输出量为解向量x。
脚本最后提供了以Julia的LinearAlgebra包中附带的QR分解法求得的验证性解向量x2，给出了x2的值、x2与x的误差向量delta、delta的二范数‖delta‖2^2。

结构：


算法原理：


使用说明:
1.安装Julia编译环境JuliaPro-1.5.2-1。可在URL：https://juliacomputing.com/products/juliapro/，下载最新版的JuliaPro-1.5.3-1。
2.安装完毕后用Juliapro打开文件Conjugate gradient.jl，点击▶即可运行。UI右侧Workspace窗口可实时显示各变量值。
3.在脚本文件Conjugate gradient.jl中可以直接更改矩阵A、尾向量b ⃗和迭代偏差ε的值，再次点击▶即可重新计算。

附上参考书【1】给出的使用建议：
1、共轭梯度法只涉及矩阵与向量的乘法和向量内积运算，在理论上对于n阶线性方程组，最多n次迭代即可求得准确解，更像是一种直接解法；
2、对与良态方程，通常所需的迭代次数远小于矩阵阶数n；
3、对于病态问题，迭代次数有所增加，且计算精度可能受限。特别地，当迭代次数=矩阵阶数n时若仍未收敛，则应继续迭代。

源码：
#conjugate gradient Method
using LinearAlgebra
A = [
3 0 1 0 0;
0 5 0 0 0;
1 0 7 0 1;
0 0 0 5 0;
0 0 1 0 3;
]
b = [1 2 3 4 5]'
ε = 1e-6 #误差
function Conjugate_Gradient(A,b,ε)
    N = length(b) #迭代步数
    x = ones(length(b),N)*ε #声明临时解向量
    r = zeros(length(b),N) #声明临时残向量
    d = zeros(length(b),N) #声明临时d向量
    α = zeros(N) #声明临时α常数
    β = zeros(N) #声明临时β常数
    r[:,1] = b - A*x[:,1] #求解r[0]
    if r[:,1] == 0#如果初值猜对了，给出x
        x[:,N] = x[:,1]
    else #如果初值不对，共轭梯度法计算x
        d[:,1] = r[:,1] #求解d[0]
        global ∑r2 = 1
        global k = 1
        while ∑r2 >= ε
            if k == N #循环次数=矩阵阶数，退出
                break
            end
             #循环次数<矩阵阶数，继续计算
            α[k] = r[:,k]'*r[:,k]/(d[:,k]'*A*d[:,k])
            x[:,k+1] = x[:,k] + α[k]*d[:,k]
            r[:,k+1] = b - A*x[:,k+1]
            #求解r的二范数平方
            global ∑r2 = 0
            for m = 1:length(b)
                global ∑r2 = ∑r2 + r[m,k+1]^2
            end
            if (∑r2 <= ε)
                break
            end
            β[k] = r[:,k+1]'*r[:,k+1]/(r[:,k]'*r[:,k])
            d[:,k+1] = r[:,k+1] + β[k]*d[:,k]
            k = k+1
        end
    end
    return x
end
x = Conjugate_Gradient(A,b,ε)[:,k] #共轭梯度计算结果
x2 = (inv(qr(A).R)*qr(A).Q'*(b)) #QR分解计算结果
delta = zeros(length(b)) #误差向量
Σdelta = 0
for i = 1:length(b) #计算误差向量
    delta = (x - x2)/x2
    global Σdelta = Σdelta + delta^2
end
Σdelta = Σdelta^0.5 #计算误差二范数

流星雨

watermelon 发表于 2020-12-26 17:10
Python写的共轭梯度法来啦！
[code]# -*- coding: utf-8 -*-
"""

numpy确实方便！

watermelon

Python写的共轭梯度法来啦！

1. # -*- coding: utf-8 -*-
   
2. """
   
3. Created on Sat Dec 26 14:54:24 2020
   
4. @copyright: Xi'an Jiaotong University - NuTHeL
   
5. @author: Zhou Xingguang
   
6. """
   
7. 
   
8. import numpy as np
   
9. import copy
   
10. 
    
11. class ConjugateGradient(object):
    
12.     def __init__(self, matrix, b, x):
    
13.         # define the little number
    
14.         self.__ERR = 1E-6  # private parameter
    
15.         # define the MAX_STEP
    
16.         self.__MAX_STEP = 1000 # private parameter
    
17.         self.matrix = np.array(copy.deepcopy(matrix))
    
18.         self.b = np.array(copy.deepcopy(b)).reshape(np.size(b), 1)
    
19.         self.x = np.array(copy.deepcopy(x)).reshape(np.size(b), 1)
    
20.         self.last_x = self.x[:, :]
    
21.         
    
22.         self.r = np.array(np.zeros((np.size(b), 1), dtype=float))
    
23.         self.last_r = self.r[:, :]
    
24.         self.d = np.array(np.zeros((np.size(b), 1), dtype=float))
    
25.         self.last_d = self.d[:, :]
    
26.         
    
27.         self.alpha = 0
    
28.         self.beta = 0
    
29.         self.time = 0
    
30.         
    
31.    
    
32.     def CG(self):
    
33.         self.r = self.b - np.matmul(self.matrix, self.x)
    
34.         self.last_r = self.r[:, :]
    
35.         self.d = self.r[:, :]
    
36.         for time in range(self.__MAX_STEP):
    
37.             self.alpha = np.linalg.norm(self.r, ord=2)**2 / np.matmul(np.matmul(np.transpose(self.d), self.matrix), self.d)
    
38.             self.last_x = self.x[:, :]
    
39.             self.x = self.x + self.alpha * self.d
    
40.             self.last_r = self.r[:, :]
    
41.             self.r = self.b - np.matmul(self.matrix, self.x)
    
42.             self.beta = (np.linalg.norm(self.r, ord=2)**2) / (np.linalg.norm(self.last_r, ord=2)**2)
    
43.             self.d = self.r + self.beta*self.d
    
44.             self.alpha = (np.linalg.norm(self.r)**2) / np.matmul(np.matmul(np.transpose(self.d), self.matrix), self.d)
    
45.             if np.linalg.norm(self.r) <= self.__ERR:
    
46.                 print("iteration has been converged...\n")
    
47.                 self.time = time + 1
    
48.                 return self.x
    
49.         print("The result is not converged, please recheck the method")
    
50.         return None
    
51. 
    
52. def GenerateExercise(matrix, b):
    
53.     cnt = np.size(b)
    
54.     b[0, 0] = -1
    
55.     b[cnt-1, 0] = -1
    
56.    
    
57.     for i in range(cnt):
    
58.         if i == 0:
    
59.             matrix[i, i] = -2
    
60.             matrix[i, i+1] = 1
    
61.         elif i == cnt-1:
    
62.             matrix[i, i] = -2
    
63.             matrix[i, i-1] = 1
    
64.         else:
    
65.             matrix[i, i] = -2
    
66.             matrix[i, i-1] = 1
    
67.             matrix[i, i+1] = 1
    
68.     return matrix, b
    
69. 
    
70. 
    
71. if __name__ == '__main__':
    
72.     cnt = 400
    
73.     matrix = np.zeros((cnt, cnt), dtype=float)
    
74.     x = np.zeros((cnt, 1), dtype=float)
    
75.     b = np.zeros((cnt, 1), dtype=float)
    
76.    
    
77.     matrix, b = GenerateExercise(matrix, b)
    
78.    
    
79.     CG = ConjugateGradient(matrix, b, x)
    
80.     result = CG.CG()
    
81.     print(result)
    
82.     print(CG.time)

复制代码

watermelon

再来一个小弟写的HouseHolder变换的QR分解：

0xAA55

这些数学公式我看了头大。。。