精华内容
下载资源
问答
  • 共轭梯度

    2019-04-25 22:45:48
    共轭梯度 标签(空格分隔): 数值优化 线性代数 这篇文章最早是在作业部落写的,但是没什么人访问,所以这里补录一下 本文介绍和总结共轭梯度(conjugate gradient)算法的如下内容: 共轭梯度算法的背景和由来 ...

    共轭梯度

    标签(空格分隔): 数值优化 线性代数


    这篇文章最早是在作业部落写的,但是没什么人访问,所以这里补录一下

    本文介绍和总结共轭梯度(conjugate gradient)算法的如下内容:

    共轭梯度算法的背景和由来
    梯度算法的几何理解和共轭梯度算法的流程
    基于SSOR预条件子的共轭梯度算法
    非线性共轭梯度算法

    背景

    共轭梯度算法的由来,源于求解如下所示的二次型问题:
    KaTeX parse error: \tag works only in display equations
    如果对于问题的求解,可以被转化为如上所示的二次型表达式(Quadratic Form)的极值,那么求解这个问题的方法也称为二次规划(Quadratic Programming)
    二次规划问题,有一个极好的性质,这个性质直接导致了该问题的求解有较好的几何解释,这个性质就是如果矩阵AA为对称矩阵,那么把这个函数对列向量xx求导取极值,得到的梯度方向向量xx满足:
    KaTeX parse error: \tag works only in display equations
    这个Ax=bAx=b是经典的线性代数求解线性方程组的问题,求极值即是求解线性方程组,而矩阵AA的性质,直接决定了梯度方向的情况。
    对一个从nn维到11维的函数映射f(x)f(x),对列向量求导即:
    (3)f(x)=[f(x)x1f(x)xn]f^{'}(x)=\begin{bmatrix} \frac{\partial f(x)}{\partial x_{1}} \\ \vdots\\ \frac{\partial f(x)}{\partial x_{n}} \end{bmatrix}\tag{3}
    如果仅仅只看二次型中的二次项,那么
    (4)g(x)=xTAx=i=0nj=0naijxixj=...+i=0naikxixk+j=0nakjxkxjakkxkxk+...=...+i=0,iknaikxixk+j=0,jknakjxkxj+akkxkxk+...    g(x)xk=i=0naikxi+j=1nakjxj g(x)=x^TAx=\sum_{i=0}^{n}\sum_{j=0}^{n}a_{ij}x_ix_j \\ =... + \sum_{i=0}^{n}a_{ik}x_ix_k+\sum_{j=0}^{n}a_{kj}x_kx_j-a_{kk}x_kx_k + ...\\ =...+\sum_{i=0,i\ne k}^{n}a_{ik}x_ix_k + \sum_{j=0,j\ne k}^{n}a_{kj}x_kx_j+a_{kk}x_kx_k+...\\ \implies \frac{\partial g(x)}{\partial x_k}=\sum_{i=0}^{n}a_{ik}x_i+\sum_{j=1}^{n}a_{kj}x_j \tag{4}
    仔细观察上述表达式,可以发现,求导结果的第kk个元素等于矩阵AA的第kk行和第kk列分别与向量点积的结果,所以
    (5)g(x)x=ATx+Ax    f(x)x=12(AT+A)xb\frac{\partial g(x)}{\partial x}=A^Tx+Ax \implies \frac{\partial f(x)}{\partial x}=\frac{1}{2}(A^T+A)x-b\tag{5}
    如果A是对称矩阵的形式,那么,通过求导求解极值相当于求解方程组
    (6)Ax=bAx=b \tag{6}
    二次规划是一些问题的最简单的建模表达形式,比如,在EDA(电子电气自动化)软件开发中,有这样的需求,就是要对各个电子元器件在电路板或者FPGA的电路平面进行自动化布局。每个元器件可以抽象为一些有面积的方型块,然后通过设定一些约束,定义损失函数,最后,通过不断迭代减少损失函数,更新各个电路元件的坐标信息,最终得到各个电子元器件在电路中的较优的全局布局结果。这些图可以有多层,而且损失函数通常也是非线性的,迭代求解通常基于各种梯度下降算法。期间需要考虑各个图层每一块面积的覆盖率,重叠程度等等参数。
    如果对问题进行全面的建模,那么自然会比较复杂,如果只考虑最简单的情况,通过把元器件视为一个点,同时通过引入锚点,那么这个问题可以被简化为二次规划问题,通常二次规划问题用于电子元器件的预布局,然后才引入非线性函数进一步处理。
    如果优化的目标可以被简化为让各个元器件之间的xxyy坐标之间的距离尽可能分散,那么此时损失函数等于
    (7)cf=i,jwij((xixj)2+(yiyj)2)+f,H(f)wfH(f)((xfxH(f))2+(yfyH(f))2)cf=\sum_{i,j}w_{ij}((x_i-x_j)^2+(y_i-y_j)^2)+ \\ \sum_{f,H(f)}w_{fH(f)}((x_f - x_{H(f)})^2+(y_f-y_{H(f)})^2) \tag{7}
    如果令x=(x1,x2,...,xn)Tx=(x_1,x_2,...,x_n)^Ty=(y1,y2,...,yn)Ty=(y_1,y_2,...,y_n)^T,那么,完全可以表达成两个二次型表达式,其中后半部分的wfH(f)w_{fH(f)}表示的是锚点到可动点的权重,实际优化该函数,锚点部分只会对AA的对角线,以及bTb^T有影响,优化的结果就是在锚点的约束下,各个点的距离之间尽可能小。
    由于问题的特殊形式,问题的求解完全也可以把xxyy左边分别分成两个部分:
    (8)Ax=bxAx=b_x \tag{8}(9)Ay=byAy=b_y \tag{9}
    也就是等价于求解线性方程组了

    几何理解和共轭梯度算法的流程

    共轭梯度算法来源于二次型问题的求解,二次型问题可以在几何上直观的理解,即使是最速梯度下降(Steepest Descent),也可以这样来理解。

    二次型的图形

    二次型的几何图解完全取决于矩阵AA,如果AA是正定的矩阵,那么图像将会是:

    也就是是一个碗状朝上的结构,因为对于正定矩阵来说,任何非零向量,都有:
    (10)xTAx>0x^TAx>0 \tag{10}
    这个定义的本质是,所有的特征向量所对应的特征值都是正的,因为:
    (11)xTAx=xTλx>0    λ>0x^TAx=x^T\lambda x>0 \implies \lambda > 0 \tag{11}
    对于二次型来说,求导的梯度方向由Ax=bAx=b决定,AA正定(positive definite),意味着各个向量xx变换到AA的列空间中得到的向量AxAx都是单调递增的,所以二次型表达式(如果是二维的)的几何图形就是中间凹陷,四周往上的碗状。同理,可以定义半正定,负定等,那么图形就是:

    对于既非负定也非正定的情况,就是马鞍型,这个不管使用什么梯度方法,理论上都可能不会收敛。对于更高维度的情况,其实也是类似的,这里边正定或者负定的属性,决定了梯度方向是不是单调递增/递减的。对二次型,后边都是基于对称正定矩阵来分析的,那些梯度方法也源于此来进行基本的几何分析。

    Steepest Descent

    对于最速梯度下降(Steepest Descent)而言,选定了初始的方向向量x0x_0后,后续每一次迭代所选择的方向向量由当前的梯度方向决定:
    (12)xi+1=xi+αrix_{i+1}=x_i+\alpha r_i \tag{12}
    其中rir_i是方向向量,默认就是取梯度的负方向。这个过程的本质,其实就是线搜索(line search),如下图所示:

    上图是二次型在变量维度为二维的时候的共势等高图,线是梯度负方向,每次迭代的本质就是在沿着这条线走,直到到达某点后停止,得到此时的向量xi+1x_{i+1},一般情况下,步长α\alpha是可以固定的,特别是在函数是非线性,且非常复杂的时候,要想每一次都计算最好的步长,可能是非常复杂且不划算的,不过对于二次型问题而言,我们却可以设法找到最优的α\alpha,通过查看上图的各个结果向量,为了更接近与最小点,ri+1r_{i+1}应该和rir_i正交,此时,才是沿着当前的搜索方向的最优解。因为这个时候,向量已经沿着该方向尽可能往这个方向走,且已经不可能贡献再多的效用了,再走就又回去了。
    对于二次型,梯度负方向等于:
    (13)fi(x)=Axib    ri=fi(x)=bAxif^{'}_i(x)=Ax_i-b \implies r_i=-f^{'}_i(x)=b-Ax_i \tag{13}
    利用正交关系进行推导,可以得到:
    ri+1Tri=(bAxi+1)Tri=(bA(xi+αri))Tri=0(bAxi)TriαriTATri=riTαriTAri=0 r_{i+1}^Tr_i=(b-Ax_{i+1})^Tr_i=(b-A(x_i+\alpha r_i))^Tr_i=0 \\ (b-Ax_i)^Tr_i-\alpha r_i^TA^Tr_i=r^T_i-\alpha r_i^TAr_i=0
    所以
    (14)α=riTririTAri\alpha=\frac{r_i^Tr_i}{r_i^TAr_i} \tag{14}
    这里边每次给出rir_i的时候都相当于直接求导了,如果根据迭代的形式给出下一个方向向量,那么
    (15)ri+1=bAxi+1=bA(xi+αri)=riαArir_{i+1}=b-Ax_{i+1}=b-A(x_i+\alpha r_i)=r_i-\alpha Ar_i \tag{15}
    用迭代的形式给出每一次的梯度方向向量,虽然也许能够节省计算开销,但是限于实际机器的浮点精度,可能会产生累积误差,所以一般也需要在迭代了几十次后,重新严格计算梯度,以避免传递误差过大。
    这样我们就可以得到最速梯度下降算法的流程是:

    import numpy as np
    
    
    def steepest_descent(A, b, x_initial, max_step, thresh=0.00001):
        assert(isinstance(A, np.matrix))
        assert (isinstance(b, np.matrix))
        assert (isinstance(x_initial, np.matrix))
        x = x_initial
        for _ in range(max_step):
            r = b - A * x
            alpha = (r.transpose() * r) / (r.transpose() * A * r)
            x = x + r * alpha
            dist = np.sqrt(np.sum(np.square(b - A * x)))
            if dist < thresh:
                break
        return x
    
    
    if __name__ == '__main__':
        N = 100
        Ar = np.mat(np.random.rand(N, N))
        As = Ar * Ar.transpose()  # get positive definite matrix
        bn = np.mat(np.random.rand(N, 1))
        xi = np.mat(np.zeros((N, 1)))
        xr = steepest_descent(As, bn, xi, 1000)
        print('1000:', np.sqrt(np.sum(np.square(bn - As * xr))))
        xr = steepest_descent(As, bn, xi, 10000)
        print('10000:', np.sqrt(np.sum(np.square(bn - As * xr))))
        xr = steepest_descent(As, bn, xi, 20000)
        print('20000:', np.sqrt(np.sum(np.square(bn - As * xr))))
    
    

    上边的Python代码体现了SD的最直接流程,梯度都是直接计算的,实际真的要应用,还要考虑稀疏矩阵,计算效率问题,否则效率太低了,也没有使用价值。当然,这个只用于二次型问题。
    对一些非线性函数问题,就需要通过别的方式求导,而且步长α\alpha也需要确定,或者设定个固定值!如下所示是该代码的一次运行结果:

    1000: 1.9453821317304896
    10000: 0.9116628408830175
    20000: 0.6950563935826073
    

    Conjugate Gradient

    运行SD的例子,会发现,随机生成的一个对称矩阵,迭代了几万次仍然不收敛,可以见得,按照SD的理论公式求解Ax=bAx=b效果有些不理想,收敛过慢。在数值优化里边,一个矩阵迭代求解方法是否收敛,取决于这个矩阵是否严格对角占优,即
    (16)A=(aij),i[1,n],aii&gt;ijaijA=(a_ij), \forall i \in [1, n], |a_{ii}| \gt \sum_{i \ne j}|a_{ij}| \tag{16}
    越是严格对角占优,越是收敛得越快。
    SD这个方法,使用梯度方向,每次选择梯度方向都和上一次的正交,这样迭代的结果取决于初始点的选择和AA,如果初始点不理想,可能一直在锯齿状的接近极值点,但是迭代多次后仍然不能命中。
    对于当前的二次型而言,不是多么复杂的问题,变量的维度也许就是nn而已,存不存在这样的迭代方法,使得每一次迭代,都消减一个维度方向,后边的方向和前边的方向都完全正交呢?这样的话,我只需要nn步迭代就可以直接命中极值点了。这个思想,我觉得就是CG背后的想法。如果选择nn个线性无关的向量u1,u2,...,unu_1,u_2,...,u_n,那么依据施密特(Gram-Schmidt Process)正交化的过程,每个方向向量为:
    (17)di=ui+k=0i1βikdkd_i=u_i+\sum_{k=0}^{i-1}\beta_{ik}d_k \tag{17}

    Conjugate Direction

    在CG中,引入了共轭方向的概念。这个概念本质上就是对正交的推广,任何向量在欧几里得空间内正交意味着:
    (18)xTy=0x^Ty=0 \tag{18}
    如果引入矩阵AA,且
    (19)xTAy=0x^TAy=0 \tag{19}
    那么称向量xxyyAA正交的,这里边AA都是方阵,看起来就像是把变量线变换到AA的行列空间中,然后这些向量在AA的行列空间中正交。这些两两AA正交的向量被称为共轭方向向量。如果给个图解,那么就会如下图所示:

    前后两个方向向量是AA正交(A-orthogonal)的,意味着diTAdj=0d_i^TAd_j=0,即使这两个向量在当前空间不正交,在AA的空间中正交也行,可以想象,当前空间中的所有向量都变换到AA的空间后,着两个向量在AA空间中就是正交的。

    二次型的标准化

    上边的过程容易让人联想起线性代数中的二次型标准化过程,二次型标准化,首先通过平移变换,让二次型变成了如下的形式:
    (20)f(x)=xTAxf(x)=x^TAx \tag{20}
    然后找寻AAnn个标准正交的特征向量C=(c1,c2,...,cn)C=(c_1,c_2,...,c_n),然后做如下变换:
    (21)Cy=x&ThickSpace;&ThickSpace;f(x)=xTAx=yTCTACy=yTΛyCy=x \implies f(x)=x^TAx=y^TC^TACy=y^T\Lambda y \tag{21}
    最终AA变成了对角矩阵的形式,这个时候对应的二次型的几何图形就是很标准的碗状了
    个人觉得,AA正交也是希望在AA标准化后的空间中正交。其实标准化前后,所代表的线性空间是不变的。在几何理解里边,都以标准化后的空间中的向量来分析理解。

    CG的推导过程

    二次型的基本迭代流程还是不变的:
    (22)xk+1=xk+αkdkx_{k+1}=x_{k}+\alpha_{k}d_k \tag{22}
    如果采用迭代求解梯度,那么下边的等式仍然成立:
    (23)rk+1=rkαkArkr_{k+1}=r_k-\alpha_{k}Ar_k \tag{23}
    对等式22向αk\alpha_{k}求导:
    (24)f(xk+1)αk=f(xk+1)Txk+1αk=rk+1Tdk=(AxAxk+1)Tdk=dkTAek+1=0&ThickSpace;&ThickSpace;dkTAek+1=0 \frac{\partial f(x_{k+1})}{\partial \alpha_{k}}=f^{&#x27;}(x_{k+1})^T\frac{\partial x_{k+1}}{\partial \alpha_k} \\ =-r_{k+1}^Td_k=-(Ax-Ax_{k+1})^Td_k=-d_k^TAe_{k+1}=0 \\ \implies d_k^TAe_{k+1}=0 \tag{24}
    可以看到,最优的步长是使得下一个误差方向和当前的方向AA正交,由于
    (25)ek+1=ek+αkdke_{k+1}=e_k+\alpha_kd_k \tag{25}
    所以
    (26)dkTA(ek+αkdk)=0&ThickSpace;&ThickSpace;αk=dkTAekdkTAdk=dkTrkdkTAdkek=xkx,Ax=bd_k^TA(e_k+\alpha_kd_k)=0 \implies \alpha_k=-\frac{d_k^TAe_k}{d_k^TAd_k}=\frac{d_k^Tr_k}{d_k^TAd_k}\\ e_k=x_k-x, Ax=b \tag{26}
    最优的步长因子和当前的梯度和方向向量均相关。从上边的表达式知道,当前反向和剩余误差方向是AA正交的,而我们也希望:
    (27)dkTAdk+1=0d_k^TAd_{k+1}=0 \tag{27}
    如果下一个方向向量和当前的也AA正交就好了,当然了,如果下一个方向向量dk+1d_{k+1}刚好选择到了剩余的误差向量ek+1e_{k+1},那么下一次迭代就会直接命中极值点了,但是一般这个是不可能的,因为ek+1e_{k+1}是依赖于αk\alpha_k的,会有循环求解的问题,除非我们已经知道解了,否则我们是得不到ek+1e_{k+1}的,而且在高维空间中,和dkd_k向量AA正交的向量是有很多个的,我们希望的是找到一个dk+1d_{k+1}满足上边表达式即可,于此同时,也满足
    (28)dkTAdm=0,m&lt;kd_k^TAd_m=0, \forall m \lt k \tag{28}
    这样每次迭代都减少了一个维度,那么,最多nn步就可以收敛了。
    这个过程,在变量为二维的情况下理解,大概如下图所示:

    这些椭圆是三维图形在二维的等势面,其实就是f(x)f(x)相等的值所构成的一个椭圆。
    一开始选择梯度作为初始方向,也就是等势椭圆的切线,下一步,由于取AA正交关系,d1d_1e1e_1是同向的,两步就命中结果了。对应到AA的空间中,这个椭圆就是个圆。变量是四维的情况,也可以用下图表示:

    变量为三维的情况,加上因变量实际上是四个维度的,由于二次型也可以描述能量关系,所以相等能量的各个三维自变量构成的图形在三维就是个椭球体。映射到AA的空间中其实就是个标准的球型。每次迭代搜寻的方向向量,都是和过去所有的方向向量AA正交的,三维的情况可以看得比较明显,那就是d1Td0=0,e1Td0=0,e1d1d_1^Td_0=0, e_1^Td_0=0, e_1 \ne d_1,也就是方向向量在高维度的情况下,即使和过去的方向向量都AA正交,但是也不能说这个方向向量就是当前的错误向量。但是从事实出发,错误向量ee随着每次迭代,都减少了各个方向向量的那部分,所以,错误向量e0e_0即,一开始,初始点到极值点的错误向量可以表示为:
    (29)e0=j=0n1δjdje_0=\sum_{j=0}^{n-1}\delta_jd_j \tag{29}
    而等式28中又限定了各个方向向量的关系,所以可以:
    (30)dkTAe0=δkdkTAdj=δkdkTAdk&ThickSpace;&ThickSpace;δk=dkTAe0dkTAdk=dkTA(e0+j=0k1αjdj)dkTAdk=dkTAekdkTAdk&ThickSpace;&ThickSpace;δk=dkTrkdkTAdk=αk d_k^TAe_0=\sum \delta_kd_k^TAd_j=\delta_kd_k^TAd_k \\ \implies \delta_k=\frac{d_k^TAe_0}{d_k^TAd_k}=\frac{d_k^TA(e_0 + \sum_{j=0}^{k-1}\alpha_jd_j)}{d_k^TAd_k}=\frac{d_k^TAe_k}{d_k^TAd_k} \\ \implies \delta_k=\frac{d_k^Tr_k}{d_k^TAd_k}=-\alpha_k \tag{30}
    所以
    (31)e0=j=0n1αjdj&ThickSpace;&ThickSpace;ek=e0+j=0k1αjdj=j=kn1αjdj e_0=-\sum_{j=0}^{n-1}\alpha_jd_j \\ \implies e_k=e_0 + \sum_{j=0}^{k-1}\alpha_jd_j=-\sum_{j=k}^{n-1}\alpha_jd_j \tag{31}
    上边的推到也证明,只要我们满足条件

    条件1:各个方向向量did_idjd_j之间两两AA正交
    条件2:通过求导让每一次迭代的步长因子α\alpha取得最优

    那么每一次迭代求解的过程,其实就是相当于减少了这个方向向量上的误差,每一步的误差向量可以被表示为方向向量的连加的形式,这样只需要迭代nn步后就会直接收敛到极值点。这个和我们在几何上边的理解是完全一致的。
    只是,我们该如何选择这些方向向量dkd_k呢?
    如果仅仅只是随便选择一组线性无关的构造基(u1,u2,...un)(u_1,u_2,...u_n),然后如同等式17那样构造方向向量,并满足AA正交,那么:
    (32)i&gt;j,diTAdj=uiTAdj+k=0i1βikdkTAdj=0&ThickSpace;&ThickSpace;βij=uiTAdjdjTAdj \forall i &gt; j, d_i^TAd_j=u_i^TAd_j+\sum_{k=0}^{i-1}\beta_{ik}d_k^TAd_j=0 \\ \implies \beta_{ij}=-\frac{u_i^TAd_j}{d_j^TAd_j} \tag{32}
    通过这一种方式,也可以得到如何去构造所有的方向向量,只是,任何方向向量did_i都和过去的i1i-1个向量相关,而且每次迭代还需要求解系数,这样就直接导致了计算所需的内存和时间开销巨大,没有实用价值。所以这一组构造基不应该随便选。而应该让等式32得出的各个系数拥有良好的关系,减少计算所需的开销。
    事实上,共轭方法,在一开始提出的时候,确实存在实现上的困难,不过通过仔细查看这些关系,我们还没有把梯度向量引入。
    在计算过程中,我们始终要得到当前点的梯度,在上边两个条件满足的情况下,观察错误向量,梯度等之间的关系:
    (33)i&lt;j,diTAej=j=pn1δjdiTAdj=diTrj&ThickSpace;&ThickSpace;diTrj=0 \forall i \lt j, -d_i^TAe_j=-\sum_{j=p}^{n-1}\delta_jd_i^TAd_j=-d_i^Tr_j \\ \implies d_i^Tr_j=0 \tag{33}
    也就是说,任何迭代后期的梯度rjr_j都是和之前的方向did_i在当前空间正交的,进一步的
    (34)diTrj=uiTrj+k=0i1βikdkTrj=0&ThickSpace;&ThickSpace;uiTrj=diTrj0,j=i&ThickSpace;&ThickSpace;uiTrj=diTrj=0,j&gt;i d_i^Tr_j=u_i^Tr_j+\sum_{k=0}^{i-1}\beta_{ik}d_k^Tr_j=0 \\ \implies u_i^Tr_j=d_i^Tr_j \ne 0, j = i \\ \implies u_i^Tr_j=d_i^Tr_j=0, \forall j \gt i \tag{34}
    结合前边的所有关系,可以得到如下图所示的关系图:

    图中所示,u2u_2d2d_2r2r_2的投影是相等的,所以末端所在的平面和d0d_0以及d1d_1所在的平面共面。通过上边的各个关系可以知道,没有必要取选择别的基uk{u_k},直接让uk=rku_k=r_k就是最优的选择,这个时候,各个向量的关系图将会变成如下所示:

    这样做的直接结果有:
    (35)αk=dkTrkdkTAdk=rkTrkdkTAdk \alpha_k=\frac{d_k^Tr_k}{d_k^TAd_k}=\frac{r_k^Tr_k}{d_k^TAd_k} \tag{35}
    (36)βij=uiTAdjdjTAdj=riTAdjdjTAdj,i&gt;jriTrj+1=riT(rjαjAdj)&ThickSpace;&ThickSpace;riTAdj=1αj(riTrjriTrj+1)&ThickSpace;&ThickSpace;riTAdj=0,βij=0,i&gt;j+1&ThickSpace;&ThickSpace;riTAdj=1αjriTri,βi,j+1=riTriri1Tri1,i=j+1 \beta_{ij}=-\frac{u_i^TAd_j}{d_j^TAd_j}=-\frac{r_i^TAd_j}{d_j^TAd_j},\forall i\gt j \\ r_i^Tr_{j+1}=r_i^T(r_j-\alpha_jAd_j) \\ \implies r_i^TAd_j=\frac{1}{\alpha_j}(r_i^Tr_j-r_i^Tr_{j+1}) \\ \implies r_i^TAd_j=0,\beta_{ij}=0,\forall i&gt;j+1 \\ \implies r_i^TAd_j=-\frac{1}{\alpha_j}r_i^Tr_i,\\ \beta_{i,j+1}=\frac{r_{i}^Tr_{i}}{r_{i-1}^Tr_{i-1}},i=j+1 \tag{36}
    综合前边所属,CG的流程可以被表示为:

    d0=r0=bAx0d_0=r_0=b-Ax_0
    for loop until reach max_step or cost_function below threshold:
    αk=rkTrkdkTAdk\qquad \alpha_k=\frac{r_k^Tr_k}{d_k^TAd_k}
    xk+1=xk+αkdk\qquad x_{k+1}=x_k+\alpha_kd_k
    rk+1=rkαkAdk or rk+1=bAxk+1\qquad r_{k+1}=r_k-\alpha_kAd_k\space or \space r_{k+1}=b-Ax_{k+1}
    βk+1=rk+1Trk+1rkTrk\qquad \beta_{k+1}=\frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k}
    dk+1=rk+1+βk+1dk\qquad d_{k+1}=r_{k+1}+\beta_{k+1}d_k
    return xnx_n

    使用Python代码实现:

    import numpy as np
    
    
    def conjugate_gradient(A, b, x_initial, max_step, threshold=0.00001):
        assert(isinstance(A, np.matrix))
        assert(isinstance(b, np.matrix))
        assert(isinstance(x_initial, np.matrix))
        r_old = b - A * x_initial
        d = r_old
        x = x_initial
        for i in range(max_step):
            alpha = (r_old.transpose() * r_old) / (d.transpose() * A * d)
            x = x + d * alpha
            # r_new = b - A * x
            r_new = r_old - A * d * alpha
            beta = (r_new.transpose() * r_new) / (r_old.transpose() * r_old)
            d = r_new + d * beta
            r_old = r_new
            cf = np.sqrt(np.sum(np.square(b - A * x)))
            if cf < threshold:
                print("Using step: ", i)
                break
        return x
    
    
    if __name__ == '__main__':
        N = 200
        Ar = np.mat(np.random.rand(N, N))
        As = Ar * Ar.transpose()
        bn = np.mat(np.random.rand(N, 1))
        xi = np.mat(np.random.rand(N, 1))
        xr = conjugate_gradient(As, bn, xi, 1000)
        print('1000:', np.sqrt(np.sum(np.square(bn - As * xr))))
        xr = conjugate_gradient(As, bn, xi, 10000)
        print('10000:', np.sqrt(np.sum(np.square(bn - As * xr))))
    
    

    运行上述程序的一个输出:

    Using step:  410
    1000: 6.43707958798e-06
    Using step:  410
    10000: 6.43707958798e-06
    

    大约需要两倍于维度的迭代次数才可以收敛到要求的范围,虽然理论上100维的变量值需要100次迭代,但是实际上各种浮点数运算误差导致不可能实现这样理想化的结果。通过对比会发现,迭代效果要比SD要好。当然了,实际使用的时候,判断损失函数其实没必要通过那样的形式。而且可以几十次迭代后才计算一下等等。

    基于SSOR预条件子的共轭梯度算法

    方程组的迭代求解

    上边的迭代算法,本质上其实还是在求解形如Ax=bAx=b的方程。毕竟最终的结果,都是让e=Axbe=Ax-b尽可能的接近于零。一般情况下,对于Ax=bAx=b的求解,通常直接通过求逆x=A1bx=A^{-1}b得到,这就存在求矩阵的逆可能不存在等等问题。如果AA不是方阵,且m&gt;&gt;nm \gt \gt n,那么就是在求解超定方程组,方程组可能没有解,求解Ax=bAx=b的过程就会等效于最小二乘,因为,为了让误差最小,其实就是让最终的的解Ax=bAx=b^{&#x27;}得到的bb^{&#x27;}bb最近,那么bb^{&#x27;}只能是bbAA的列空间的投影,所以它们之间的差和ATA^T的列空间正交,所以有
    (37)AT(Axb)=0&ThickSpace;&ThickSpace;x=(ATA)1ATb A^T(Ax-b)=0 \implies x = (A^TA)^{-1}A^Tb \tag{37}
    即使AA是非方阵,但是ATAA^TA却是对称方阵,所以求解逆是可行的,这个时候其实还是转换成了
    (38)ATAx=ATb,b=ATb,A=ATA&ThickSpace;&ThickSpace;Ax=b A^TAx=A^Tb,b^{&#x27;}=A^Tb, A^{&#x27;}=A^TA \implies A^{&#x27;}x=b^{&#x27;} \tag{38}
    所以还是可以使用上边的梯度迭代方法来求解。所以可以见得,线性系统问题,即可以表达为求解Ax=bAx=b的系统问题,还是比较多的。
    求解的方法,除了直接求逆,以及高斯的LU分解,其它的基本就是迭代方法了。最早的迭代有高斯-赛德尔迭代,SOR等,它们都是首先把AA拆分:
    (39)A=D+L+UA=D+L+U \tag{39}
    也就是把矩阵拆解成对角,上三角,下三角的形式,对于高斯-赛德尔方法,把(D+L+U)x=b(D+L+U)x=b变成了:
    (40)(L+D)xk+1=Uxk+bxk+1=D1(bUxkLxk+1) (L+D)x_{k+1}=-Ux_k+b \\ x_{k+1}=D^{-1}(b-Ux_k-Lx_{k+1})\tag{40}
    提取成对角,上下三角矩阵的原因,恐怕还是因为逆矩阵好算。上边的过程,每次迭代其实都是循环回代的过程,所以等式40右边有xk+1x_{k+1},其实是说要用当前解出来的那些新的变量值的意思。
    对于SOR来说,则引入了松弛系数ww,迭代过程是
    (41)xk+1=(wL+D)1[(1w)DxkwUxk]+w(D+wL)1b x_{k+1} =(wL+D)^{-1}[(1-w)Dx_k-wUx_k]+w(D+wL)^{-1}b \tag{41}
    上边的等式其实是用了w(D+L+U)x=wbw(D+L+U)x=wb转换而来。
    这些方法都有其收敛条件,在数值优化中有分析。
    对于稀疏矩阵,或者是对称正定矩阵来说,还是一般使用本篇总结的梯度迭代方法。只不过值得留意的却是,Ax=bAx=b竟然对应于一个二次型的极值求解过程

    Precondition

    在数值分析里边,矩阵迭代求解的收敛速度取决于矩阵条件数:
    (42)cond(A)=AA1cond(A)=||A|| \cdot ||A||^{-1} \tag{42}
    其中A||A||是矩阵范数,数值分析里边研究了矩阵迭代的收敛很大程度上取决于条件数,而预条件方法就是为了减少条件数,快速收敛。详细的分析很复杂,详情可看数值分析教材,这里只阐述用于共轭梯度算法的预条件过程。
    预条件是希望找到一个MM,将问题转化为求解:
    (43)M1Ax=M1bM^{-1}Ax=M^{-1}b \tag{43}
    其中MMn×nn \times n的可逆矩阵,称为预条件矩阵。在改进的共轭梯度方法中,这个矩阵MM是对称正定的。这个矩阵试图对AA逆转,从而让AA的条件数降低。此外,还引入了一种广义MM内积(v,w)M=vTMw(v,w)_M=v^TMw来取代欧几里得内积,所有内积默认都改成MM内积,原始的共轭梯度方法仍然成立,因为矩阵M1AM^{-1}A相对于新的内积仍然是对称正定矩阵:
    (44)(M1Av,w)M=vTAM1Mw=vTAw=vTMM1Aw=(v,M1Aw)M(M^{-1}Av,w)_M=v^TAM^{-1}Mw\\ =v^TAw=v^TMM^{-1}Aw=(v,M^{-1}Aw)_M \tag{44}
    那么,由于CG的推到基本都是基于vTAwv^TAw的形式,如果把等式44的内积结果应用到所有推导的过程,就有:
    (45)zk=M1bM1Axk=M1rkαk=(zk,zk)M(dk,M1Adk)Mxk+1=xk+αkdkzk+1=zkαkM1Adkβk=(zk+1,zk+1)M(zk,zk)Mdk+1=zk+1+βkdk z_k=M^{-1}b-M^{-1}Ax_k=M^{-1}r_k \\ \alpha_k=\frac{(z_k, z_k)_M}{(d_k, M^{-1}Ad_k)_M} \\ x_{k+1}=x_k+\alpha_kd_k \\ z_{k+1}=z_{k}-\alpha_kM^{-1}Ad_k \\ \beta_k=\frac{(z_{k+1},z_{k+1})_M}{(z_k,z_k)_M} \\ d_{k+1}=z_{k+1}+\beta_kd_k \tag{45}
    这里边一个有意思的地方在于,这个方法,把vTM1wv^TM^{-1}w的欧几里得的内积,替换成了(v,M1Aw)M(v, M^{-1}Aw)_M,定义内积在MM的逆空间里边,因为MM内积仍然满足普通内积该有的性质,所以可以完全替代欧几里得内积,公式仍然不发生变化。
    进一步转换可以得到:
    (46)(zk,zk)M=zkTMzk=zkTrk(dk,M1Adk)M=dkTMM1Adk=dkTAdk(zk+1,zk+1)M=zk+1Trk+1 (z_k,z_k)_M=z_k^TMz_k=z_k^Tr_k \\ (d_k,M^{-1}Ad_k)_M=d_k^TMM^{-1}Ad_k=d_k^TAd_k \\ (z_{k+1},z_{k+1})_M=z_{k+1}^Tr_{k+1} \tag{46}
    对于SSOR(对称连续过松弛)预条件子来说:
    (47)M=(D+wL)D1(D+wU),w(0,2)&ThickSpace;&ThickSpace;M1=(D+wU)1D(D+wL)1 M=(D+wL)D^{-1}(D+wU),w \in (0,2) \\ \implies M^{-1}=(D+wU)^{-1}D(D+wL)^{-1} \tag{47}
    总结上边的过程可以得到,SSOR预条件子的共轭梯度算法的流程是:

    r0=bax0,z0=M1r0,d0=z0r_0=b-ax_0,z_0=M^{-1}r_0,d_0=z_0
    for loop until reach max_step or cost_function below threshold:
    αk=zkTrkdkTAdk\qquad \alpha_k=\frac{z_k^Tr_k}{d_k^TAd_k}
    xk+1=xk+αkdk\qquad x_{k+1}=x_k+\alpha_kd_k
    rk+1=rkαkAdk\qquad r_{k+1}=r_k - \alpha_kAd_k
    zk+1=M1rk+1\qquad z_{k+1}=M^{-1}r_{k+1}
    βk+1=zk+1Trk+1zkTrk\qquad \beta_{k+1}=\frac{z_{k+1}^Tr_{k+1}}{z_k^Tr_k}
    dk+1=zk+1+βk+1dk\qquad d_{k+1}=z_{k+1}+\beta_{k+1}d_k
    return xnx_n

    对应的Python代码:

    import numpy as np
    
    def get_ssor_precondition_matrix(A, w):
        UD = np.triu(A)
        LD = np.tril(A)
        dim = A.shape[0]
        D = np.mat(np.zeros((dim, dim)))
        for i in range(dim):
            D[i, i] = A[i, i]
        for i in range(dim):
            for j in range(i+1, dim):
                UD[i, j] = w * UD[i, j]
        for i in range(dim):
            for j in range(0, i):
                LD[i, j] = w * LD[i, j]
        # 对上下三角矩阵求逆矩阵,其实不必用通用的求逆方法,不停回代即可
        return np.mat(np.linalg.inv(UD)) * D * np.mat(np.linalg.inv(LD))
    
    
    def precondition_conjugate_gradient(A, b, x_initial, max_step,
                                        threshold=0.00001, w=0.2):
        assert(isinstance(A, np.matrix))
        assert(isinstance(b, np.matrix))
        assert(isinstance(x_initial, np.matrix))
        r_old = b - A * x_initial
        M_inv = get_ssor_precondition_matrix(A, w)
        z_old = M_inv * r_old
        d = z_old
        x = x_initial
        for i in range(max_step):
            alpha = (z_old.transpose() * r_old) / (d.transpose() * A * d)
            x = x + d * alpha
            # r_new = b - A * x
            r_new = r_old - A * d * alpha
            z_new = M_inv * r_new
            beta = (z_new.transpose() * r_new) / (z_old.transpose() * r_old)
            d = z_new + d * beta
            r_old = r_new
            z_old = z_new
            cf = np.sqrt(np.sum(np.square(b - A * x)))
            if cf < threshold:
                print("Using step: ", i)
                break
        return x
    
    
    if __name__ == '__main__':
        N = 200
        Ar = np.mat(np.random.rand(N, N))
        As = Ar * Ar.transpose()
        bn = np.mat(np.random.rand(N, 1))
        xi = np.mat(np.random.rand(N, 1))
        xr = precondition_conjugate_gradient(As, bn, xi, 1000, 0.00001, 0.05)
        print('w=0.05:', np.sqrt(np.sum(np.square(bn - As * xr))))
        xr = precondition_conjugate_gradient(As, bn, xi, 1000, 0.00001, 0.5)
        print('w=0.5:', np.sqrt(np.sum(np.square(bn - As * xr))))
        xr = precondition_conjugate_gradient(As, bn, xi, 1000, 0.00001, 1)
        print('w=1:', np.sqrt(np.sum(np.square(bn - As * xr))))
    

    运行三次的结果:

    runfile('C:/Users/zczx1/.spyder-py3/temp.py', wdir='C:/Users/zczx1/.spyder-py3')
    Using step:  378
    w=0.05: 6.71578034331e-06
    Using step:  405
    w=0.5: 9.67223926338e-06
    Using step:  573
    w=1: 8.09438315554e-06
    
    runfile('C:/Users/zczx1/.spyder-py3/temp.py', wdir='C:/Users/zczx1/.spyder-py3')
    Using step:  371
    w=0.05: 6.1508910381e-06
    Using step:  401
    w=0.5: 8.51261715479e-06
    Using step:  602
    w=1: 7.57385104633e-06
    
    runfile('C:/Users/zczx1/.spyder-py3/temp.py', wdir='C:/Users/zczx1/.spyder-py3')
    Using step:  373
    w=0.05: 6.44798434081e-06
    Using step:  406
    w=0.5: 6.55294163482e-06
    Using step:  580
    w=1: 8.13783429322e-06
    

    虽然随机生成的矩阵不太能说明问题,但是这个只能说,在松弛系数比较小的时候,PCG和CG相比,相对来说减少了迭代的次数。但是当松弛系数较大的时候,还不如CG。

    非线性共轭梯度算法

    算法流程

    上边的总结,都是在二次型问题上进行的分析。可以看到,虽然算法本身的实现比较简单,但是算法背后的数学思路如果进行严谨地推导的话,还是有些复杂的,不过因为有了几何的理解,所以也还是比较直观。
    但是如果这个问题的表达式函数f(x)f(x)是非线性的,共轭梯度的思想仍然可以被推广,只是形式上不那么精确。而事实上,就算严格按照上边的算法流程,通过上边的实践也会发现,由于round off误差,其实不能达到理论的计算效果,所以在某种程度上近似求解也是可以的。回顾共轭梯度算法和核心,其实就两个:

    1. 依据步长和方向向量更新解向量,即xk+1=xk+αkdkx_{k+1}=x_k+\alpha_kd_k
    2. 利用梯度更新方向向量,即dk+1=rk+1+βk+1dkd_{k+1}=r_{k+1}+\beta_{k+1}d_{k}
      其中,rk=f(xk),d0=f(x0)r_k=-f^{&#x27;}(x_k),d_0=-f^{&#x27;}(x_0),我们只需要考虑怎么选择αk\alpha_kβk+1\beta_{k+1}即可
      αk\alpha_k,是希望找到这样的步长因子使得f(xk+αkdk)Tdk=0f^{&#x27;}(x_k+\alpha_kd_k)^Td_k=0,本质上还是希望下一个梯度方向和当前的方向向量正交,这个时候应该就是最优的解,和SD中的线搜索是类似的。不过,在实际的实现中,可能只是设置一个固定值,然后乘以一个比例,这个比例随着迭代的进行可能会变得越来越小。所以:
      (48)αk=argmin(f(xk+αkdk) or αk=C×radio \alpha_k=argmin(f(x_k + \alpha_kd_k) \space or \space \alpha_k=C\times radio \tag{48}
      而对于βk+1\beta_{k+1},不同的学者提出了使用不同的公式,比如:
      (49)FletcherReeves:βk+1FR=ri+1Tri+1riTri Fletcher-Reeves: \beta_{k+1}^{FR}=\frac{r_{i+1}^Tr_{i+1}}{r_i^Tr_i} \tag{49}
      (50)PolakRibiere:βk+1PR=max(0,rk+1T(rk+1rk)rkTrk) Polak-Ribiere: \beta_{k+1}^{PR}=max(0, \frac{r^T_{k+1}(r_{k+1}-r_{k})}{r_k^Tr_k}) \tag{50}
      也还有很多其它的计算方法,不过显然,这些只是一种近似求解的方法。因为非线性的损失函数各种各样,精确的求解不现实。这样,非线性CG的流程如下:

    d0=r0=f(x0)d_0=r_0=-f^{&#x27;}(x_0)
    for loop until max step or f(x)&lt;threshf^{&#x27;}(x) \lt thresh:
    αk=argmin(f(xk+αkdk))\qquad \alpha_k=argmin(f(x_k+\alpha_kd_k))
    xk+1=xk+αkdk\qquad x_{k+1}=x_k+\alpha_kd_k
    rk+1=f(xk+1)\qquad r_{k+1}=-f^{&#x27;}(x_{k+1})
    βk+1PR=max(0,rk+1T(rk+1rk)rkTrk)\qquad \beta_{k+1}^{PR}=max(0, \frac{r^T_{k+1}(r_{k+1}-r_{k})}{r_k^Tr_k})
    dk+1=rk+1+βk+1dk\qquad d_{k+1}=r_{k+1}+\beta_{k+1}d_k

    非线性函数中最优α\alpha的选择

    寻求最优的步长因子的过程其实就是广义线搜索,根据泰勒级数展开,有:
    (51)f(x+αd)f(x)+α[ddαf(x+αd)]α=0+α22[d2dα2f(x+αd)]α=0=f(x)+α[f(x)]Td+α22dTf(x)d f(x+\alpha d) \approx f(x)+\alpha [\frac{d}{d\alpha}f(x+\alpha d)]_{\alpha=0}+\frac{\alpha^2}{2}[\frac{d^2}{d\alpha^2}f(x+\alpha d)]_{\alpha=0} \\ = f(x)+\alpha[f^{&#x27;}(x)]^Td+\frac{\alpha^2}{2}d^Tf^{&#x27;&#x27;}(x)d \tag{51}
    (52)ddαf(x+αd)[f(x)]Td+αdTf(x)d=0&ThickSpace;&ThickSpace;α=f(x)TddTf(x)d \frac{d}{d\alpha}f(x+\alpha d) \approx [f^{&#x27;}(x)]^Td+\alpha d^Tf^{&#x27;&#x27;}(x)d=0 \\ \implies \alpha=-\frac{f^{&#x27;}(x)^Td}{d^Tf^{&#x27;&#x27;}(x)d} \tag{52}
    从等式52中,可以看到二次型中求解最优α\alpha的影子,如果是二次型,那么就是α=rTddTAd\alpha=\frac{r^Td}{d^TAd}的形式了,不过对于非线性问题来说,就需要计算二阶导数。二阶导数为:
    (52)f(x)=[2fx1x12fx1x22fx1xn2fx2x12fx2x22fx2xn2fxnx12fxnx22fxnxn] f^{&#x27;&#x27;}(x)= \begin{bmatrix} \\ \frac{\partial^2f}{\partial x_1 \partial x_1} &amp; \frac{\partial^2f}{\partial x_1 \partial x_2} &amp; \cdots &amp; \frac{\partial^2f}{\partial x_1 \partial x_n}\\ \frac{\partial^2f}{\partial x_2 \partial x_1} &amp; \frac{\partial^2f}{\partial x_2 \partial x_2} &amp; \cdots &amp; \frac{\partial^2f}{\partial x_2 \partial x_n}\\ \vdots &amp; &amp; \ddots &amp; \\ \frac{\partial^2f}{\partial x_n \partial x_1} &amp; \frac{\partial^2f}{\partial x_n \partial x_2} &amp; \cdots &amp; \frac{\partial^2f}{\partial x_n \partial x_n}\\ \end{bmatrix} \tag{52}
    其实就是求解Hessian矩阵。这个计算量可谓真的大,而且随着迭代,这个矩阵的具体数值,可能还是需要不断计算的,这个开销几乎是不可忍受的,所以严格的计算步长因子基本是不可取的。何况有些时候,可能都求不出二阶导,也无怪乎实际应用中,直接给个固定步长乘以比例了事。与其算这个矩阵那么复杂,不如多迭代几次。

    非线性函数求导简介

    非线性CG需要对函数进行求导,而求导的方法,目前有:

    1. 公式法
    2. 符号微分法
    3. 自动微分法
      这些根据具体场景来实现,只是非线性CG迭代流程就是如上所示的过程。
      值得注意的是,非线性函数的这个过程,最终收敛到的可能只是局部极值,和初始值的选择相关。如果初始值在一个类二次型的区域,那么就可以收敛到这块区域的极值点。

    总结

    本文总结了本人学习EDA软件中布局算法中使用到的共轭梯度相关的内容。共轭梯度算法,应用到实际工程当中,可能还是使用非线性的部分内容,只是由于也是近似的求解,所以和随机梯度下降相比也没有太多的优势,毕竟都是和初始点相关,所以引入更多的随机性,也许可以取得更好的效果!

    参考资料

    1. An Introduction to the Conjugate Gradient Method Without the Agonizing Pain
    2. Numerical Analysis

    @fsfzp888
    2019 年 01月

    展开全文
  • 共轭梯度法:共轭梯度
  • 共轭梯度法.py
  • 共轭梯度
  • 共轭梯度共轭梯度法原理共轭梯度法原理共轭梯度法原理共轭梯度法原理
  • matlab共轭梯度法求目标函数的最小极值-共轭梯度-王.rar 我是地球物理专业的一名学生,把自己实习的作业发上来大家分享下
  • 共轭梯度的BP算法.rar共轭梯度的BP算法.rar共轭梯度的BP算法.rar共轭梯度的BP算法.rar共轭梯度的BP算法.rar共轭梯度的BP算法.rar共轭梯度的BP算法.rar
  • GCG广义共轭梯度法X = GCG(A,B) 尝试求解线性方程组 A*X=B 对于 X。 N×N 系数矩阵 A 必须对称且为正确定且右侧列向量 B 的长度必须为 N。 X = GCG(A,B,TOL) 指定方法的容差。 如果 TOL 是 [] 然后 GCG 使用默认值 ...
  • MATLAB实现共轭梯度

    2018-04-17 21:17:59
    MATLAB实现共轭梯度法MATLAB实现共轭梯度法MATLAB实现共轭梯度
  • 共轭梯度算法Matlab

    2020-05-03 22:04:06
    基于共轭梯度法的详细案例,共轭梯度法是最优化方法的其中一种优化方案。通过变分法求解线性方程组。方向是在求出梯度方向的前提下,添加正则项,使得前后两次方向互为共轭所得出的方向向量。
  • 共轭梯度求解

    2018-10-24 16:53:36
    共轭梯度方法求解优化问题。目标函数是高维二次函数。
  • 共轭梯度法MATLAB程序%conjugate gradient methods%method:FR,PRP,HS,DY,CD,WYL,LS%精确线搜索,梯度终止准则function [ m,k,d,a,X,g1,fv] conjgradme G,b,c,X,e,methodif nargin 6 error '输入参数必须为6' ;...

    共轭梯度法MATLAB程序

    %conjugate gradient methods

    %method:FR,PRP,HS,DY,CD,WYL,LS

    %精确线搜索,梯度终止准则

    function [ m,k,d,a,X,g1,fv] conjgradme G,b,c,X,e,method

    if nargin 6 error '输入参数必须为6' ;

    end

    n length G ;

    if n 2

    format long e %rat

    syms x1 x2

    f 1/2*[x1,x2]*G*[x1;x2]+b'*[x1;x2]+c;

    g [diff f,x1 ;diff f,x2 ];

    g1 subs subs g,x1,X 1,1 ,x2,X 2,1 ;

    d -g1;

    a - d'*g1 / d'*G*d ;% a - X :,1 '*G*d+b'*d / d'*G*d ; a g1 :,1 '*g1 :,1 / d :,1 '*G*d :,1 ;

    X :,2 X :,1 +a*d;

    g1 [g1 subs subs g,x1,X 1,2 ,x2,X 2,2 ];

    m1 norm g1 :,1 ;

    m norm g1 :,2 ;

    i 2;

    k zeros 1 ;

    switch method case 'FR' while m e k i-1 m/m1 ^2; d :,i -g1 :,i +k i-1 *d :,i-1 ; a i - d :,i '*g1 :,i / d :,i '*G*d :,i ; %a1 i - X :,i '*G*d :,i +b'*d :,i / d :,i '*G*d :,i ;a i g1 :,i '*g1 :,i / d :,i '*G*d :,i ; X :,i+1 X :,i +a i *d :,i ; g1 [g1 subs subs g,x1,X 1,i+1 ,x2,X 2,i+1 ]; m1 m; m norm g1 :,i+1 ; i i+1; end case 'PRP' while m e k i-1 g1 :,i '* g1 :,i -g1 :,i-1 / norm g1 :,i-1 ^2; d :,i -g1 :,i +k i-1 *d :,i-1 ; a i - d :,i '*g1 :,i / d :,i '*G*d :,i ; X :,i+1 X :,i +a i *d :,i ; g1 [g1 subs subs g,x1,X 1,i+1 ,x2,X 2,i+1 ]; m norm g1 :,i+1 ; i i+1; end case 'HS' while m e k i-1 g1 :,i '* g1 :,i -g1 :,i-1 / d :,i-1 '* g1 :,i -g1 :,i-1 ; d :,i -g1 :,i +k i-1 *d :,i-1 ; a i - d :,i '*g1 :,i / d :,i '*G*d :,i ; X :,i+1 X :,i +a i *d :,i ; g1 [g1 subs subs g,x1,X 1,i+1 ,x2,X 2,i+1 ]; m norm g1 :,i+1 ; i i+1; end case 'DY' while m e k i-1 g1 :,i '*g1 :,i / d :,i-1 '* g1 :,i -g1 :,i-1 ; d :,i -g1 :,i +k i-1 *d :,i-1 ; a i - d :,i '*g1 :,i / d :,i '*G*d :,i ; X :,i+1 X :,i +a i *d :,i ; g1 [g1 subs subs g,x1,X 1,i+1 ,x2,X 2,i+1 ]; m norm g1 :,i+1 ; i i+1; end case 'LS' while m e k i-1 g1 :,i '* g1 :,i -g1 :,i-1 / d :,i-1 '* -g1 :,i-1 ; d :,i -g1 :,i +k i-1 *d :,i-1 ; a i - d :,i '*g1 :,i / d :,i '*G*d :,i ; %a i - X :,i '*G*d :,i +b'*d :,i / d :,i '*G*d :,i ; X :,i+1 X :,i +a i *d :,i ; g1 [g1 subs subs g,x1,X 1,i+1 ,x2,X 2,i+1 ]; m norm g1 :,i+1 ; i i+1; end

    展开全文
  • 利用共轭梯度和投影梯度对函数的优化
  • python实现共轭梯度

    2021-01-20 06:05:02
    共轭梯度法是介于最速下降法与牛顿法之间的一个方法,它仅需利用一阶导数信息,但克服了最速下降法收敛慢的缺点,又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点,共轭梯度法不仅是解决大型线性方程组最有用的...
  • matlab共轭梯度

    2020-12-24 20:19:39
    适合初学用matlab编写共轭梯度
  • 共轭梯度

    2014-05-24 12:51:09
    了解学习共轭梯度法,英文
  • 对称佩里共轭梯度

    2021-03-18 04:02:21
    对称佩里共轭梯度
  • 共轭梯度法 共轭方向法中,最大的问题在于寻找一组彼此独立的向量u1,...,unu1,...,unu_1,...,u_n,如果选取不当,那么和高斯消元法就没有区别了。共轭梯度法 (Conjugate Gradient)实际上是一种特殊的共轭方向法,它...

    共轭梯度法

    共轭方向法中,最大的问题在于寻找一组彼此独立的向量u1,...,un,如果选取不当,那么和高斯消元法就没有区别了。共轭梯度法 (Conjugate Gradient)实际上是一种特殊的共轭方向法,它取ui=r(i)

    首先为何残差是彼此独立的?首先,由于共轭方向法每次(A正交地)消除了一个维度上的误差,所以有

    diTAe(j)=0(i<j)

    由于r(j)=Ae(j),于是有diTr(j)=0。另外,由于现在的ui=r(i),所以

    span{d1,...,dn}=span{ui,...,un}=span{r(1),...,r(n)}

    r(j)di(i<j)张成的子空间正交,因此也和所有的r(i)(i<j)正交。

    共轭梯度法和最速下降法的区别

    这里我产生的一个疑问是,残差其实就是梯度的(反)方向,既然当次迭代的梯度方向总是和之前的梯度方向垂直,那么它和最速下降法有什么区别?为何最速下降法会出现“之”字形的迭代路线?我认为是这样的。首先,在最速下降法中,当次迭代的梯度方向也是和上次迭代梯度方向垂直,但和再之前的梯度方向就不垂直了,所以会有“之”形路线。然后,共轭梯度法要求的是关于矩阵正交,并非直接正交。在上一节的图中,我们可以看到关于矩阵正交的形状,通过拉伸,可以看到,它才能真正的分离维度,普通的正交反而是有偏差的。共轭方向法将本来普通正交的梯度向量,调整为新的关于矩阵正交,虽然张成的空间是一样的,但方向更加合理,所以能保证在n步内迭代收敛。

    共轭梯度法的优势

    于是我们可以看到,起码共轭梯度法作为一种特殊的共轭方向法,是没有错误的。那共轭梯度法有什么优势呢?

    先来看一个有趣的性质。在上一节的末尾,我们提到了迭代更新残差的方法:

    r(i+1)=r(i)α(i)Adi

    注意到新的残差(或梯度)引入了一个新的维度(因为新的搜索方向和之前所有梯度张成的空间关于矩阵正交)。即假设span{r(0),...,r(i+1)}=D(i+1),那么有rank(D(i+1))=rank(D(i))+1。显然r(i)D(i),新维度只能来源于Adi,而diD(i),所以A带来的变换引入了新维度。引申开来,不难得到

    D(i)=span{d1,...,dn}=span{d1,Ad1,A2d1,...,Ai1d1}=span{r(1),Ar(1),A2r(1),...,Ai1r(1)}

    诸如此类的,由一个向量反复乘以同样的矩阵变换张成的子空间被成为克雷洛夫子空间 (Krylov Subspace)。这样的子空间有一个很好的性质,那就是施密特正交化法会变得非常容易。

    在上一节中,我们给出了施密特正交化法的参数表达式,

    βij=r(i)AdjdjTAdj

    由于

    r(i)Tr(j+1)=r(i)T(r(j)α(j)Adj)=r(i)Tr(j)α(j)r(i)TAdjr(i)TAdj=1α(j)[r(i)Tr(j)r(i)Tr(j+1)]

    由于i>j,所以

    (i=j+1):βij=1α(j)r(i)Tr(i)djTAdj(i>j+1):βij=0

    这样,共轭梯度法作为一种特殊的共轭方向法,解决了施密特正交化法速度慢的问题。当然α(j)的表达式和之前是一致的,即

    α(j)=djTr(j)djTAdj

    代入