无约束二次规划——共轭梯度法

二次型

二次型可以理解为初中学习的二次函数$y=a{x}^2+bx+c$的矩阵版本：
$$f(x)=\frac{1}{2}{x}^{T}Ax-b^{T}x+c$$
因为$c$对$x$的最优值没有影响，所以一般情况将其写成：
$$f(x)=\frac{1}{2}{x}^{T}Ax-b^{T}x$$
一般我们解决的是凸二次规划，这就要求$A$是一个半正定矩阵，如果$A$是一个正定矩阵，那么其就是一个严格凸二次规划，有唯一的最小值。这里举一个例子：
$$f(x)=\frac{1}{2}\left(\begin{array}{cc}{x}_1& x_2\end{array}\right)\left(\begin{array}{cc}2& 1\\ 1& 3\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)-\left(\begin{array}{cc}1& 2\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$$

当$A$是一个对称正定矩阵的时候，$f(x)$的最优解可以通过求导等于零计算出，即：
$$Ax=b$$
但是当维数很大的时候求解这个方程组在以前是个很困难的事情，所以人们想用迭代的方式去求解，也就产生了共轭梯度法。

对正定的理解

首先需要理解为什么是正定的，带有正定矩阵的函数是什么样子的。
这里我们使用MATLAB画图，做出$A$分别为正定，半正定，负定和不定矩阵的图像，函数为了简单，统一形式为：
$$f(x_1,x_2)=x^{T}Ax$$
1.正定矩阵
$$A=\left(\begin{array}{cc}2& 0\\ 0& 3\end{array}\right)$$

2.半正定矩阵
$$A=\left(\begin{array}{cc}2& 1\\ 0& 0\end{array}\right)$$

3.负定矩阵
$$A=\left(\begin{array}{cc}-2& -1\\ -1& -1\end{array}\right)$$

4.不定矩阵
$$A=\left(\begin{array}{cc}2& 0\\ 0& -2\end{array}\right)$$

记不记得高中的圆锥曲线，对于正定矩阵，其等高线是一个椭圆，当然高维空间就是椭球或者超椭球，负定就是翻下去。可以从图像中看到正定矩阵就像是一张被完全抬起来的桌布，拥有一个全局最小值。我们目前所处理的二次规划问题的函数基本都长成这个样子。

二次型的两种等高线图

1.正椭圆
第一种等高线图是比较中规中矩的，椭圆的长短轴和坐标轴平行。具体来说就是没有$x_1$和$x_2$的耦合项，表现在$A$矩阵中就是只有对角元素有值：
$$A=\left(\begin{array}{cc}3& 0\\ 0& 1\end{array}\right)$$
这里可以做如下函数的图像参考：
$$f(x)=\left(\begin{array}{cc}{x}_1& x_2\end{array}\right)\left(\begin{array}{cc}3& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)-\left(\begin{array}{cc}1& 2\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$$

2.斜椭圆
这个等高线中椭圆的长短轴和坐标轴就不是平行的了，也就是在矩阵中的其他位置有数字了：
$$A=\left(\begin{array}{cc}3& 1\\ 1& 1\end{array}\right)$$
这里可以做如下函数的图像参考：
$$f(x)=\left(\begin{array}{cc}{x}_1& x_2\end{array}\right)\left(\begin{array}{cc}3& 1\\ 1& 1\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)-\left(\begin{array}{cc}1& 2\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$$

梯度下降法的缺点

一般我们刚接触最优化，肯定会学到梯度下降法，那这个问题很显然可以用梯度下降的方法来求解，为什么还要开发其他算法呢？其实一张图就可以明白：

这里可以看到，梯度下降法需要走的路比共轭梯度要曲折的多，随着维数的增加，椭圆的倾斜等等，梯度下降寻找的路会更加复杂，而对于共轭梯度法，如果自变量有$n$个，那么它能在不多于$n$步之内找到最小值。所以梯度下降法可以用吗，当然可以，但是对于二次规划问题，有比它更快更好更直接的方法去求解。梯度下降法使用的很广泛，是一种万能算法，而共轭梯度在提出来的时候就是专门解决这类问题的，普通西装肯定没有定制的西装合身，所以可以把共轭梯度法看成一种私人定制。

坐标轮换法的思想

从正椭圆出发

搜索方向使用坐标轴的方向，这样可以轻易的看出，这种简化的共轭梯度法能够很快的找到最小点，而且有$n$个坐标轴，就能在至多$n$步之后找到最小值。

斜椭圆与特征向量

如果是一个斜椭圆,怎么处理，很简单，我们使用坐标变换将坐标系旋转到椭圆长短轴的位置不就可以了，对于这个函数：
$$f(x)=\left(\begin{array}{cc}{x}_1& x_2\end{array}\right)\left(\begin{array}{cc}3& 1\\ 1& 1\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)-\left(\begin{array}{cc}1& 2\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$$
可以先找到特征值和特征向量,将其变成：
$$f(x)=\left(\begin{array}{cc}{x}_1& x_2\end{array}\right)\left(\begin{array}{cc}0.3827& -0.9239\\ -0.9239& -0.3827\end{array}\right)\left(\begin{array}{cc}0.5858& 0\\ 0& 3.4142\end{array}\right)\left(\begin{array}{cc}0.3827& -0.9239\\ -0.9239& -0.3827\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)-\left(\begin{array}{cc}1& 2\end{array}\right)\left(\begin{array}{cc}0.3827& -0.9239\\ -0.9239& -0.3827\end{array}\right)\left(\begin{array}{cc}0.3827& -0.9239\\ -0.9239& -0.3827\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$$
设新的变量为$a=\left(\begin{array}{cc}0.3827& -0.9239\\ -0.9239& -0.3827\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$,则就可以变成一个新的二次规划问题:
$$f(a)=\left(\begin{array}{cc}{a}_1& a_2\end{array}\right)\left(\begin{array}{cc}0.5858& 0\\ 0& 3.4142\end{array}\right)\left(\begin{array}{c}{a}_1\\ a_2\end{array}\right)-\left(\begin{array}{cc}1& 2\end{array}\right)\left(\begin{array}{cc}0.3827& -0.9239\\ -0.9239& -0.3827\end{array}\right)\left(\begin{array}{c}{a}_1\\ a_2\end{array}\right)$$
图像上呈现出来的就是将坐标系旋转，然后将斜椭圆变成了正椭圆：

在各自的坐标系中的搜索方向都是$$e_1=[1,0{]}^T，e_2=[0,1{]}^T$$如果在原方向看处理过后的坐标系，那么搜索方向为：
$$e_1=[0.3827,-0.9239{]}^T,e_2=[-0.9239,-0.3827{]}^T$$

写到这里可能会有疑惑，这和共轭梯度法有什么关系?

这里的主要问题是特征向量是如何计算出来的？当特征值求解出来之后，需要求解一个线性方程组：
$$Ax=0$$
这里其实也是同样的道理，如果维数很大的话，这些特征向量找起来很麻烦，因为还是要求解一个大规模的方程组。所以共轭方向法希望找到一个向量组$S$，替代难求的特征向量组，同时也能对原矩阵进行对角化。

共轭方向法

如何理解共轭

首先看定义：
考虑正定矩阵$Q$和非零向量$d^i$和$d^j$，若满足：
$$(d^i{)}^{T}Q{d}^j=0$$
则称$d^i$和$d^j$关于Q共轭。类比于正交，共轭向量组$S=[d_0,d_1,...,d_{n-1}]$就是向量组中的向量两两共轭。

我们这里使用共轭梯度法先求解出一个方向组大概看一下。这里还是使用这个函数：
$$f(x)=\left(\begin{array}{cc}{x}_1& x_2\end{array}\right)\left(\begin{array}{cc}3& 1\\ 1& 1\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)-\left(\begin{array}{cc}1& 2\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$$
先使用共轭梯度找到最小点，然后看一看共轭方向：

共轭方向框架

对于给定的无约束二次规划问题$$minf(x)=\frac{1}{2}{x}^{T}Qx+c^{T}x$$给定初始点$x_0$以及一组关于正定矩阵$Q$的共轭方向$d_0,d_1,...,d_{n-1}$，这样我们就可以进行如下的迭代：
$$x_{k+1}=x_k+{\alpha }_k{d}_k$$
这里的步长选择很简单，因为是一个凸的，所以有：
$$\alpha =argmi{n}_{\alpha }f(x_k+\alpha d_k)$$
这样求导等于0，所以有：
$${\alpha }_k=-\frac{\nabla f(x_k)d_k}{d_k^{T}Q{d}_k}$$
这里最主要的问题就是共轭方向$d_0,d_1,...,d_{n-1}$的选取，所以共轭方向法是一大类方法，而共轭梯度法是在构造共轭方向的时候采用梯度信息。

共轭向量代替特征向量实现对角化

那么我们先不管如何计算方向，默认方向已经给出了，来看看能否将Q对角化呢？类似于特征向量，我们将共轭向量组进行如下操作：
$$\left(\begin{array}{c}{d}_0^T\\ d_1^T\\ ⋮\\ d_{n-1}^T\end{array}\right)Q\left(\begin{array}{cccc}{d}_0& d_1& \dots & d_{n-1}\end{array}\right)$$
下标不同的因为共轭的性质都变成了0.所以最后计算出来必然为对角矩阵，也即实现了对$Q$矩阵的对角化：
$$\left(\begin{array}{cccc}{d}_0^{T}Q{d}_0& 0& \dots & 0\\ 0& d_1^{T}Q{d}_1& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & d_{n-1}^{T}Q{d}_{n-1}\end{array}\right)$$
那类比于特征向量，这个时候的求解问题就变了，将原变量替换$x=S\hat{x}$,这样原问题变成：
$$minf(x)=\frac{1}{2}{\hat{x}}^T{S}^{T}QS\hat{x}+c^{T}S\hat{x}$$
这样就达到了我们想要的效果，二次项对角，完全变成了一个和坐标轴平行的正椭圆。这里还是看一下图像：

原图像函数为：
$$f(x)=\left(\begin{array}{cc}{x}_1& x_2\end{array}\right)\left(\begin{array}{cc}3& 1\\ 1& 1\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)-\left(\begin{array}{cc}1& 2\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$$
使用共轭向量组变换之后的函数为：
$$f(x)=\left(\begin{array}{cc}{x}_1& x_2\end{array}\right)\left(\begin{array}{cc}0.5682& 0\\ 0& 0.5568\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)-\left(\begin{array}{cc}1.1364& 1.1136\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$$
那搜索方向在另一个空间是如何进行的呢？
我们在这个空间的搜索是按照$$x^{\ast }=x_0+{\alpha }_0{d}_0+{\alpha }_1{d}_1+{\alpha }_2{d}_2+...+{\alpha }_{n-1}{d}_{n-1}$$
进行的，由于有$\hat{x}=S^{-1}x$,所以能得到：
$$S^{-1}{x}^{\ast }=S^{-1}{x}_0+{\alpha }_0{S}^{-1}{d}_0+{\alpha }_1{S}^{-1}{d}_1+{\alpha }_2{S}^{-1}{d}_2+...+{\alpha }_{n-1}{S}^{-1}{d}_{n-1}$$
这里需要注意:
$$S^{-1}{d}_i={\left(\begin{array}{cccc}{d}_0^T& d_1^T& \dots & d_{n-1}^T\end{array}\right)}^{-1}{d}_i=\left(\begin{array}{cccc}0& 0& \dots 1& \dots 0\end{array}\right)=e_{i+1}$$
所提最优解的寻找过程变成：
$$S^{-1}{x}^{\ast }=S^{-1}{x}_0+{\alpha }_0{e}_1+{\alpha }_1{e}_2+{\alpha }_2{e}_3+...+{\alpha }_{n-1}{e}_n$$
这说明，在另一个空间是按照坐标轴搜索的，和特征向量的效果是一样的。从上图也可以明白，右边蓝色和橙色的搜索方向的向量通过空间变换，变到了左边的x轴和y轴

一个重要的性质 $\nabla f(x_k{)}^T{d}_i=0$

因为这个公式很重要，所以需要单独拿出来说一下。
首先$\nabla f(x_k{)}^T{d}_{k-1}=0$是根据最优步长计算出来的条件，所以这里我们把其当作默认条件，这里可以简单的理解为，一个人要下山，他选好了方向，走多少步呢，他选择走到某一个点，这个点能让他转九十度就可以得到速度下降最快的方向。那剩余的方向为什么和$\nabla f(x_k{)}^T$垂直呢：
$$\begin{array}{rl}\nabla f(x_k{)}^T{d}_i& =(Q{x}_k+c{)}^T{d}_i\\ & =(Q(x_{i+1}+{\alpha }_{i+1}{d}_{i+1}+...+{\alpha }_{k-1}{d}_{k-1})+c{)}^T{d}_i\\ & =(Q{x}_{i+1}+c{)}^T{d}_i\\ & =\nabla f(x_{i+1}{)}^T{d}_i=0\end{array}$$

共轭梯度法

下面我们着重解决如何构造共轭方向。共轭梯度不仅是使用梯度信息，它还有一个特点是一边迭代，一边生成共轭向量，并不像特征向量一样可以一次性直接生成特征向量组。这里先给出伪代码框架：

• 给定初始点，令$d_0=-\nabla f(x_0)$,给定一个误差限$e$
• 判断$||\nabla f(x_0)||\le e$是否成立，是则终止
• 计算步长：${\alpha }_k=-\frac{\nabla f(x_k)d_k}{d_k^{T}Q{d}_k}$
• $x_{k+1}=x_k+{\alpha }_k{d}_k$
• 计算方向$d_{k+1}=-\nabla f(x_{k+1})+somedirection$
• $k=k+1$，to Step 1
  这里我们可以顺便看一下共轭梯度的搜索路径：
  

回到公式这里主要关注倒数第二项，以下一时刻点的梯度方向为基础，再进行一些修正，希望下一个方向与前面所有方向都共轭。这让我们想到了什么？因为共轭可以看成另一个空间的正交，而正交的情况下想做这件事就是施密特正交化的思想，这里思想一样。

快速回顾施密特正交化

给定一个向量$v$和一个初始方向$\alpha$，如何构造一组正交基呢？

首先需要清楚一个向量到另一个向量的投影是如何用数学公式表示的也就是图中红色向量是怎样通过橘色的向量$v$映射到$\alpha$方向上的。
$$v_{red}=\frac{<v,\alpha >}{<\alpha ,\alpha >}\alpha$$
其中$<>$表示内积，也就是点乘。其实就是先映射再给单位向量的过程，映射体现在$$\frac{v\alpha cos\theta }{\alpha }$$给出单位向量表现在$$\frac{\alpha }{||\alpha ||}$$。
那么第一个方向设置为${\beta }_1=\alpha$，下一个方向自然可以设置成蓝色虚线的方向，根据向量的求和，可以计算出$${\beta }_2=v-v_{red}=v-\frac{<v,{\beta }_1>}{<{\beta }_1,{\beta }_1>}{\beta }_1$$
如果是三维，那么公式就要再加一项：
$${\beta }_3=v-\frac{<v,{\beta }_1>}{<{\beta }_1,{\beta }_1>}{\beta }_1-\frac{<v,{\beta }_2>}{<{\beta }_2,{\beta }_2>}{\beta }_2$$
以此类推

从正交化到共轭化

那么回到共轭我们应该怎么处理呢？这里我们的初始方向是当前位置的梯度的负方向$-\nabla f(x_0)$，这里我们做一个类比，在正交化的过程中，我们的投影是两个向量做点集，如果两个向量正交的话，那么会满足$d_i\ast d_j=0$，这里共轭的映射操作可以看作$d_{i}Q{d}_j$，如果满足共轭，那么会有$d_{i}Q{d}_j=0$,那么按照上述的类比，下一时刻的共轭方向应该是
$$d_{k+1}=-\nabla f(x_{k+1})+\sum _{i=0}^k\frac{\nabla f(x_{k+1})Q{d}_i}{d_i^{T}Q{d}_i}{d}_i$$
也可以换一种思路求出这个公式。根据定义希望下一个方向与之前的方向都共轭，就会要满足$d_{k+1}^{T}Q{d}_i=0,i=0,1,2,...,k$，类比于施密特正交的思想，我们有：
$$(-\nabla f(x_{k+1})+{\beta }_0{d}_0+{\beta }_1{d}_1+...+{\beta }_k{d}_k{)}^{T}Q{d}_i=0,i=0,1,2,...,k$$
因为其他方向都和其满足共轭，所以只剩下第：
$$-\nabla f(x_{k+1}{)}^{T}Q{d}_i+{\beta }_i{d}_i^{T}Q{d}_i=0,i=0,1,2,...,k$$
这样就得到了${\beta }_i$为：
$${\beta }_i=\frac{\nabla f(x_{k+1}{)}^{T}Q{d}_i}{d_i^{T}Q{d}_i},i=0,1,2,...,k$$
我们用图片看一下这里共轭和正交的关系：

左边是变换空间之后的搜索方向，右边是在我们这个空间的视角之下的搜索方向。

一些化简

其实写到这里，共轭梯度法基本介绍完了，但是还需要一些数学形式上的化简才能到最终版本，目的就是能使用梯度表示的就用梯度表示，最终目的是将其用到非线性当中。

求和化简

首先我们先看一下每个共轭放向前面的这些项$\nabla f(x_{k+1}{)}^{T}Q{d}_i$,首先要清楚一点，我们处理的是二次规划，所以有一个关系是:
$$\nabla f(x)=Qx+c$$
同时还有一个步长关系是：
$$x_{k+1}=x_k+{\alpha }_k{d}_k$$
现在对分子进行一些处理：
$$\begin{array}{rl}\nabla f(x_{k+1}{)}^{T}Q{d}_i& =\nabla f(x_{k+1}{)}^{T}Q(x_{i+1}-x_i)\frac{1}{{\alpha }_i}\\ & =\frac{1}{{\alpha }_i}\nabla f(x_{k+1}{)}^T(\nabla f(x_{i+1})-\nabla f(x_i))\\ & =\frac{1}{{\alpha }_i}\nabla f(x_{k+1}{)}^T(({\beta }_0{d}_0+{\beta }_1{d}_1+...+{\beta }_i{d}_i-{\beta }_{i+1}{d}_{i+1})-({\beta }_0{d}_0+{\beta }_1{d}_1+...+{\beta }_{i-1}{d}_{i-1}-{\beta }_i{d}_i))\\ & =\frac{1}{{\alpha }_i}\nabla f(x_{k+1}{)}^T(2{\beta }_i{d}_i-{\beta }_{i+1}{d}_{i+1})\end{array}$$
这里根据前面的重要的性质可以看出，当$i=0,1,...,k-1$的时候$\nabla f(x_{k+1}{)}^{T}Q{d}_i=0$所以求和项中除最后一项，剩下的都变成了零(其实这里还可以得到一个性质就是$\nabla f(x_{k+1}{)}^T\nabla f(x_i)=0$)。这样我们就能得到一个简单的公式，即：
$$\begin{array}{rl}{d}_{k+1}& =-\nabla f(x_{k+1})+\sum _{i=0}^k\frac{\nabla f(x_{k+1})Q{d}_i}{d_i^{T}Q{d}_i}{d}_i\\ & =-\nabla f(x_{k+1})+\frac{\nabla f(x_{k+1})Q{d}_k}{d_k^{T}Q{d}_k}{d}_k\end{array}$$
之后的两个公式简化就很简单了

步长公式简化

原来的步长公式为：
$$\alpha =-\frac{\nabla f(x_k)d_k}{d_k^{T}Q{d}_k}$$
将$d_k$替换掉：
$$\begin{array}{rl}\alpha & =-\frac{\nabla f(x_k)d_k}{d_k^{T}Q{d}_k}\\ & =-\frac{\nabla f(x_k)(-\nabla f(x_k)+{\beta }_{k-1}{d}_{k-1})}{d_k^{T}Q{d}_k}\\ & =\frac{\nabla f(x_k{)}^T\nabla f(x_k)}{d_k^{T}Q{d}_k}\end{array}$$
可以看到这些化简频繁的用到了之前说的重要性质。

系数简化

根据前面的推导，我们能得到
$$\nabla f(x_{k+1}{)}^{T}Q{d}_k=\frac{1}{{\alpha }_i}\nabla f(x_{k+1}{)}^T(\nabla f(x_{k+1})-\nabla f(x_k))$$
这样我们将其带入到系数中得到：
$${\beta }_k=\frac{1}{{\alpha }_i}\frac{\nabla f(x_{k+1}{)}^T(\nabla f(x_{k+1})-\nabla f(x_k))}{d_k^{T}Q{d}_k}$$
将刚刚简化好的步长带进去：
$$\begin{array}{rl}{\beta }_k& =\frac{\nabla f(x_{k+1}{)}^T(\nabla f(x_{k+1})-\nabla f(x_k))}{\nabla f(x_k{)}^T\nabla f(x_k)}\\ & =\frac{\nabla f(x_{k+1}{)}^T\nabla f(x_{k+1})}{\nabla f(x_k{)}^T\nabla f(x_k)}\end{array}$$
其实不化简也是可以的，但是化简之后有利于我们从线性方法转变到非线性方法。观察化简之后的式子能够发现，这里面的系数其实是没有任何二次项的系数存在的(矩阵Q)，而只有梯度，下面我们简单看看共轭梯度法是如何转变到非线性的。

非线性共轭梯度法(PRP,FR)

秉持着用到再了解的思想，这里只给出框架，具体原理等以后需要了再探究，不过这些公式都眼熟，前面都推导过：

• 给定初始点，令$d_0=-\nabla f(x_0)$,给定一个误差限$e$
• 判断$||\nabla f(x_0)||\le e$是否成立，是则终止
• 利用线性搜索计算步长${\alpha }_k$
• $x_{k+1}=x_k+{\alpha }_k{d}_k$
• 计算方向$d_{k+1}=-\nabla f(x_{k+1})+{\beta }_k{d}_k$
• $k=k+1$，to Step 1

这里面有两种常见的${\beta }_k$算法：
第一种是PRP算法
$${\beta }_k=\frac{\nabla f(x_{k+1}{)}^T(\nabla f(x_{k+1})-\nabla f(x_k))}{\nabla f(x_k{)}^T\nabla f(x_k)}$$
第二种是FR算法：
$${\beta }_k=\frac{\nabla f(x_{k+1}{)}^T\nabla f(x_{k+1})}{\nabla f(x_k{)}^T\nabla f(x_k)}$$

附录：二次规划的MATLAB代码

这里参考了MATLAB - 共轭梯度法
注意，这里的格式中的$A$和$b$对应如下函数：
$$f(x)=\frac{1}{2}{x}^{T}Ax-b^{T}x$$
更改相应的$A$和$b$即可。同时MATLAB对于二次规划问题有quadprog求解器，使用起来非常方便，具体请见：quadprog

function func()
A = [6 2;
    2 2];
b=[1 2]';
x0=[0 0]';
max_iter=10000;
 
fprintf('\n');
fprintf('共轭梯度法:\n');
fprintf('==========================\n');
[y,iter]=cgm (A,b,x0,max_iter);
fprintf('\n');
fprintf('迭代次数:\n   %d \n',iter);
fprintf('方程的解: \n');
fprintf('%10.6f',y);
fprintf('\n\n=========================\n\n');
end
 
function [x,iter] = cgm (A,b,x0,max_iter)
x=x0;
epsilon=1.0e-6;
fprintf('\n x0= ');
fprintf('   %10.6f',x0);
r=b-A*x;
d=r;
for k=0:max_iter
    alpha=(r'*r)/(d'*A*d);
    xx=x+alpha*d;
    rr=b-A*xx;
    if (norm(rr,2)/norm(b,2))<= epsilon
        fprintf('\n 找到啦~');
        iter = k+1;
        x=xx;
        r=rr;
        fprintf('\n x%d = ',k+1);
        fprintf('   %10.6f',x);
        fprintf('\n r%d = ',k+1);
        fprintf('   %10.6f',r);
        
        return
    end
    beta=(rr'*rr)/(r'*r);
    d=rr+beta*d;
    x=xx;
    r=rr;
    fprintf('\n x%d = ',k+1);
    fprintf('   %10.6f',x);
end
iter = max_iter;
return 
end