精华内容
下载资源
问答
  • 应用于正定对称稀疏系数矩阵的方程求解,对特定类型方程求解进行加速
  • 共轭梯度法简介

    2019-11-12 18:39:21
    共轭梯度法简介 关键词:共轭梯度法;预处理;收敛分析 文章目录共轭梯度法简介introductionnotationThe Quadratic Form(二次型)The Method of Steepest Descent introduction Conjugate Gradient Method...

    共轭梯度法简介

    关键词:共轭梯度法;预处理;收敛分析

    1.introduction

    Conjugate Gradient Method(henceforth,CG),CG这样的迭代方法适合用于稀疏矩阵

    共轭梯度法是求解稀疏线性方程组最突出的迭代方法。CG是求解大型线性方程组最常用的迭代方法。CG是有效的系统的形式 A x = b Ax=b Ax=b ; A A A: 对称,正定的方形矩阵。

    引入二次型的概念,用于推导最速下降法、共轭方向法和共轭梯度法。对特征向量进行了解释,并用它们来检验雅可比矩阵法、陡降法和共轭梯度法的收敛性。其他主题包括预处理和非线性共轭梯度法。

    文章结构
    文章结构

    2.notation

    用大写字母表示矩阵,用小写字母表示向量,用希腊字母表示标量。 A A A N × N N\times N N×N的矩阵, x x x b b b n × 1 n \times 1 n×1的矩阵。满矩阵方程1为:

    两个向量的内积可写成 x T y x^Ty xTy表示标量和: ∑ i = 1 n x i y i \sum^{ n}_{i=1}{x_iy_i} i=1nxiyi,其中 x T y = y T x x^Ty=y^Tx xTy=yTx.

    采样二维线性系统。解位于这两条直线的交点上。图一
    在这里插入图片描述
    对于每个非零向量 x x x, x T A x > 0 x^TAx>0 xTAx>0,

    3.The Quadratic Form(二次型)

    二次型只是一个标量,一个向量的二次函数 f ( x ) = 1 2 x T A x − b T x + c f(x)=\frac{1}{2}x^TAx-b^Tx +c f(x)=21xTAxbTx+c这里 A A A是一个矩阵, x x x b b b是向量, c c c是标量常数。如果A是对称正定的,则 f ( x ) f(x) f(x) A x = b Ax=b Ax=b的极小化解。

    我将用简单的示例问题来演示思想:
    在这里插入图片描述
    这是求上面两条直线交点的 A x = b Ax=b Ax=b。一般,解 x x x位于 n n n维超平面的交点上,每个超平面都有 n − 1 n-1 n1维。对于此问题,解是 x = [ 2 , − 2 ] T x=[2,-2]^T x=[2,2]T.对应的二次型 f ( x ) f(x) f(x)如下图所示:这个曲面的最小点是 A x = b Ax=b Ax=b
    图二:
    在这里插入图片描述
    下面是二次型 f ( x ) f(x) f(x)等高线图,每个椭圆曲线都有一个常数图三
    在这里插入图片描述
    梯度 f ′ ( x ) f'(x) f(x),对于每个 x x x,梯度点在 f ( x ) f(x) f(x)最陡上升方向,与等高线正交.如下图:图四

    在这里插入图片描述
    因为 A A A是正定的, f ( x ) f(x) f(x)的表面形状是形状像抛物面碗的。

    二次型的梯度定义为
    在这里插入图片描述
    f ′ ( x ) f'(x) f(x)等于零时, f ( x ) f(x) f(x)可以极小化。将上式带入
    f ( x ) = 1 2 x T A x − b T x + c f(x)=\frac{1}{2}x^TAx - b^Tx+c f(x)=21xTAxbTx+c
    过程
    得到: f ′ ( x ) = 1 2 A T x + 1 2 A x − b f'(x)=\frac{1}{2}A^Tx + \frac{1}{2}Ax-b f(x)=21ATx+21Axb
    如果 A A A是对称矩阵,方程化为:
    f ′ ( x ) = A x − b f'(x)=Ax-b f(x)=Axb
    将梯度化为零,就得到我们要解的线性方程 A X = b . AX=b. AX=b.因此 A x = b Ax=b Ax=b的解是 f ( x ) f(x) f(x)的一个重要点。如果 A A A正定且对称,解就是 f ( x ) f(x) f(x)的最小值,也就是说 A x = b Ax=b Ax=b可以通过找 f ( x ) f(x) f(x)的最小值来解出.(如果 A A A不是对称的,方程 f ′ ( x ) = 1 2 A T x + 1 2 A x − b f'(x)=\frac{1}{2}A^Tx + \frac{1}{2}Ax-b f(x)=21ATx+21Axb表明CG将要找到关于系统 1 2 ( A T + A ) x = b \frac{1}{2}(A^T+A)x=b 21(AT+A)x=b的解.可以知道 1 2 ( A T + A ) \frac{1}{2}(A^T+A) 21(AT+A)是对称矩阵。)

    为什么对称正定矩阵有比较好的特征?考虑 f f f在任意点 p p p和在解点 x = A − 1 b x=A^{-1}b x=A1b之间的关系.如果 A A A是对称的(正定与否),方程 f ( x ) = 1 2 x T A x − b T x + c f(x)=\frac{1}{2}x^TAx-b^Tx +c f(x)=21xTAxbTx+c可以看出
    在这里插入图片描述
    f ( p ) = f ( x ) + 1 2 ( p − x ) T A ( p − x ) f(p)=f(x)+\frac{1}{2}(p-x)^TA(p-x) f(p)=f(x)+21(px)TA(px)
    如果 A A A是正定的,通过不等式 x T A x > 0 x^TAx>0 xTAx>0,对于所有 p ≠ x p\neq x p=x后面一项都是正的。于是出现了 x x x f f f的一个全局最小化。

    f ( x ) f(x) f(x)是一个抛物面这是我们对于什么是正定矩阵最好的解释。如果 A A A不是正定的,还有其他几种可能性:
    1、 A A A可能是负定: 负的正定矩阵。
    2、 A A A可能是奇异的: 这种情况下解不唯一;对于 f f f解集是具有一致值的直线或超平面。
    如果 A A A以上几种情况都不是, x x x是一个鞍点,像最速下降和CG这样的技术可能会失败。下图展示了可能性。 b b b c c c的值确定抛物面的最小点在哪里,但不影响抛物面的形状。
    图五在这里插入图片描述
    (a)正定矩阵的二次型。(b)负定矩阵。©一个奇异(和正不定)矩阵。穿过山谷底部的一条线是一组解。(d)不定矩阵。因为解是鞍点,所以最陡下降和重心不成立。在三维或更高的空间中,奇异矩阵也可以有鞍。

    为什么要麻烦地把线性系统转换成一个看起来更复杂的问题呢?最小化问题比如图二,而不是关于相交的超平面比如图一。

    4.The Method of Steepest Descent

    在最陡下降法中,我们从任意一点 x ( 0 ) x_{(0)} x(0)开始,然后滑到抛物面的底部.我们采取一系列迭代步骤 x ( 1 ) , x ( 2 ) , . . . x_{(1)},x_{(2)},... x(1),x(2),...直到我们已经足够接近解 x x x.

    第一步,我们选择 x x x下降最快的方向,就是 f ′ ( x ) f'(x) f(x)的反方向。根据方程 f ′ ( x ) = A x − b f'(x)=Ax-b f(x)=Axb,则这里的方向为 − f ′ ( x ( i ) ) = b − A x ( i ) -f'(x_{(i)})=b-Ax_{(i)} f(x(i))=bAx(i)

    定义 :
    1.误差error : e ( i ) = x ( i ) − x e_{(i)}=x(i)-x e(i)=x(i)x表示距离解有多远
    2.残差residual: r ( i ) = b − A x ( i ) r_{(i)}=b-Ax_{(i)} r(i)=bAx(i) 表示离正确值 b b b的距离。很容易观察到 r ( i ) = − A e ( i ) r_{(i)}=-Ae_{(i)} r(i)=Ae(i),你应该把残差看成是由 A A A转换成与 b b b相同空间的误差.
    3. r ( i ) = − f ′ ( x ( i ) ) r_{(i)}=-f'(x_{(i)}) r(i)=f(x(i)).
    应该把残差看作是最陡下降的方向。对于非线性问题,在第14节,只有后者的定义适用。所以记住,当你读到“残余”时,想想“最陡下降的方向”。

    假设从 x ( 0 ) = [ − 2 , − 2 ] T x_{(0)}=[-2,-2]^T x(0)=[2,2]T开始迭代。第一步:沿着最陡的下降方向,将落在图6(a)中的实线上。换句话说,我们选择一个点 x ( 1 ) = x ( 0 ) + α r ( 0 ) x_{(1)}=x_{(0)}+\alpha r_{(0)} x(1)=x(0)+αr(0)
    应该取多大一步?

    行搜索是选择 α \alpha α 使 f f f沿直线最小化的过程。图6(b)说明了这个任务:我们被限制在垂直面和抛物面的交点上选择一个点。图6©是由这些曲面的交点定义的抛物线。 α \alpha α在抛物线底部的值是多少?

    从基本的微积分,当方向导数 d d α f ( x ( 1 ) ) = 0 \frac{d}{d\alpha}f(x_{(1)})=0 dαdf(x(1))=0时, 通过链式法则,
    d d α f ( x ( 1 ) ) = f ′ ( x ( 1 ) ) T d d α x ( 1 ) = f ′ ( x ( 1 ) ) T r ( 0 ) \frac{d}{d\alpha}f(x_{(1)})=f'(x_{(1)})^T\frac{d}{d\alpha}x_{(1)}=f'(x_{(1)})^Tr_{(0)} dαdf(x(1))=f(x(1))Tdαdx(1)=f(x(1))Tr(0)
    为了上式为零, α \alpha α应该选择 r ( 0 ) r_{(0)} r(0) f ′ ( x ( 1 ) ) f'(x_{(1)}) f(x(1))正交的。(图6(d))

    有一个直观的原因为什么我们期望这些向量在最小值处是正交的。图7显示了搜索线上各点的梯度向量。抛物线的斜率
    (图6©))在任意点上等于梯度投影到直线上的大小(图7,虚线箭头)。这些投影表示当一个人穿过搜索线时 f f f的增长率. f f f在投影为0的地方最小化——在梯度垂直于搜索线的地方。

    为确定 α \alpha α,通过 f ′ ( x ( 1 ) ) = − r ( 1 ) f'(x_{(1)})=-r_{(1)} f(x(1))=r(1),得到
    r ( 1 ) T r ( 0 ) = 0 r_{(1)}^Tr_{(0)} = 0 r(1)Tr(0)=0 ( b − A x ( 1 ) ) T r ( 0 ) = 0 (b-Ax_{(1)})^Tr_{(0)}=0 (bAx(1))Tr(0)=0 ( b − A ( x ( 0 ) + α r ( 0 ) ) ) T r ( 0 ) = 0 (b-A(x_{(0)}+\alpha r_{(0)}))^Tr_{(0)}=0 (bA(x(0)+αr(0)))Tr(0)=0 ( b − A x ( 0 ) ) T r ( 0 ) − α ( A r ( 0 ) ) T r ( 0 ) = 0 (b-Ax_{(0)})^Tr_{(0)}-\alpha (Ar_{(0)})^Tr_{(0)}=0 (bAx(0))Tr(0)α(Ar(0))Tr(0)=0 ( b − A x ( 0 ) ) T r ( 0 ) = α ( A r ( 0 ) ) T r ( 0 ) (b-Ax_{(0)})^Tr_{(0)}=\alpha (Ar_{(0)})^Tr_{(0)} (bAx(0))Tr(0)=α(Ar(0))Tr(0) r ( 0 ) T r ( 0 ) = α r ( 0 ) T ( A r ( 0 ) ) r_{(0)}^Tr_{(0)}=\alpha r_{(0)}^T(Ar_{(0)}) r(0)Tr(0)=αr(0)T(Ar(0)) α = r ( 0 ) T r ( 0 ) r ( 0 ) T A r ( 0 ) \alpha=\frac{r_{(0)}^Tr_{(0)}}{r_{(0)}^TAr_{(0)}} α=r(0)TAr(0)r(0)Tr(0)

    在这里插入图片描述
    图6: 最陡下降法。(a)从 [ − 2 , 2 ] T [-2,2]^T [2,2]T开始,在 f f f的最陡下降方向迭代一步。(b)求出这两个曲面的交点处f最小的点。( c)这个抛物线是曲面的交点,最下面点是我们所要的。(d)最下面一点的梯度与前一步的梯度是正交的。
    在这里插入图片描述
    图7: 梯度 f ′ f' f显示在搜索线上的几个位置(实箭头)。每个梯度在直线上的投影也显示出来(虚线箭头)。梯度向量表示 f f f最陡增长的方向,这些投影表示当一个人穿过搜索线时的增长率。 f f f是当梯度垂直于搜索线时最小化。
    在这里插入图片描述
    图8: [ − 2 , 2 ] T [-2,2]^T [2,2]T开始使用最陡下降法并收敛在 [ 2 , − 2 ] T [2,-2]^T [2,2]T.

    综上所述,最陡下降法是:
    r ( i ) = b − A x ( i ) , r_{(i)}=b-Ax_{(i)}, r(i)=bAx(i), α ( i ) = r ( i ) T r ( i ) r ( i ) T A r ( i ) \alpha_{(i)}=\frac{r_{(i)}^Tr_{(i)}}{r_{(i)}^TAr_{(i)}} α(i)=r(i)TAr(i)r(i)Tr(i) x ( i + 1 ) = x ( i ) + α ( i ) r ( i ) x_{(i+1)}=x_{(i)}+\alpha_{(i)}r_{(i)} x(i+1)=x(i)+α(i)r(i)
    运行该示例,直到它在图8中收敛。注意之字形路径,它的出现是因为每个梯度与前一个梯度是正交的。

    如上所述,该算法每次迭代需要两次矩阵-向量乘法。最速下降的计算成本主要由矩阵向量乘积决定;幸运的是,有一个是可以消除的。
    将方程 x ( i + 1 ) = x ( i ) + α ( i ) r ( i ) x_{(i+1)}=x_{(i)}+\alpha_{(i)}r_{(i)} x(i+1)=x(i)+α(i)r(i)两边左乘 A A A再加 b , b, b,得到: r ( i + 1 ) = r ( i ) − α ( i ) A r ( i ) . r_{(i+1)}=r_{(i)}-\alpha_{(i)}Ar_{(i)}. r(i+1)=r(i)α(i)Ar(i).
    虽然方程 r ( i ) = b − A x ( i ) , r_{(i)}=b-Ax_{(i)}, r(i)=bAx(i),还需要计算 r ( 0 ) , r_{(0)}, r(0),方程 r ( i + 1 ) = r ( i ) − α ( i ) A r ( i ) . r_{(i+1)}=r_{(i)}-\alpha_{(i)}Ar_{(i)}. r(i+1)=r(i)α(i)Ar(i).可用于此后的每次迭代。 A r Ar Ar的乘积出现在方程 α ( i ) = r ( i ) T r ( i ) r ( i ) T A r ( i ) \alpha_{(i)}=\frac{r_{(i)}^Tr_{(i)}}{r_{(i)}^TAr_{(i)}} α(i)=r(i)TAr(i)r(i)Tr(i) r ( i + 1 ) = r ( i ) − α ( i ) A r ( i ) r_{(i+1)}=r_{(i)}-\alpha_{(i)}Ar_{(i)} r(i+1)=r(i)α(i)Ar(i),只需要被计算一次。使用这种递归的缺点是,由方程 r ( i + 1 ) = r ( i ) − α ( i ) A r ( i ) . r_{(i+1)}=r_{(i)}-\alpha_{(i)}Ar_{(i)}. r(i+1)=r(i)α(i)Ar(i).定义的序列生成时没有任何来自值 x ( i ) x_{(i)} x(i)的反馈,因此,浮点舍入误差的累积可能导致 x ( i ) x_{(i)} x(i)收敛到附近的某个点 x x x.这种影响可以通过周期性地使用方程 r ( i ) = b − A x ( i ) r_{(i)}=b-Ax_{(i)} r(i)=bAx(i)重新计算正确的余差来避免。

    在分析最陡下降的收敛性之前,我必须离题,以确保你对特征向量有一个坚实的理解。

    5.Thinking with Eigenvectors and Eigenvalues

    特征向量主要用作分析工具;最速下降和CG不计算任何特征向量的值作为算法的一部分.

    5.1 Eigen do it if I try

    当矩阵 B B B作用于他的非零向量 v v v时,特征向量 v v v不旋转被应用于它(除了可能指向恰好相反的方向)。特征向量 v v v可能会改变长度或改变方向,但不会转向侧面。也就是说,有一个标量常数 λ \lambda λ使得 B v = λ v . Bv=\lambda v. Bv=λv. λ \lambda λ也就是 B B B的特征值。对于任意常数 α \alpha α,向量 α v \alpha v αv也是特征值 λ \lambda λ的一个特征向量.因为 B ( α v ) = α B v = λ α v . B(\alpha v)=\alpha Bv=\lambda \alpha v. B(αv)=αBv=λαv.也就是说,如果你缩放一个特征向量,它仍然是一个特征向量。

    迭代方法通常依赖于一次又一次地对向量应用 B B B。当 B B B重复地应用于特征向量 v v v时,两种情况之一会发生:
    1、如果 ∣ λ ∣ < 1 |\lambda|<1 λ<1,当 i i i趋近于无穷, B i v = λ i v B^iv=\lambda^iv Biv=λiv将会消失(图九)。
    2、如果 ∣ λ ∣ > 1 |\lambda|>1 λ>1,当 i i i趋于无穷时, B i v = λ i v B^iv=\lambda^iv Biv=λiv将会趋近与无穷(图十).
    每次应用 B B B时,向量都会根据 λ \lambda λ的值增长或收缩.
    在这里插入图片描述
    图九: v v v是矩阵 B B B特征值为 − 0.5 -0.5 0.5的一个特征向量。当 i i i增大, B i v B^iv Biv收敛为0.

    在这里插入图片描述
    图十: v v v对应的特征值为2,当 i i i增大, B i v B^iv Biv趋近与无穷。

    如果 B B B是对称矩阵,存在一组 n n n维线性无关的特征向量 v 1 , v 2 , . . . , v n . v_1,v_2,...,v_n. v1,v2,...,vn.这个集合不是唯一的,因为每个特征向量可以被一个任意的非零常数缩放.每个特征向量有一个对应的特征值 λ 1 , λ 2 , . . . , λ n \lambda_1,\lambda_2,...,\lambda_n λ1,λ2,...,λn.对于给定的矩阵,这些是唯一定义的。特征值可能相等,也可能不相等;例如,单位矩阵 I I I的特征值都是1,每一个非零向量都是 I I I的一个特征向量.

    如果 B B B作用于一个不是特征向量的向量呢?理解线性代数的一个非常重要的技巧——本节要教的技巧——是把一个向量看作是一个其行为被理解的其他向量的和。考虑 R n R^n Rn的一组基 { v i } \lbrace v_i\rbrace {vi}(因为对称矩阵 B B B n n n个线性无关的特征向量).任意 n n n维向量可以表示为这些特征向量的线性组合,由于矩阵-向量乘法是分配的,我们可以分别考察 B B B对每个特征向量的影响。

    在图11中,向量 x x x表示为两个特征向量 v 1 v_1 v1 v 2 v_2 v2的和.将 B B B作用于 x x x,并对结果求和。重复作用得到 B i x = B i v 1 + B i v 2 = λ 1 i v 1 + λ 2 i v 2 B^ix=B^iv_1+B^iv_2=\lambda^i_1v_1+\lambda^i_2v_2 Bix=Biv1+Biv2=λ1iv1+λ2iv2.如果所有特征值的大小都小于1, B i x B^ix Bix将收敛于0,因为构成 x x x的特征向量在重复作用于 B B B时收敛于0。如果有一个特征值的大于1, x x x会发散到无穷远.这就是为什么数值分析重视矩阵的谱半径: ρ = m a x ∣ λ i ∣ , λ 是 A 的 特 征 值 . \rho=max|\lambda_i| ,\lambda 是A的特征值. ρ=maxλi,λA.如果要 x x x很快的趋近于0, ρ ( B ) \rho(B) ρ(B)应该远小于1.在这里插入图片描述
    图11: 向量 x x x(实线箭头)可以表示为特征向量的线性组合(虚线箭头),特征值分别是 λ 1 = 0.7 \lambda_1=0.7 λ1=0.7 λ 2 = − 2. \lambda_2=-2. λ2=2.反复作用B对x的影响最好通过检查B对每个特征向量的影响来理解。当 B B B反复作用,一个特征向量收敛于零,而其他是发散的;则 B i x B^ix Bix是发散的。

    值得一提的是存在非对称矩阵族它没有一组完整的 n n n个独立的特征向量。这些矩阵被称为“缺陷矩阵”,这个名字背叛了它们在线性代数受挫时所应得的敌意。本文所描述的细节过于复杂,但缺陷矩阵的行为可以用广义特征向量和广义特征值来分析。当且仅当所有的广义特征值的量值都小于1时, B i x B^ix Bix收敛于0的规则仍然成立,但较难证明。

    这里有一个有用的事实:一个正的定义矩阵的特征值都是正的。这个事实可以从特征值的定义中得到证明: B v = λ v Bv=\lambda v Bv=λv v T B v = λ v T v v^TBv=\lambda v^Tv vTBv=λvTv根据正定的定义,对于任意非零向量 v v v, v T B v v^TBv vTBv 是正的,因此, λ \lambda λ肯定是正的。

    5.2 Jacobi iterations

    使用雅可比方法去求解 A x = b . Ax=b. Ax=b. 矩阵 A A A分成两部分: D D D,对角线为 A A A的对角元,其他非对角线的都为零;E,他的对角元全为零,非对角元的元素全为 A A A的非对角元。因此 A = D + E . A=D+E. A=D+E.得到雅可比法: A x = b Ax=b Ax=b D x = − E x + b Dx=-Ex+b Dx=Ex+b x = − D − 1 E x + D − 1 b x=-D^{-1}Ex+D^{-1}b x=D1Ex+D1b x = B x + z       w h e r e   B = − D − 1 E , z = D − 1 b x=Bx+z ~~~~~where ~ B = -D^{-1}E,z=D^{-1}b x=Bx+z     where B=D1E,z=D1b因为 D D D对角矩阵,很容易求逆。这个恒等式可以通过形成递归式转化为迭代法: x ( i + 1 ) = B x ( i ) + z x_{(i+1)}=Bx_{(i)}+z x(i+1)=Bx(i)+z给初始向量 x ( 0 ) x_{(0)} x(0),这个公式产生了一个向量序列,我们希望每一个连续的向量都比上一个向量更接近解 x x x x x x叫做方程 x = B x + z       w h e r e   B = − D − 1 E , z = D − 1 b x=Bx+z ~~~~~where ~ B = -D^{-1}E,z=D^{-1}b x=Bx+z     where B=D1E,z=D1b的稳定点,因为如果 x ( i ) = x x_{(i)}=x x(i)=x, x ( i + 1 ) x_{(i+1)} x(i+1)也会等于 x x x.

    这个推导看起来很随意,你是对的。对于 x x x替代方程 x = B x + z       w h e r e   B = − D − 1 E , z = D − 1 b x=Bx+z ~~~~~where ~ B = -D^{-1}E,z=D^{-1}b x=Bx+z     where B=D1E,z=D1b,我们可以形成任意数量的恒等式。事实上,通过分裂 A A A不同的是,通过选择一个不同的 D D D E E E——我们可以推导出高斯赛德尔法,或者连续超松弛(SOR).我们的希望是我们选择了一个 B B B的光谱半径较小的分裂。这里,为了简单,我任意选择了雅可比分解。

    假设我们从任意向量 x ( 0 ) x_{(0)} x(0)开始.对于每一次迭代,我们作用 B B B到它的向量,然后再加 z z z.每个迭代做什么?

    再一次,应用将一个向量看作是另一个被广泛理解的向量的和的原则。表达每个迭代 x ( i ) x_{(i)} x(i)看作精确解 x x x和误差项 e ( i ) e_{(i)} e(i)的和。则方程 x ( i + 1 ) = B x ( i ) + z x_{(i+1)}=Bx_{(i)}+z x(i+1)=Bx(i)+z变成: x ( i + 1 ) = B x ( i ) + z x_{(i+1)}=Bx_{(i)}+z x(i+1)=Bx(i)+z = B ( x + e ( i ) ) + z =B(x+e_{(i)})+z =B(x+e(i))+z = B x + z + B e ( i ) =Bx+z+Be_{(i)} =Bx+z+Be(i) = x + B e ( i )     b y   x = B x + z =x+Be_{(i)}~~~by~x=Bx+z =x+Be(i)   by x=Bx+z e ( i + 1 ) = B e ( i ) e_{(i+1)}=Be_{(i)} e(i+1)=Be(i)
    每次迭代不影响 x ( i ) x_{(i)} x(i)的”正确部分“(因为 x x x是稳定点);但是每次迭代都会影响误差项。从方程 e ( i + 1 ) = B e ( i ) e_{(i+1)}=Be_{(i)} e(i+1)=Be(i)中很显然知道,如果 ρ ( B ) < 1 \rho_{(B)}<1 ρ(B)<1,当 i i i趋近于无穷时,误差项 e ( i ) e_{(i)} e(i)接近于0.因此,初始向量 x ( 0 ) x_{(0)} x(0)对必然结果没有影响.

    当然, x ( 0 ) x_{(0)} x(0)的选择影响收敛到 x x x附近所需的迭代的次数.然而,他并没有谱半径 ρ ( B ) \rho(B) ρ(B)的影响大,谱半径决定收敛速度。假设 v j v_j vj是矩阵 B B B最大特征值的特征向量(即 ρ ( B ) = λ ( j ) \rho(B)=\lambda_{(j)} ρ(B)=λ(j)).初始误差 e ( 0 ) e_{(0)} e(0)表示为特征向量的线性组合,在 v j v_j vj方向上,包含了一个分量,这个分量是收敛最慢的。
    在这里插入图片描述
    图12: 矩阵 A A A的特征向量是沿着由二次型 f ( x ) f(x) f(x)定义的抛物面的坐标轴方向.每个特征向量用它的相关特征值来标记。每个特征值与相应斜率的陡度成正比。

    B B B一般不是对称的(即使 A A A是),甚至可能是有缺陷的。然而雅可比的收敛速度是主要依靠 ρ ( B ) \rho(B) ρ(B),这也依靠 A . A. A.不幸的是,雅可比不对每个矩阵 A A A收敛,甚至是正定矩阵 A . A. A.

    6.Convergence Analysis of Steepest Descent

    6.1 Instant Results

    为了解最陡下降的收敛性,第一步考虑特征向量 e ( i ) e_{(i)} e(i)的特征值 λ e \lambda_e λe.然后,残量 r ( i ) = − A e ( i ) = − λ e e ( i ) r_{(i)}=-Ae_{(i)}=-\lambda_ee_{(i)} r(i)=Ae(i)=λee(i)也是一个特征向量。通过方程 x ( i + 1 ) = x ( i ) + α ( i ) r ( i ) x_{(i+1)}=x_{(i)}+\alpha_{(i)}r_{(i)} x(i+1)=x(i)+α(i)r(i)得到: e ( i + 1 ) = e ( i ) + r ( i ) T r ( i ) r ( i ) T A r ( i ) r ( i )                       = e ( i ) + r ( i ) T r ( i ) λ e r ( i ) T r ( i ) ( − λ e e ( i ) ) = 0. e_{(i+1)}=e_{(i)}+\frac{r_{(i)}^Tr_{(i)}}{r_{(i)}^TAr_{(i)}}r_{(i)} \\ ~~~~~~~~~~~~~~~~~~~~~=e_{(i)}+\frac{r_{(i)}^Tr_{(i)}}{\lambda_er_{(i)}^Tr_{(i)}}(-\lambda_ee_{(i)}) \\ =0. e(i+1)=e(i)+r(i)TAr(i)r(i)Tr(i)r(i)                     =e(i)+λer(i)Tr(i)r(i)Tr(i)(λee(i))=0.

    图14演示了为什么只需一步就可以收敛到精确的解决方案。点 x ( i ) x_{(i)} x(i)位于椭球体的一个轴上,所以剩余的点直接指向椭球体的中心.选择 α ( i ) = λ e − 1 \alpha_{(i)}=\lambda_e^{-1} α(i)=λe1得到了瞬时收敛。

    以作更全面的分析,我们必须将 e ( i ) e_{(i)} e(i)作为特征向量的线性组合,我们进一步要求这些特征向量是正交的,可以知道如果 A A A是对称矩阵,则 A A A可以相似对角化,则 A A A存在 n n n维线性无关的特征向量。
    在这里插入图片描述
    图13: 雅可比方法的收敛性(a) B B B特征向量及其对应的特征值。不同于 A A A的特征向量,这些特征向量不是抛物面的坐标轴。(b)雅可比法开始于 [ − 2 , − 2 ] T [-2,-2]^T [2,2]T收敛于 [ 2 , − 2 ] T [2,-2]^T [2,2]T.。(c,d,e)误差 e ( 0 ) , e ( 1 ) , e ( 2 ) e_{(0)},e_{(1)},e_{(2)} e(0),e(1),e(2)(实箭头)及其特征向量分量(虚线箭头)。(f)箭头表示前四个错误向量的特征向量分量。误差的每个特征向量分量都以基于其特征值的期望速率收敛于零。
    在这里插入图片描述
    图14: 如果误差项是一个特征向量,最陡下降在第一次迭代时收敛于精确解。

    由于我们可以任意缩放特征向量,让我们选择使得每个特征向量都是单位长度的。这个选择给了我们有用的性质。
    v j T v k = { 1 ,     j = k 0 ,     j ≠ k v_j^Tv_k=\left\{ \begin{aligned} 1 ,& ~~~ j=k\\ 0, & ~~~j\neq k \end{aligned} \right. vjTvk={1,0,   j=k   j=k把误差项表示成特征向量的线性组合: e ( i ) = ∑ j = 1 n ξ j v j , e_{(i)}=\sum_{j=1}^{n}\xi_jv_j, e(i)=j=1nξjvj,这里 ξ j \xi_j ξj e ( i ) e_{(i)} e(i)的每个分量的长度,从上上式和上式我们得到以下特征:
    r ( i ) = − A e ( i ) = − ∑ j ξ j λ j v j , r_{(i)}=-Ae_{(i)}=-\sum_j\xi_j\lambda_jv_j, r(i)=Ae(i)=jξjλjvj, ∣ ∣ e ( i ) ∣ ∣ 2 = e ( i ) T e ( i ) = ∑ j ξ j 2 , ||e_{(i)}||^2=e_{(i)}^Te_{(i)}=\sum_j\xi_j^2, e(i)2=e(i)Te(i)=jξj2, e ( i ) T A e ( i ) = ( ∑ i ξ j v j T ) ( ∑ j ξ j λ j v j ) = ∑ j ξ j 2 λ j , e_{(i)}^TAe_{(i)}=(\sum_{i}\xi_jv_j^T)(\sum_j\xi_j\lambda_jv_j)=\sum_j\xi_j^2\lambda_j, e(i)TAe(i)=(iξjvjT)(jξjλjvj)=jξj2λj, ∣ ∣ r ( i ) ∣ ∣ 2 = r ( i ) T r ( i ) = ∑ j ξ j 2 λ j 2 , ||r_{(i)}||^2=r_{(i)}^Tr_{(i)}=\sum_j\xi_j^2\lambda_j^2, r(i)2=r(i)Tr(i)=jξj2λj2, r ( i ) T A r ( i ) = ∑ j ξ j 2 λ j 3 . r_{(i)}^TAr_{(i)}=\sum_j\xi_j^2\lambda_j^3. r(i)TAr(i)=jξj2λj3.
    在这里插入图片描述
    图15: 当特征值相等时,最陡下降在第一次迭代收敛于精确解。

    从方程 r ( i ) = − A e ( i ) = − ∑ j ξ j λ j v j r_{(i)}=-Ae_{(i)}=-\sum_j\xi_j\lambda_jv_j r(i)=Ae(i)=jξjλjvj可以看出 r ( i ) r_{(i)} r(i)也可以表示为特征向量分量的和以及这些分量的长度是 − ξ j λ j . -\xi_j\lambda_j. ξjλj.方程 ∣ ∣ e ( i ) ∣ ∣ 2 = e ( i ) T e ( i ) = ∑ j ξ j 2 , ||e_{(i)}||^2=e_{(i)}^Te_{(i)}=\sum_j\xi_j^2, e(i)2=e(i)Te(i)=jξj2, ∣ ∣ r ( i ) ∣ ∣ 2 = r ( i ) T r ( i ) = ∑ j ξ j 2 λ j 2 , ||r_{(i)}||^2=r_{(i)}^Tr_{(i)}=\sum_j\xi_j^2\lambda_j^2, r(i)2=r(i)Tr(i)=jξj2λj2,就是毕达哥拉斯定律(Pythagoras’ Law)。

    现在我们可以继续分析了。从 x ( i + 1 ) = x ( i ) + α ( i ) r ( i ) x_{(i+1)}=x_{(i)}+\alpha_{(i)}r_{(i)} x(i+1)=x(i)+α(i)r(i)得到: e ( i + 1 ) = e ( i ) + r ( i ) T r ( i ) r ( i ) T A r ( i ) r ( i ) e_{(i+1)}=e_{(i)}+\frac{r_{(i)}^Tr_{(i)}}{r_{(i)}^TAr_{(i)}}r_{(i)} e(i+1)=e(i)+r(i)TAr(i)r(i)Tr(i)r(i) = e ( i ) + ∑ j ξ j 2 λ j 2 ∑ j ξ j 2 λ j 3 r ( i ) =e_{(i)}+\frac{\sum_j\xi_j^2\lambda_j^2}{\sum_j\xi_j^2\lambda_j^3}r_{(i)} =e(i)+jξj2λj3jξj2λj2r(i)

    我们看最后一个例子,如果 e ( i ) e_{(i)} e(i)只有一个特征向量组合,则选择 α ( i ) = λ e − 1 , \alpha_{(i)}=\lambda_e^{-1}, α(i)=λe1然后一步就可以实现收敛。当 e ( i ) e_{(i)} e(i)是任意的,但所有特征向量有共同的的特征值 λ . \lambda. λ.由上面的方程得到: e ( i + 1 ) = e ( i ) + λ 2 ∑ j ξ j 2 λ 3 ∑ j ξ j 2 ( − λ e ( i ) ) = 0 e_{(i+1)}=e_{(i)}+\frac{\lambda^2\sum_j\xi_j^2}{\lambda^3\sum_j\xi_j^2}(-\lambda e_{(i)})=0 e(i+1)=e(i)+λ3jξj2λ2jξj2(λe(i))=0图15再次说明了为什么会出现瞬时收敛.因为所有特征值是相等的椭球体为球形;因此,无论我们从哪一点开始,残量必须指向球体的中心。正如前面所说,选择 α ( i ) = λ − 1 . \alpha_{(i)}=\lambda^{-1}. α(i)=λ1.

    在这里插入图片描述

    图16: 这两个向量的能量模量相等。

    然而,如果有几个不相等的,非零的特征值,那么没有 α ( i ) \alpha_{(i)} α(i)的选择会消除所有的特征向量分量,我们的选择变成了一种妥协。事实上, e ( i + 1 ) = e ( i ) + ∑ j ξ j 2 λ j 2 ∑ j ξ j 2 λ j 3 r ( i ) e_{(i+1)}=e_{(i)}+\frac{\sum_j\xi_j^2\lambda_j^2}{\sum_j\xi_j^2\lambda_j^3}r_{(i)} e(i+1)=e(i)+jξj2λj3jξj2λj2r(i)式中的分数最好被认为是 λ j − 1 \lambda_j^{-1} λj1这些值的加权平均值.权重 ξ j 2 \xi_j^2 ξj2确保 e ( i ) e_{(i)} e(i)长的组件是给定优先级。因此,在任何给定的迭代中, e ( i ) e_{(i)} e(i)的一些较短的组件实际上可能会增加长度(尽管不会很长)。因此,最陡下降法和共轭梯度法称为粗法。因此,最陡下降法和共轭梯度法称为粗法。相比之下,雅可比矩阵法更平滑,因为每个特征向量分量在每次迭代中都被约简。最陡下降和共轭梯度并不是平滑的,尽管它们在数学文献中经常被错误地认为是平滑的。

    6.2 General Convergence

    为了约束一般情况下最陡下降的收敛性,我们将定义能量范数 ∣ ∣ e ∣ ∣ A = ( e T A e ) 1 / 2 ||e||_A=(e^TAe)^{1/2} eA=(eTAe)1/2(见图16).这个范数比欧几里德范数更容易处理,而且在某种意义上是更自然的范数; 对方程 f ( p ) = f ( x ) + 1 2 ( p − x ) T A ( p − x ) f(p)=f(x)+\frac{1}{2}(p-x)^TA(p-x) f(p)=f(x)+21(px)TA(px)的检验表明极小化 ∣ ∣ e ( i ) ∣ ∣ A ||e_{(i)}||_A e(i)A等于极小化 f ( x ( i ) ) f(x_{(i)}) f(x(i)),由能量范数得到: ∣ ∣ e ( i + 1 ) ∣ ∣ A 2 = e ( i + 1 ) T A e ( i + 1 ) ||e_{(i+1)}||^2_A=e_{(i+1)}^TAe_{(i+1)} e(i+1)A2=e(i+1)TAe(i+1) = ( e ( i ) T + α ( i ) r ( i ) T ) A ( e ( i ) + α ( i ) r ( i ) ) =(e_{(i)}^T+\alpha_{(i)}r_{(i)^T)}A(e_{(i)}+\alpha_{(i)}r_{(i)}) =(e(i)T+α(i)r(i)T)A(e(i)+α(i)r(i)) = e ( i ) T A e ( i ) + 2 α ( i ) r ( i ) T A e ( i ) + α ( i ) 2 r ( i ) T A r ( i )     ( A 的 对 称 性 ) =e_{(i)}^TAe_{(i)}+2\alpha_{(i)}r_{(i)}^TAe_{(i)}+\alpha_{(i)}^2r_{(i)}^TAr_{(i)}~~~(A的对称性) =e(i)TAe(i)+2α(i)r(i)TAe(i)+α(i)2r(i)TAr(i)   (A) ∣ ∣ e ( i ) ∣ ∣ A 2 + 2 r ( i ) T r ( i ) r ( i ) A r ( i ) ( − r ( i ) T r ( i ) ) + ( r ( i ) T r ( i ) r ( i ) T A r ( i ) ) 2 r ( i ) T A r ( i ) ||e_{(i)}||^2_A +2\frac{r_{(i)}^Tr_{(i)}}{r_{(i)}Ar_{(i)}}(-r_{(i)}^Tr_{(i)})+(\frac{r_{(i)}^Tr_{(i)}}{r_{(i)}^TAr_{(i)}})^2r_{(i)}^TAr_{(i)} e(i)A2+2r(i)Ar(i)r(i)Tr(i)(r(i)Tr(i))+(r(i)TAr(i)r(i)Tr(i))2r(i)TAr(i) = ∣ ∣ e ( i ) ∣ ∣ A 2 − ( r ( i ) T r ( i ) ) 2 r ( i ) A r ( i ) =||e_{(i)}||^2_A-\frac{(r_{(i)}^Tr_{(i)})^2}{r_{(i)}Ar_{(i)}} =e(i)A2r(i)Ar(i)(r(i)Tr(i))2 = ∣ ∣ e ( i ) ∣ ∣ A 2 ( 1 − ( r ( i ) T r ( i ) ) 2 ( r ( i ) T A r ( i ) ( e ( i ) T A e ( i ) ) ) ) =||e_{(i)}||^2_A(1-\frac{(r_{(i)}^Tr_{(i)})^2}{(r_{(i)}^TAr_{(i)}(e_{(i)}^TAe_{(i)}))}) =e(i)A2(1(r(i)TAr(i)(e(i)TAe(i)))(r(i)Tr(i))2) = ∣ ∣ e ( i ) ∣ ∣ A 2 ( 1 − ( ∑ j ξ j 2 λ j 2 ) 2 ( ∑ j ξ j 2 λ j 3 ) ( ∑ j ξ j 2 λ j ) ) ( 由 e ( i ) T A e ( i ) = ( ∑ i ξ j v j T ) ( ∑ j ξ j λ j v j ) = ∑ j ξ j 2 λ j , ∣ ∣ r ( i ) ∣ ∣ 2 = r ( i ) T r ( i ) = ∑ j ξ j 2 λ j 2 , r ( i ) T A r ( i ) = ∑ j ξ j 2 λ j 3 . ) =||e_{(i)}||^2_A(1-\frac{(\sum_j\xi_j^2\lambda_j^2)^2}{(\sum_j\xi_j^2\lambda_j^3)(\sum_j\xi_j^2\lambda_j)})\\(由e_{(i)}^TAe_{(i)}=(\sum_{i}\xi_jv_j^T)(\sum_j\xi_j\lambda_jv_j)=\sum_j\xi_j^2\lambda_j,\\||r_{(i)}||^2=r_{(i)}^Tr_{(i)}=\sum_j\xi_j^2\lambda_j^2,\\ r_{(i)}^TAr_{(i)}=\sum_j\xi_j^2\lambda_j^3.) =e(i)A2(1(jξj2λj3)(jξj2λj)(jξj2λj2)2)e(i)TAe(i)=(iξjvjT)(jξjλjvj)=jξj2λj,r(i)2=r(i)Tr(i)=jξj2λj2,r(i)TAr(i)=jξj2λj3. = ∣ ∣ e ( i ) ∣ ∣ A 2 w 2 ,       w 2 = 1 − ( ∑ j ξ j 2 λ j 2 ) 2 ( ∑ j ξ j 2 λ j 3 ) ( ∑ j ξ j 2 λ j ) =||e_{(i)}||_A^2w^2,~~~~~w^2=1-\frac{(\sum_j\xi_j^2\lambda_j^2)^2}{(\sum_j\xi_j^2\lambda_j^3)(\sum_j\xi_j^2\lambda_j)} =e(i)A2w2,     w2=1(jξj2λj3)(jξj2λj)(jξj2λj2)2
    分析依赖于找到一个 w w w上界.为展示权重和特征值影响收敛性,我们应该可以从 n = 2 n=2 n=2得到一个结果。假设 λ 1 ≥ λ 2 . \lambda_1\ge\lambda_2. λ1λ2. A A A的谱条件数定义为 κ = λ 1 / λ 2 ≥ 1. \kappa=\lambda_1/\lambda_2\ge1. κ=λ1/λ21. e ( i ) e_{(i)} e(i)的斜率(相对于由特征向量定义的坐标系),取决于起点,是表示 μ = ξ 2 / ξ 1 \mu=\xi_2/\xi_1 μ=ξ2/ξ1,有: w 2 = 1 − ( ∑ j ξ j 2 λ j 2 ) 2 ( ∑ j ξ j 2 λ j 3 ) ( ∑ j ξ j 2 λ j ) w^2=1-\frac{(\sum_j\xi_j^2\lambda_j^2)^2}{(\sum_j\xi_j^2\lambda_j^3)(\sum_j\xi_j^2\lambda_j)} w2=1(jξj2λj3)(jξj2λj)(jξj2λj2)2 = 1 − ( κ 2 + μ 2 ) 2 ( κ + μ 2 ) ( κ 3 + μ 2 ) . =1-\frac{(\kappa^2+\mu^2)^2}{(\kappa+\mu^2)(\kappa^3+\mu^2)}. =1(κ+μ2)(κ3+μ2)(κ2+μ2)2.
    w w w的值决定最陡下降的收敛速度, μ \mu μ κ \kappa κ函数用图17表示。这个图表证实了我的两个例子。如果 e ( 0 ) e_{(0)} e(0)是特征向量,则斜率 μ \mu μ是零(或无穷);从图可以看出 w w w等于零时,收敛瞬间的。如果特征值相等,那么条件数 κ = 1 \kappa=1 κ=1;我们再次看到 w w w等于0.图18展示了图17四个角附近的示例。这些二次型在由它们的特征向量定义的坐标系中作图。图18(a)和18(b)是带有较大条件数的示例。如果选择一个幸运的起点,最陡的下降可以迅速收敛(图18(a)),但通常在 κ \kappa κ较大时最糟糕(图18(b)).后一个数字给了我们最好的直觉,为什么一个大的条件数可能是坏的: f ( x ) f(x) f(x)形成一个低谷,最陡的下降在低谷的两条线之间来回反弹,而几乎没有变化。在图18©和图18(d)中,条件数很小,所以二次形式近似于球面,无论起点如何,收敛速度都很快。
    在这里插入图片描述
    图17: 最陡下降的 w w w作为 μ \mu μ( e ( i ) e_{(i)} e(i)的斜率)和 κ \kappa κ( A A A的条件数)的函数的收敛性。当 μ \mu μ κ \kappa κ比较小时收敛是很快的。当 μ = ± κ \mu=\pm \kappa μ=±κ时,对于一个固定的矩阵,收敛性最差。

    在这里插入图片描述
    图18: 这四个例子表示图中对应的四个角附近的点
    图17。(a) κ \kappa κ大, μ \mu μ小.(b)收敛性差的一个例子, κ \kappa κ μ \mu μ都大.© κ \kappa κ μ \mu μ都小.(d)小 κ \kappa κ,大 μ \mu μ.

    在这里插入图片描述
    图19: 实线表示在最陡下降时收敛性最差的起始点。
    虚线表示收敛的步骤。如果第一次迭代从最坏的情况开始,那么所有后续的迭代都要这样做。每一步都与抛物面轴(灰色箭头)相交于45度角。这里 κ = 3.5. \kappa=3.5. κ=3.5.

    保持 κ \kappa κ常数(因为 A A A是固定的),当 μ = ± κ \mu=\pm\kappa μ=±κ时,简单的微积分运算就能得出方程 w 2 = 1 − ( κ 2 + μ 2 ) 2 ( κ + μ 2 ) ( κ 3 + μ 2 ) w^2=1-\frac{(\kappa^2+\mu^2)^2} {(\kappa+\mu^2)(\kappa^3+\mu^2)} w2=1(κ+μ2)(κ3+μ2)(κ2+μ2)2是最大化的.在图17中,人们可以看到这条线形成的模糊的山脊。图19显示了矩阵 A A A最坏情况的起始点。这些起始点落在由 ξ 2 / ξ 1 = ± κ \xi_2/\xi_1=\pm\kappa ξ2/ξ1=±κ定义的线上。对于 w w w的一个最大值(对应于最坏情况的起始点)是通过设置 μ 2 = κ 2 \mu^2=\kappa^2 μ2=κ2得到:
    w 2 ≤ 1 − 4 κ 4 κ 5 + 2 κ 4 + κ 3 w^2\le1-\frac{4\kappa^4}{\kappa^5+2\kappa^4+\kappa^3} w21κ5+2κ4+κ34κ4 = k 5 − 2 κ 4 + κ 3 κ 5 + 2 κ 4 + κ 3 =\frac{k^5-2\kappa^4+\kappa^3}{\kappa^5+2\kappa^4+\kappa^3} =κ5+2κ4+κ3k52κ4+κ3 = ( κ − 1 ) 2 ( κ + 1 ) 2 =\frac{(\kappa-1)^2}{(\kappa+1)^2} =(κ+1)2(κ1)2 w ≤ κ − 1 κ + 1 w\le \frac{\kappa-1}{\kappa+1} wκ+1κ1
    上式不等式绘制在图20中,条件越差的矩阵(即其条件数 κ \kappa κ越大),最速下降的收敛速度越慢.在9.2节中,证明了对于 n ≥ 2 n\ge 2 n2时,上式也是有效的,如果一个对称的正定矩阵的条件数被定义为 κ = λ m a x / λ m i n , \kappa=\lambda_{max}/\lambda_{min}, κ=λmax/λmin,最大特征值与最小特征值之比。最陡下降的收敛结果是 ∣ ∣ e ( i ) ∣ ∣ A ≤ ( κ − 1 κ + 1 ) i ∣ ∣ e ( 0 ) ∣ ∣ A ||e_{(i)}||_A\le(\frac{\kappa-1}{\kappa+1})^i||e_{(0)}||_A e(i)A(κ+1κ1)ie(0)A f ( x ( i ) ) − f ( x ) f ( x ( 0 ) ) − f ( x ) = 1 2 e ( i ) T A e ( i ) 1 2 e ( 0 ) T A e ( 0 ) \frac{f(x_{(i)})-f(x)}{f(x_{(0)})-f(x)}=\frac{\frac{1}{2}e_{(i)}^TAe_{(i)}}{\frac{1}{2}e_{(0)}^TAe_{(0)}} f(x(0))f(x)f(x(i))f(x)=21e(0)TAe(0)21e(i)TAe(i) ≤ ( κ − 1 κ + 1 ) 2 i . \le(\frac{\kappa-1}{\kappa+1})^{2i}. (κ+1κ1)2i.
    在这里插入图片描述
    图20: 最速(陡)下降的收敛性(每次迭代)随着矩阵的条件数的增加而恶化。

    7. The Method of Conjugate Directions

    7.1 Conjugacy

    如果陡坡下降,我就会像前面的步骤一样往下走(见图8)。如果我们每迈出一步,第一次就能成功不是更好吗?这里有一个想法:让我们选择一组正交搜索方向 d ( 0 ) , d ( 1 ) , . . . , d ( n − 1 ) . d_{(0)},d_{(1)},...,d_{(n-1)}. d(0),d(1),...,d(n1).每个搜索方向上,我们都要精确地走一步,这一步的长度刚好与 x x x相等。 n n n步之后,我们就完成了。

    图21说明了这个想法,使用坐标轴作为搜索方向。第一步(水平方向)指向正确的 x 1 x_1 x1-坐标;第二步(垂直)将击中要害。注意 e ( 1 ) e_{(1)} e(1)正交于 d ( 0 ) . d_{(0)}. d(0).一般来说,每一步我们选择一个点: x ( i + 1 ) = x ( i ) + α ( i ) d ( i ) x_{(i+1)}=x_{(i)}+\alpha_{(i)}d_{(i)} x(i+1)=x(i)+α(i)d(i)

    未找到 α ( i ) \alpha_{(i)} α(i)的值,利用 e ( i + 1 ) e_{(i+1)} e(i+1)正交于 d ( i ) , d_{(i)}, d(i),所以我们永远不需要再介入 d ( i ) d_{(i)} d(i)方向.利用这个条件得到: d ( i ) T e ( i + 1 ) = 0 d^T_{(i)}e_{(i+1)}=0 d(i)Te(i+1)=0 d ( i ) T ( e ( i ) + α ( i ) d ( i ) ) = 0 d_{(i)}^T(e_{(i)}+\alpha_{(i)}d_{(i)})=0 d(i)T(e(i)+α(i)d(i))=0 α ( i ) = − d ( i ) T e ( i ) d ( i ) T d ( i ) \alpha_{(i)}=-\frac{d_{(i)}^Te_{(i)}}{d_{(i)}^Td{(i)}} α(i)=d(i)Td(i)d(i)Te(i)
    如果不知道 e ( i ) , e_{(i)}, e(i),我们无法计算 α ( i ) \alpha_{(i)} α(i),这样就会导致不我们什么也没做成;如果我们知道 e ( i ) e_{(i)} e(i),则问题也就会解决了.

    解决的办法是使搜索方向 A A A-正交而不是正交。 d ( i ) d_{(i)} d(i) d ( j ) d_{(j)} d(j)两个向量都是 A A A-正交或者共轭,如果: d ( i ) T A d ( j ) = 0. d_{(i)}^TAd_{(j)}=0. d(i)TAd(j)=0.

    图22(a)显示了 A A A-正交向量的样子.想象一下,如果这篇文章是印在泡泡糖上的,你抓住图22(a)的末端,拉伸它,直到椭圆变成圆形。然后这些向量看起来是正交的,如图22(b)所示。

    我们新的需求是 e ( i + 1 ) e_{(i+1)} e(i+1) A A A-正交于 d ( i ) d_{(i)} d(i)(见图23(a)).在最陡的下降中,并非巧合的是,这个正交条件等价于在搜索方向 d ( i ) d_{(i)} d(i)上找到最小点,。为了看到这一点,将方向导数设为0: d d α f ( x ( i + 1 ) ) = 0 \frac{d}{d\alpha}f(x_{(i+1)})=0 dαdf(x(i+1))=0 f ′ ( x ( i + 1 ) ) T d d α x ( i + 1 ) = 0 f'(x_{(i+1)})^T\frac{d}{d\alpha}x_{(i+1)}=0 f(x(i+1))Tdαdx(i+1)=0 − r ( i + 1 ) T d ( i ) = 0 -r_{(i+1)}^Td_{(i)}=0 r(i+1)Td(i)=0 d ( i ) T A e ( i + 1 ) = 0. d_{(i)}^TAe_{(i+1)}=0. d(i)TAe(i+1)=0.

    根据方程 α ( i ) = − d ( i ) T e ( i ) d ( i ) T d ( i ) \alpha_{(i)}=-\frac{d_{(i)}^Te_{(i)}}{d_{(i)}^Td{(i)}} α(i)=d(i)Td(i)d(i)Te(i)的推导,当搜索方向是 A A A-正交时,下面就是关于 α ( i ) \alpha_{(i)} α(i)的表达式: α ( i ) = − d ( i ) T A e ( i ) d ( i ) T A d ( i ) \alpha_{(i)}=-\frac{d_{(i)}^TAe_{(i)}}{d_{(i)}^TAd_{(i)}} α(i)=d(i)TAd(i)d(i)TAe(i) = d ( i ) T r ( i ) d ( i ) T A d ( i ) =\frac{d_{(i)}^Tr_{(i)}}{d_{(i)}^TAd_{(i)}} =d(i)TAd(i)d(i)Tr(i)

    在这里插入图片描述
    图22: 这些向量对是 A A A-正交的…因为这些向量对是正交的。

    与方程 α ( i ) = − d ( i ) T e ( i ) d ( i ) T d ( i ) \alpha_{(i)}=-\frac{d_{(i)}^Te_{(i)}}{d_{(i)}^Td{(i)}} α(i)=d(i)Td(i)d(i)Te(i)不同,我们可以计算这个表达式。注意,如果搜索向量是残差,这个公式将与最速下降所用的公式相同。

    在这里插入图片描述
    图23: 共轭方向法在 n n n步内收敛。(a)第一步是沿着某个方向 d ( 0 ) d_{(0)} d(0)走.限制 e ( 1 ) e_{(1)} e(1)必须是 A A A-正交于 d ( 0 ) 下 选 择 最 小 值 点 d_{(0)}下选择最小值点 d(0)x_{(1)} 。 ( b ) 初 始 误 差 。(b)初始误差 (b)e_{(0)} 能 被 表 达 成 能被表达成 A$-正交组合的和(灰色箭头)。共轭方向的每一步都会消除其中一个分量。

    为了证明这个过程确实是在 n n n步中计算 x x x,将误差项表示为搜索方向的线性组合;即 e ( 0 ) = ∑ j = 0 n − 1 δ j d ( j ) . e_{(0)}=\sum_{j=0}^{n-1}\delta_jd_{(j)}. e(0)=j=0n1δjd(j).

    δ j \delta_j δj的值可以用数学技巧找到。因为搜索方向是 A A A-正交,通过对表达式进行预先乘以 d ( k ) T A d_{(k)}^TA d(k)TA,可以从上式中除去除一个之外的所有 δ j \delta_j δj值: d ( k ) T A e ( 0 ) = ∑ j δ ( j ) d ( k ) T A d ( j ) d_{(k)}^TAe_{(0)}=\sum_j\delta_{(j)}d_{(k)}^TAd_{(j)} d(k)TAe(0)=jδ(j)d(k)TAd(j) d ( k ) T A e ( 0 ) = δ ( k ) d ( k ) T A d ( k )     ( d 向 量 的 A 正 交 ) d_{(k)}^TAe_{(0)}=\delta_{(k)}d_{(k)}^TAd_{(k)}~~~(d向量的A正交) d(k)TAe(0)=δ(k)d(k)TAd(k)   (dA) δ ( k ) = d ( k ) T A e ( 0 ) d ( k ) T A d ( k ) \delta_{(k)}=\frac{d_{(k)}^TAe_{(0)}}{d_{(k)}^TAd_{(k)}} δ(k)=d(k)TAd(k)d(k)TAe(0) = d ( k ) T A ( e ( 0 ) + ∑ i = 0 k − 1 α ( i ) d ( i ) ) d ( k ) T A d ( k )     ( d 向 量 的 A 正 交 ) =\frac{d_{(k)}^TA(e_{(0)}+\sum_{i=0}^{k-1}\alpha_{(i)}d_{(i)})}{d_{(k)}^TAd_{(k)}}~~~(d向量的A正交) =d(k)TAd(k)d(k)TA(e(0)+i=0k1α(i)d(i))   (dA) = d ( k ) T A e ( k ) d ( k ) T A d ( k )     由 ( x ( i + 1 ) = x ( i ) + α ( i ) d ( i ) ) =\frac{d_{(k)}^TAe_{(k)}}{d_{(k)}^TAd_{(k)}}~~~由(x_{(i+1)}=x_{(i)}+\alpha_{(i)}d_{(i)}) =d(k)TAd(k)d(k)TAe(k)   (x(i+1)=x(i)+α(i)d(i))
    由上式方程和方程 α ( i ) = − d ( i ) T A e ( i ) d ( i ) T A d ( i ) \alpha_{(i)}=-\frac{d_{(i)}^TAe_{(i)}}{d_{(i)}^TAd_{(i)}} α(i)=d(i)TAd(i)d(i)TAe(i),我们发现 α ( i ) = − δ ( i ) . \alpha_{(i)}=-\delta_{(i)}. α(i)=δ(i).这个事实给了我们一种看待误差项的新方法。如下式所示,按组件构建 x x x组件的过程也可以看作是按组件减少错误项的过程(参见图23(b))。 e ( i ) = e ( 0 ) + ∑ j = 0 i − 1 α ( j ) d ( j ) e_{(i)}=e_{(0)}+\sum_{j=0}^{i-1}\alpha_{(j)}d_{(j)} e(i)=e(0)+j=0i1α(j)d(j) = ∑ j = 0 n − 1 δ ( j ) d ( j ) − ∑ j = 0 i − 1 δ ( j ) d ( j ) =\sum_{j=0}^{n-1}\delta_{(j)}d_{(j)}-\sum_{j=0}^{i-1}\delta_{(j)}d_{(j)} =j=0n1δ(j)d(j)j=0i1δ(j)d(j) = ∑ j = i n − 1 δ ( j ) d ( j ) . =\sum_{j=i}^{n-1}\delta_{(j)}d_{(j)}. =j=in1δ(j)d(j).经过 n n n次迭代,每个组件都被砍掉,以及 e ( n ) = 0 ; e_{(n)}=0; e(n)=0;证毕。

    7.2 Gram-Schmidt Conjugation

    现在只需要 A A A-正交 { d ( i ) } \{d_{(i)}\} {d(i)}搜索方向的一个集合。幸运的是,有一种生成它们的简单方法,称为共轭Gram-Schmidt过程。

    假设我们由一组 n n n维线性无关的向量 u 0 , u 1 , . . . , u n − 1 . u_0,u_1,...,u_{n-1}. u0,u1,...,un1.坐标轴会在紧要关头发挥作用,尽管可能会有更明智的选择。为了构造 d ( i ) d_{(i)} d(i),取 u i u_i ui并减去与之前 d d d向量不 A A A-正交的任何分量(见图24)。也就是说,令 d ( 0 ) = u 0 d_{(0)}=u_0 d(0)=u0,以及对 i > 0 i>0 i>0时,令 d ( i ) = u i + ∑ k = 0 i − 1 β i k d ( k ) d_{(i)}=u_i+\sum_{k=0}^{i-1}\beta_{ik}d_{(k)} d(i)=ui+k=0i1βikd(k)这里 β i k 是 定 义 在 \beta_{ik}是定义在 βik i > k i>k i>k的,为求出他的值,我们可以用于求 δ j \delta_j δj同样的技巧: d ( i ) T A d ( j ) = u i T A d ( j ) + ∑ k = 0 i − 1 β i k d ( k ) T A d ( j ) d_{(i)}^TAd_{(j)}=u_i^TAd_{(j)}+\sum_{k=0}^{i-1}\beta_{ik}d_{(k)}^TAd_{(j)} d(i)TAd(j)=uiTAd(j)+k=0i1βikd(k)TAd(j) 0 = u i T A d ( j ) + β i j d ( j ) T A d ( j ) ,     i > j 0=u_i^TAd_{(j)}+\beta_{ij}d_{(j)}^TAd_{(j)},~~~i>j 0=uiTAd(j)+βijd(j)TAd(j),   i>j β i j = − u i T A d ( j ) d ( j ) T A d ( j ) \beta_{ij}=-\frac{u_i^TAd_{(j)}}{d_{(j)}^TAd_{(j)}} βij=d(j)TAd(j)uiTAd(j)

    在共轭方向法中使用gramm - schmidt共轭的困难在于,必须将所有旧的搜索向量保存在内存中,以构造每个新向量,还有 O ( n 3 ) O(n^3) O(n3)操作。
    在这里插入图片描述
    图24: 两个向量的共轭Gram-Schmidt 。

    展开全文
  • 基于复变量求导法的共轭梯度法及其在边界条件辨识中的应用,崔苗,端维伟,为了利用共轭梯度法的计算精度高和收敛速度快的优点,避免传统共轭梯度法在求解非线性热传导反问题中的微分处理、复杂的推导过程
  • 优化算法-共轭梯度法

    2017-06-28 22:32:38
    共轭梯度法的matlab代码
  • 预处理共轭梯度法解病态问题及在GPS中的应用.pdf
  • 用matlab实现共轭梯度法求解实例.doc 用MATLAB 实现共轭梯度法求解实例 康福 201103710031 1.无约束优化方法 1.1 无约束优化方法的必要性 一般机械优化设计问题,都是在一定的限制条件下追求某一指标为最小,它 们...

    41528d3028836879cd698677c3999917.gif用matlab实现共轭梯度法求解实例.doc

    用MATLAB 实现共轭梯度法求解实例 康福 201103710031 1.无约束优化方法 1.1 无约束优化方法的必要性 一般机械优化设计问题,都是在一定的限制条件下追求某一指标为最小,它 们都属于约束优化问题。但是为什么要研究无约束优化问题?(1)有些实际问题,其数学模型本身就是一个无约束优化问题。(2)通过熟悉它的解法可以为研究约束优化问题打下良好的基础。(3)约束优化问题的求解可以通过一系列无约束优化方法来达到。所以无约束 优化问题的解法是优化设计方法的基本组成部分,也是优化方法的基础。 (4)对于多维无约束问题来说,古典极值理论中令一阶导数为零,但要求二阶 可微,且要判断海赛矩阵为正定才能求得极小点,这种方法有理论意义, 但无实用价值。和一维问题一样,若多元函数F(X)不可微,亦无法求解。 但古典极值理论是无约束优化方法发展的基础。 1.2共轭梯度法 目前已研究出很多种无约束优化方法,它们的主要不同点在于构造搜索方向 上的差别。 (1)间接法——要使用导数,如梯度法、 (阻尼)牛顿法、变尺度法、共轭梯度 法等。 (2)直接法——不使用导数信息,如坐标轮换法、鲍威尔法单纯形法等。 用直接法寻找极小点时,不必求函数的导数,只要计算目标函数值。这类方 法较适用于解决变量个数较少的(n ≤20)问题,一般情况下比间接法效率低。 间接法除要计算目标函数值外,还要计算目标函数的梯度,有的还要计算其海赛 矩阵。 搜索方向的构成问题乃是无约束优化方法的关键。 共轭梯度法是沿着共轭方向进行搜索,属于共轭方向法中的一种,该方法中 每一个共轭向量都是依赖于迭代点处的负梯度而构造出来。共轭梯度法作为一种 实用的迭代法,它主要有下面的优点: (1)算法中,系数矩阵A的作用仅仅是用来由已知向量P产生向量W=AP,这不仅 可充分利用A的稀疏性,而且对某些提供矩阵A较为困难而由已知向量P产 生向量W=AP又十分方便的应用问题是很有益的。 (2)不需要预先估计任何参数就可以计算,这一点不像SOR等; (3)每次迭代所需的计算,主要是向量之间的运算,便于并行化。 共轭梯度法原理的知识较多,请详见《机械优化设计》第四章的第四、五节。 图1为共轭梯度法的程度框图 1 ( 0,1,2, ) k k k k s k       x x图1为共轭梯度法的程度框图 2.设计题目及要求 2.1设计题目 用共轭梯度法求二次函数 2 2 1 2 1 2 1 1 2 ( , ) 2 4 2 f x x x x x x x    的极小点及极小值。 2.2设计要求 (1)使用matlab编写程序,熟练撑握matlab编程方法。(2)学习并撑握共轭梯度法的原理、方法及应用,并了解不同无约束优化方法的 区别、优缺点及特殊要求。 (3)编写程序,计算出二次函数的极小点及极小值,并适当选取不同的初始点及 迭代精度精度,分析比较结果。 三.计算步骤 3.1计算求解 解:已知初始点[1,1] T迭代精度 0.001   1)第一次沿负梯度方向搜寻 计算初始点处的梯度: 为一维搜索最佳步长,应满足 得:2)第二次迭代代入目标函数 由 得 从而有: 因 收敛。 0 1 2 0 2 1 2 2 4 4 ( ) 4 2 2 x x f x x                    x x 0 1 0 0 0 0 0 1 4 1 4 1 2 1 2                               x x d 1 0 0 2 ( ) min ( ) min(40 20 3) f f           x x d 0 0.25   1 2 0.5        x 1 1 ( ) 2 f           x 2 1 2 0 0 ( ) 5 0.25 20 ( ) f f       x x 1 1 0 0 2 ( ) 1.5 f             d x d 2 1 1 2 2 2 2 0.5 1.5 0.5 1.5                              x x d 2 2 ( ) (2 2 ) 2(0.5 1.5 ) 2(2 2 )(0.5 1.5 ) 4(2 2 ) ( ) f x                  ( ) 0     1   2 2 2 4 0 , ( ) 8, ( ) 2 0 f f                  x x x 2 ( ) 0 f     x3.2运行与程序 运行:打开matlab,确定conjugate_grad_2d.m文件夹为当前目录。在命令窗中输入:f=conjugate_grad_2d([1,1],0.001)选择不同的初始点坐标[0,0],[0,1],[1,0],和迭代精度0.01,0.0001,进 行运行时,需要多次调用conjugate_grad_2d函数。 程序及说明: function f=conjugate_grad_2d(x0,t) %用共轭梯度法求已知函数f(x1,x2)=x1^2+2*x2^2-4*x1-2*x1*x2的极值点 %已知初始点坐标:x0 %已知收敛精度:t %求得已知函数的极值:f x=x0; syms xi yi a; %定义自变量,步长为符号变量 f=xi^2+2*yi^2-4*xi-2*xi*yi; %创建符号表达式f fx=diff(f,xi); %求表达式f对xi的一阶求导 fy=diff(f,yi); %求表达式f对yi的一阶求导 fx=subs(fx,{xi,yi},x0); %代入初始点坐标计算对xi的一阶求导实值 fy=subs(fy,{xi,yi},x0); %代入初始点坐标计算对yi的一阶求导实值 fi=[fx,fy]; %初始点梯度向量 count=0; %搜索次数初始为0 while double(sqrt(fx^2+fy^2))>t %搜索精度不满足已知条件s=-fi; %第一次搜索的方向为负梯度方向if count<=0s=-fi;elses=s1;endx=x+a*s; %进行一次搜索后的点坐标f=subs(f,{xi,yi},x); %构造一元搜索的一元函数φ(a)f1=diff(f); %对函数φ(a)进行求导f1=solve(f1); %得到最佳步长aif f1~=0ai=double(f1); %强制转换数据类型为双精度数值elsebreak %若a=0,则直接跳出循环,此点即为极值点endx=subs(x,a,ai); %得到一次搜索后的点坐标值f=xi^2+2*yi^2-4*xi-2*xi*yi;fxi=diff(f,xi);fyi=diff(f,yi);fxi=subs(fxi,{xi,

    展开全文
  • 针对目前焦炭质量预测线性方法的现状及不足,建立了共轭梯度法优化的BP神经网络焦炭质量预测...模型的实例验证分析表明,共轭梯度法优化的BP神经网络焦炭质量预测模型有较好的适应性和预测精度,具有一定的实际应用价值。
  • 一个非线性共轭梯度法C++源码,含一个简单的非线性方程组,可以应用扩展到M*N情况下。
  • 共轭梯度法

    千次阅读 2015-06-13 17:58:21
    概述作为一种迭代的优化方法,共轭梯度(Conjugate Gradient,cg)由Hestenes和Stiefe于1951年提出。cg是针对形如式(1-1)的优化方法。 Ax=b(1-1) Ax=b \tag{1-1}\label{1-1} 需要指出的是,式(1-1)有着广泛的应用...

    概述

    作为一种迭代的优化方法,共轭梯度(Conjugate Gradient,cg)由Hestenes和Stiefe于1951年提出。cg是针对形如式(1-1)的优化方法。

    Ax=b(1-1)

    需要指出的是,式(1-1)有着广泛的应用场景。例如,令 x 为二次问题(如式(1-2)所示)的最小值:
    f(x)=12xTAxbx(1-2)

    x=argminxf(x)(1-3)