数值分析03 木灵的炼金工作室

残量与误差控制

我们继续进行误差分析的讨论:

为表征线性方程组$Ax = b$之解$x^$的相对误差,记$x^$的残量$r=b-Ax^*$,则相对误差有如下关系:

\[\frac{\Vert x-x^* \Vert}{\Vert x \Vert} \leqslant cond(A)\frac{\Vert r \Vert}{\Vert b \Vert}\tag{3.1}\]

证明借用了矩阵误差的相容性条件,推导不难,请见讲义.

式(3.1)表明,即使求解出的残量很小,但如果方程的状态数较大,仍可能会导出一个较大的相对误差,这表明$cond(A)$导致的误差增大是一个不可回避的问题.

残量方程

根据残量的定义,有第$n$次迭代的残量:$b-Ax_n=r_n$,解方程$Ad_n=r_n$可以解出残量$r_n$反射至$x_n$的数值,迭代:$x_{n+1}=x_n+d_n$即可减小$x_{n+1}$的残量.

超定线性方程组的最小二乘解

看一个引例:对物理量$x$执行三次测量,结果分别是$x_1, x_2, x_3$,则这个真实量最有可能是? 记函数$f(x) = \frac{1}{3}[(x-x_1)^2+(x-x_2)^2+(x-x_3)^2]$,称其为方差函数. 我们可以后验地认为,当方差函数取最小值时的$x=(x_1+x_2+x_3)/3$是最有可能取到的. 这样的思想被称为“极大似然”,这种方法被称为“最小二乘”. 我们将上述的情况转化为矩阵形式:

\[\begin{pmatrix} 1\\ 1\\ 1\\ \end{pmatrix} \begin{pmatrix} x \end{pmatrix} = \begin{pmatrix} x_1\\ x_2\\ x_3\\ \end{pmatrix}\]

上式看作是一个线性方程组,而若$x_1, x_2, x_3$不完全相同,则上述方程组的增广矩阵之秩比系数矩阵秩大,故这样的方程组属”无解”方程组,称为超定线性方程组. 但我们仍可以通过某一优化目标来将这个问题转化为一个最优化问题. 当这一优化目标是”均方误差最小化”,或称为”2-范数最小化”时,这样的解法称为”超定线性方程组的最小二乘解法”.

关于超定方程最小二乘解的导函数推导需要用到矩阵分析领域的矩阵求导术,较为复杂且不易理解,详细内容请参见课本或讲义. 下面提出一个较为容易理解的推导方法:

解的推导:几何理解

对于线性方程组$Ax =b$,将$A$列向量记为$(a_1, a_2, a_3, …, a_n)$,其中$a_i\in \mathbb{R}^{n\times1}$,考虑优化条件:

\[\min_{x^*}\Vert b - Ax^* \Vert_2\tag{3.2}\]

它的意义是:$b$到$Ax^$的距离最小. 其中,$Ax^$可以由矩阵乘法法则展开为:

\[Ax^* = \sum_{i=0}^na_ix^*_i\tag{3.3}\]

式子(3.3)代表着向量组$(a_1, a_2, a_3, …, a_n)$的任意线性组合,即向量空间$(a_1, a_2, a_3, …, a_n)$中的所有向量.

现在,我们需要知道,既然$Ax=b$是一个超定线性方程组,那么矩阵$A$的秩一定比增广矩阵$(A,b)$小,那么向量$b$一定不可以被列向量组$(a_1, a_2, a_3, …, a_n)$线性表出,这就意味着,在向量$b$中至少有一个分量在向量空间$(a_1, a_2, a_3, …, a_n)$之外,而最小化条件要求$b$到$Ax^$的距离最小,而$Ax^$可以是向量空间$(a_1, a_2, a_3, …, a_n)$中的所有向量,这就必然意味着向量$b - Ax^*$与所有的$a_i,i\in[1,n]$正交.

如果不理解的话可以在纸上画以下三个向量:

\[a_1 = (1,0,0)\] \[a_2 = (0,1,0)\] \[b = (1,1,1)\]

此处,与$b$距离最近且在向量空间$(a_1, a_2)$中的向量必然是$Ax^* = (1,1,0)$,此时$b - Ax^*$与所有的$a_i,i\in[1,n]$正交.

那么,我们就可以使用向量内积将问题转化为:

