（一）共轭梯度算法

共轭梯度算法是干什么？

共轭梯度算法是一种迭代算法，在一次次的迭代中最终求得结果，可以类比牛顿迭代法。共轭梯度算法主要用在求解矩阵方程，也就是求解n元一次方程组,如Ax=b的解x。比一般的迭代算法都要快，最多迭代n次就能出结果。

基本思想

共轭梯度法是属于最小化类的迭代方法。为了求解Ax = b这样的一次函数，可以先构造一个二次齐函数
$$f(x)=\frac{1}{2}{x}^{T}Ax-b^{T}x$$
这样求解Ax = b的值可以转换为求解f(x)的最小值。（证明略，且矩阵A要对称且正定）。
为什么一般都要求矩阵A对阵正定？

1. 实际上正定是必须的,正定可以保证f(x)有最小值，也就是Ax=b有且只有一个解。
2. 但其实并不需要矩阵一定是对称的，要求对称可能是好构造？只要主对角线的值全部不为0且对称的矩阵就是一定正定了。可能是这个原因。

关于对称和正定参考知乎解答

共轭梯度算法的步骤

1. 初始化,x(k)表示第k次迭代的解向量。d(k)表示第k次迭代的方向向量。r(k)表示第k次迭代的残差向量。

$$x(0)=0,d(0)=0,r(0)=-b_n$$

进行第k次迭代，主要分为四个步骤

1. 计算残差向量
   $$r(k)=Ax(k-1)-b$$

2. 计算方向向量
   $$d(k)=-r(k)+\frac{r^{\mathrm{T}}(k)r(k)}{r^{\mathrm{T}}(k-1)r(k-1)}d(k-1)$$

3. 计算步长
   $$\alpha (k)=-\frac{d^{\mathrm{T}}(k)r(k)}{d^{\mathrm{T}}(k)Ad(k)}$$

4. 更新解向量
   $$x(k)=x(k-1)+\alpha (k)d(k)$$

伪代码

输入：A(n*n矩阵)，b(n*1)矩阵
输出：x(n*1)矩阵

x = [0,0...0],d = [0,0...0],r = -b;
for (i = 1:n){     //n次迭代
    denom1 = r^T*r //计算r的转置乘以r
    r = A*x-b //计算残差r = Ax-b
    num1 = r^T*r; //计算r的转置乘以r
    if(num1 < 精度)
        break;
    //计算方向向量d
    d = -r + num1/denom1*d;//d和r都是向量
    
    num2 = d^T*r;//计算d的转置乘以r
    denom2 = d^T*(A*d) //计算d的转置乘以A乘以d最后
    a = -num2/denom2; //计算步长
    x = x+Ad // 修正x
    
}
print(x);

主要有以下八步

1. 计算denom1：$$denom1=r^T\ast r$$

2. 计算残差：$$r=A\ast x-b$$

3. 计算内积：$$num1=r^T\ast r$$

4. 计算方向向量：$$d=-r+\frac{num1}{denom1}d$$

5. 计算内积：$$num2=d^T\ast r$$

6. 计算denom2：$$denom2=d^T\ast A\ast d$$

7. 计算步长：$$length=\frac{-num2}{denom2}$$

8. 修正x：$$x=x+length\ast d$$

c++ 代码

#include<iostream>
#include<vector>
#include<cmath>
#include<stdio.h>
using namespace std;
int N =10;
//初始化A为一个N*N的矩阵的对称矩阵
vector<vector<double> > A(N,vector<double>(N,0));

//初始化数组b
vector<double> b(N,1);

 //初始化残差r,结果x,计算方向向量d
vector<double> r(N,-1);
vector<double> d(N,0);
vector<double> x(N,0);

void displayA(vector<vector<double> > &a){
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){
            cout<<a[i][j]<<"    ";
        }
        cout<<endl;
    }
}

void displayb(vector<double> &b){
    for(int i=0;i<b.size();i++){
        cout<<b[i]<<" ";
    }
    cout<<endl;
}

//计算内积，也就是模的平方
double INNER_PRODUCT(vector<double> &a,vector<double>&b){
    double res = 0;
    for(int i=0;i<N;i++){
        res+=a[i]*b[i];
    }
    return res;
}

//更新残差 r = A*x-b
vector<double>&  MATRIX_VECTOR_PRODUCT(vector<vector<double> > &a,vector<double> &x,vector<double>& b){
    double temp = 0;
    for(int i=0;i<N;i++){
        temp = 0;
        for(int j=0;j<N;j++){
            temp += a[i][j]*x[j];
        }
        r[i] = temp - b[i];
    }
    return r;
}

//计算dtAd
double MATRIX_PRODUCT(vector<vector<double> >&a,vector<double> &d){
    double res = 0;
    double temp = 0;
    for(int i=0;i<N;i++){
        temp = 0;
        for(int j = 0;j<N;j++){
            temp += d[j]*A[i][j];
        }
        res += temp*d[i];
    }
    return res;
}


int main(){

    //初始化A
    for(int i=0;i<N;i++){
        for(int j =0;j<N;j++){
            if(i==j){
                A[i][j] = 2;
            }
            if(abs(i-j) == 1){
                A[i][j] = -1;
            }
        }
    }
    //displayA(A);
    //displayb(b);
    //开始迭代
    int count = 0;
    for(int i =0;i<N;i++){
        count++;

        //计算r^Tr,
        double denom1 = INNER_PRODUCT(r,r);
        r = MATRIX_VECTOR_PRODUCT(A,x,b);
    
        double num1 = INNER_PRODUCT(r,r);
        if(num1 < 0.000001){
            break;
        }
        double temp = num1/denom1;
        //计算方向向量d
        for(int j = 0;j<N;j++){
            d[j] = -r[j]+temp*d[j];
        }
        double num2 = INNER_PRODUCT(d,r);
        double denom2 = MATRIX_PRODUCT(A,d);

        //计算步长
        double  length = -num2/denom2;

        //修正x
        for(int j=0;j<N;j++){
            x[j] = x[j]+ length*d[j];
        }
    }

    cout<<"迭代次数: "<<count<<endl;
    displayb(x);
    return 0;
}

直接编译执行，输出迭代的次数和方程的解。
共轭梯度算法的主要参考