\[\left\{ \begin{aligned} &a_1^T(b-Ax^*)=0\\ &a_2^T(b-Ax^*)=0\\ &a_3^T(b-Ax^*)=0\\ &...\\ &a_n^T(b-Ax^*)=0\\ \end{aligned} \right.\]

上式可以转写为一个线性方程:

\[\begin{pmatrix} a_1^T\\ a_2^T\\ ...\\ a_n^T\\ \end{pmatrix} \begin{pmatrix} b-Ax^* \end{pmatrix} = \begin{pmatrix} 0\\ 0\\ ...\\ 0 \end{pmatrix}\]

而上式中左侧的矩阵恰好为$A^T$,即:

\[A^T(b-Ax^*)=0\] \[A^TAx^*=A^Tb\tag{3.4}\]

$\blacksquare$

注意,由于$A$不一定为方阵,即使为方阵也不一定满秩,故式(3.4)不可继续使用矩阵的标准逆元来化简,我们须使用其他的技术:

(以下内容复制自我之前的笔记,不代表目前水平)

矩阵的Moore–Penrose广义逆运算

https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse

定义

对于一个$m\times n$矩阵$A$,它的Moore–Penrose广义逆元被定义为:满足如下条件(Moore–Penrose条件)的一个$n\times m$矩阵$A^+$:

  1. $AA^+A=A$
  2. $A^+AA^+=A^+$
  3. $(AA^+)^T=AA^+$
  4. $(A^+A)^T=A^+A$

若矩阵$A$可逆,那么$A^{-1}$是其Moore–Penrose广义逆元.

我们仅需要证明其性质其一:

\[AB = ABB^+A^+AB\] \[AB = (AB)(AB)^+(AB)\]

得到

\[ABB^+A^+AB = AB(AB)^+AB\]

\[(AB)^+ = B^+A^+\]

通过SVD计算Moore–Penrose广义逆元

由SVD运算的拆分式:

\[M = U\Sigma V^T\]

若矩阵$A$可逆,那么$A^{-1}$是其Moore–Penrose广义逆元和上述证明的性质得到

\[M^+ = V\Sigma^+ U^T\]

而$\Sigma$是对角矩阵,其Moore–Penrose广义逆的算法明显是将其对角线元素取逆元(倒数).

由此,我们可以通过SVD计算出任意矩阵的Moore–Penrose广义逆元.

Moore–Penrose广义逆元用于超定方程组求解

我们回到超定方程组求解时的倒数第二步:

\[A^Tb-A^TAx^* = 0\]

此时,我们的$A$不是列满秩矩阵,因此不可用$A^TA$的一般逆来计算其值,因此我们应用Moore–Penrose广义逆运算(此时我们可以继续化简,因为即使是非方阵也具有Moore–Penrose广义逆元):

\[x^* = (A^TA)^+A^Tb = A^+A^{T+}A^Tb = A^+b\]

如此,我们解决了当属性矩阵列不满秩时的情况. 由于计算Moore–Penrose广义逆元需要更长的时间,所以我们还是建议在设计属性值的时候避免列不满秩的情况发生.

线性方程组迭代法,数值解

向量和矩阵的极限

记${ x^{(k)} }$是$\mathbb{R}^{n\times1}$上的向量序列,若:

\[\lim_{n\to\infin}\Vert x{(k)} - x \Vert = 0\]

我们就称向量序列${ x^{(k)} }$收敛至$x$.

向量收敛的充分必要条件是向量列表中的每一个元素$x_i^{(k)}$均收敛至$x_i$.

矩阵同理.

迭代法的一般形式

定理3.1:对恰定线性方程组$Ax=b$,构造与其同解的线性方程组

\[x = Mx+g\tag{3.6}\]

则$\forall x^{(0)}\in \mathbb{R}^{n\times 1}$,代入上述方程:

\[x^{(1)} = Mx^{(0)}+g\] \[x^{(k+1)} = Mx^{(k)}+g\]

若:当$k$充分大时$x^{(k)}$收敛至某一点$x^$,则$x^$是方程组$Ax=b$的解.

证明:

不妨设$x^{(0)}=b$,有$x^{(1)}=(E-A)b+b$,$x^{(2)}=(E-A)^2b+(E-A)b+b$,则

\[x^{(n)} = (E-A)^nb+(E-A)^{n-1}b+...+b\tag{3.6}\]

又,矩阵的乘法运算是结合的,而加法运算是交换的,故可以对式(3.6)应用等比数列求和公式:

\[x^{(n)} = \frac{b-(E-A)^{n+1}}{E-(E-A)}=\frac{b-(E-A)^{n+1}}{A}\tag{3.7}\]

上述分数线是方便表达的记法,指”分母的逆元左乘分子”.

已知上式收敛,故$\lim_{n\to\infin}x^{(n)} = A^{-1}b$

$\blacksquare$

常用迭代法

Jacobi迭代法

将方程组$Ax = b$展开:

\[\left\{ \begin{aligned} &a_{11}x_1 + a_{12}x_2 + a_{13}x_3 + ... + a_{1n}x_n=b_1\\ &a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + ... + a_{2n}x_n=b_2\\ &a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + ... + a_{3n}x_n=b_3\\ &...\\ &a_{n1}x_1 + a_{n2}x_2 + a_{n3}x_3 + ... + a_{nn}x_n=b_n\\ \end{aligned} \right.\]

移项:第$i$个方程左边只留下$x_i$:

\[\left\{ \begin{aligned} &x_1 = (b_1-a_{12}x_2 - a_{13}x_3 - ... + a_{1n}x_n)/a_{11}\\ &x_2 = (b_2-a_{21}x_1 - a_{23}x_3 - ... + a_{2n}x_n)/a_{22}\\ &x_3 = (b_3-a_{31}x_1 - a_{32}x_2 - ... + a_{3n}x_n)/a_{33}\\ &...\\ &x_n = (b_n-a_{n1}x_1 - a_{n2}x_3 - ... + a_{n,n-1}x_{n-1})/a_{nn}\\ \end{aligned} \right.\]

上式可以转写为$x = Mx+g$,其中

\[g = (b_1/a_{11}, b_2/a_{22}, ..., b_i/a_{ii}, ..., b_n/a_{nn})^T\] \[M = \begin{pmatrix} 0&-a_{12}/a_{11}& ... & -a{ij}/a_{ii} &... & -a{1n} / a_{11}\\ -a_{21}/a_{22}&0& ... & -a{ij}/a_{ii} &... & -a{2n} / a_{22}\\ -a_{31}/a_{33}&-a{32}/a_{33}& ... & -a{ij}/a_{ii} &... & -a{2n} / a_{22}\\ ...&...&-a{ij}/a_{ii}&...&...&...\\ -a_{n1}/a_{nn}&-a{n2}/a_{nn}& ... & -a{ij}/a_{ii} &... & 0\\ \end{pmatrix}\]

易知$x = Mx+g$与$Ax=b$同解.

将上述$M$与$g$重新表示:

定义$D=diag(a_{11}, a_{22}, …, a_{ii}, …, a_{nn})$

则$D^{-1} = diag(a^{-1}{11}, a^{-1}{22}, …, a^{-1}{ii}, …, a^{-1}{nn})$

则$M = E-D^{-1}A$,$g = D^{-1}b$

Gauss-Seidel迭代法

Jacobi迭代的运算过程可以描述如下:

\[\left\{ \begin{aligned} &x^{(k+1)}_1 = (b_1-a_{12}x^{(k)}_2 - a_{13}x^{(k)}_3 - ... + a_{1n}x^{(k)}_n)/a_{11}\\ &x^{(k+1)}_2 = (b_2-a_{21}x^{(k)}_1 - a_{23}x^{(k)}_3 - ... + a_{2n}x^{(k)}_n)/a_{22}\\ &x^{(k+1)}_3 = (b_3-a_{31}x^{(k)}_1 - a_{32}x^{(k)}_2 - ... + a_{3n}x^{(k)}_n)/a_{33}\\ &...\\ &x^{(k+1)}_n = (b_n-a_{n1}x^{(k)}_1 - a_{n2}x^{(k)}_3 - ... + a_{n,n-1}x^{(k)}_{n-1})/a_{nn}\\ \end{aligned} \right.\]

为节省存储空间和加快迭代进度,在运算过程中我们可以使用前步求出的结果代入后步:

\[\left\{ \begin{aligned} &x^{(k+1)}_1 = (b_1-a_{12}x^{(k)}_2 - a_{13}x^{(k)}_3 - ... + a_{1n}x^{(k)}_n)/a_{11}\\ &x^{(k+1)}_2 = (b_2-a_{21}x^{(k+1)}_1 - a_{23}x^{(k)}_3 - ... + a_{2n}x^{(k)}_n)/a_{22}\\ &x^{(k+1)}_3 = (b_3-a_{31}x^{(k+1)}_1 - a_{32}x^{(k+1)}_2 - ... + a_{3n}x^{(k)}_n)/a_{33}\\ &...\\ &x^{(k+1)}_n = (b_n-a_{n1}x^{(k+1)}_1 - a_{n2}x^{(k+1)}_3 - ... + a_{n,n-1}x^{(k+1)}_{n-1})/a_{nn}\\ \end{aligned} \right.\]

在实践过程中,只需要使用同一块存储单元,将求出的结果直接覆盖在旧的结果上即可.

数学上,我们将Jacobi迭代的$M$矩阵一分为二:

\[x^{(k+1)} = Lx^{(k+1)}+Ux^{(k)}+g\]

其中:$L$是下三角矩阵,$U$是上三角矩阵

\[L = \begin{pmatrix} 0&0& ... & 0 &... & 0\\ -a_{21}/a_{22}&0& ... &0 &... & 0\\ -a_{31}/a_{33}&-a{32}/a_{33}& ... & 0 &... & 0\\ ...&...&...&...&...&0\\ -a_{n1}/a_{nn}&-a{n2}/a_{nn}& ... & -a{ij}/a_{ii} &... & 0\\ \end{pmatrix}\] \[U = \begin{pmatrix} 0&-a_{12}/a_{11}& ... & -a{ij}/a_{ii} &... & -a{1n} / a_{11}\\ 0&0& ... & -a{ij}/a_{ii} &... & -a{2n} / a_{22}\\ 0&0& ... & -a{ij}/a_{ii} &... & -a{2n} / a_{22}\\ ...&...&...&...&...&...\\ 0&0& ... & 0 &... & 0\\ \end{pmatrix}\]

松弛迭代法

继续改进,将每一次迭代的Gauss-Seidel近似解$x^{(k+1)’}$与上一次迭代的$x^{(k)}$作加权平均:

\[x^{(k+1)} = \omega x^{(k+1)'} + (1-\omega )x^{(k)}\]

这样的方法称为松弛迭代.

其中$\omega$称为松弛系数,$\omega >1$称为高松弛,$\omega <1$称为低松弛,$\omega =1$为Gauss-Seidel法.

三种迭代方式的通式

观察得到,三种迭代方式的数学表达拥有共同的形式:

\[x^{(k+1)} = (L,U)(ax^{(k+1)}, bx^{(k)})^T\]

通过调整不同的$(a, b)$即可得到不同的迭代方式.

迭代法收敛的判断

谱半径

矩阵的谱指一个矩阵的特征值的集合. 矩阵的谱半径$\rho(A)$是指模最大的特征值.

由特征值的性质容易导出,$\rho(A^k)=\rho(A)^k$

关于谱半径有如下定理(证明见讲义):

  1. 矩阵谱半径不大于矩阵的任意范数
  2. 存在某矩阵范数$\Vert\Vert_i$,对于任意小的正数$\epsilon$,有:$\Vert A\Vert_i \leqslant \rho(A) + \epsilon$
  3. 若$A\in \mathbb{R}^{n\times n}$,则$\lim_{k\to \infin}A^k=0$当且仅当$\rho(A) < 1$

由上述定义及定理,立即推得:迭代式$x^{(k+1)} = Mx{(k)}+g$收敛,当且仅当$\rho(M) < 1$

有若干推论:

  1. 若$M$的任意范数$\Vert M\Vert < 1$,则收敛
  2. 松弛法收敛的必要条件是$\omega \in (0,2)$

快速判断收敛的若干结论

对角占优矩阵

若一个方阵某一行的对角元之绝对值大于该行所有其余元素绝对值之和,则称这个方阵是弱对角占优的. 若所有行都满足上述条件,则称这个方阵是严格对角占优的.

例:下述方阵就是严格对角占优的:

\[\begin{pmatrix} 100&1& 2 &3 \\ 1&90& 3 &6 \\ 0&1& 120 &-9 \\ 3&1& 2 & -203 \\ \end{pmatrix}\]

有若干结论:

  1. 若方阵是严格对角占优阵或者不可约的弱对角占优阵,则Jacobi和G-S必收敛
  2. 若方阵是严格对角占优阵且$\omega\in(0,1]$,松弛法收敛
  3. 若方阵是正定二次型且$\omega\in(0,2)$,松弛法收敛

误差估计

设有迭代式$x^{(k+1)} = Mx^{(k)}+g$收敛,则误差:

\[\Vert x^{(k)} - x^* \Vert \leqslant \frac{\Vert M \Vert^k}{1 - \Vert M \Vert}\Vert x^{(1)} - x^{(0)} \Vert\]

证明不难,见讲义.

借此,若定义了误差容许范围$\epsilon$,我们可以将其代入上式求出迭代次数的范围:

\[k \geqslant \frac{\ln\frac{\epsilon(1 - \Vert M \Vert)}{\Vert x^{(1)} - x^{(0)} \Vert}}{\ln \Vert M \Vert}\]

对于正定二次型的迭代方法

对于正定二次型有梯度下降法和共轭梯度法,本质上是非线性方程问题,不应在此处叙述.

2021.10.23 Hautbois


Copyright AmachiInori 2017-2021. All Right Reserved.
Powered By Jekyll.
amachi.com.cn