线性代数数值计算注解

线性代数数值计算注解

Chp1 线性代数中的方程组

1. 浮点数，因为储存原因存在“四舍五入”误差
2. 对于方阵，一般使用消元法
3. 部分主元法：通常选择一列中绝对值最大的元素作为主元，减少舍入误差
4. 矩阵先前化简需要$O(n^3)$计算，向后化简只需要次$O(n^2)$次计算
5. C语言编写计算$AX$ 会考虑使用$X^T{A}^T$

Chp2 矩阵代数

1. 计算$AB$的速度依赖于储存方法，按列计算，可以使用并行算法，或按行计算

2. 实际中很少计算$A^{-1}$，除非需要其中的元素，计算$A^{-1}$次数为解$Ax=b$的3倍，且行变换可能更准确

3. 实际中会遇到“接近奇异”病态的矩阵，该矩阵可逆，但是稍微改变其中的元素就会变成奇异矩阵，此时行变换可能由于舍入误差产生少于n个主元位置，有时舍入误差也可能奇异矩阵变成可逆

4. 某些程序计算条件数，条件数越大，越接近奇异。单位阵条件数为1，奇异矩阵条件数无穷大，极端情况下，矩阵程序可能无法区分奇异矩阵和病态矩阵

5. 当方程$Ax=b$的右端有一点变化，例如$Ax=b+\Delta b$，方程的解变为$x+\Delta x$，此处$\Delta x$满足$A\Delta x=\Delta b$，商$||\Delta b||/||b||$称为$b$的相对改变(误差)，解的相对改变为$||\Delta x||/||x||$，当A可逆时，记A的条件数为$cond(A)$，不可逆时A的条件数为无穷大，有如下不等式
   $$\frac{||\Delta x||}{||x||}\le cond(A)\cdot \frac{||\Delta b||}{||b||}$$

6. 例：P116， Hilbert矩阵$A_{ij}=\frac{1}{i+j-1}$，
   $$\left[\begin{array}{cccc}1& 1/2& 1/3& 1/4\\ 1/2& 1/3& 1/4& 1/5\\ 1/3& 1/4& 1/5& 1/6\\ 1/4& 1/5& 1/6& 1/7\end{array}\right]$$

7. MatLab命令:hilb()，高度病态。11阶或以上Hilbert矩阵MatLab不能正常求逆

8. MatLab Linsolve使用LU分解，或QR分解，解线性方程组,但是位于$m\times n$的矩阵只返回一个特解$m>n$

9. 矩阵过大的时候，使用分块矩阵会更有效

10. 对于解方程$Ax=b$，可以把A写成$I-C$，$x=(I-C{)}^{-1}b$，其中$(I-C{)}^{-1}$可以如下计算
    $$(I-C)(I+C+C^2+C^3+\cdots +C^m)=I-C^{m+1}$$
    当$m\to \infty$时，若有$C^m\to 0$，则可以使用$I+C+C^2+C^3+\cdots +C^m$作为$(I-C{)}^{-1}$近似。

矩阵LU分解

对于$n\times n$的稠密矩阵，$n$足够大$n\ge 30$，有

1. 计算$A$的LU分解大约需要$2n^3/3$浮点运，近似于化简$[Ab]$相同,而求$A^{-1}$需要$2n^3$次运算
2. 解 $Ly=b$ 和 $Ux=y$，即 $LUx=b$ 大约需要$2n^2$次浮点运算，任意$n\times n$三角方程组可以大约$n^2$用浮点算解出
3. 把b乘以$A^{-1}$也需要$2n^2$次浮算，而且有可能不如LU分解计算来的精确，由于计算$A^{-1}$和$A^{-1}b$存在舍入误差
4. 若$A$是稀疏矩阵，如带状，有两到三条对角线，则L和U也可能是稀疏的然而$A^{-1}$可能是稠密的，此时LU分解解方程$AX=b$可能会快很多

Chp3 行列式

1. 行列式的计算不使用代数余子式展开，而是用消元法化简成上三角形式，$n\times n$行列式大约需要$2n^3/3$次计算
2. 克拉默法则仅限于理论分析稳定性，并不实际应用，对于逆矩阵公式夜市同样情况

Chp4 向量空间

1. 确定矩阵的秩通常使用SVD分解，考虑$\Sigma$对角矩阵中非零元个数，使用行化简形式来确定秩会因为舍入误差而导致得出错误的结果，例如
   $$\left[\begin{array}{cc}5& 7\\ 5& 6.99999999\end{array}\right]$$
   这取决于计算机储存及计算精度。SVD分解还用于得到四个子空间的基向量。

2. MatLab使用rref来确定NullSpace，null命令

3. 计算$P^k{x}_0$，先计算$P^k$再计算$P^k{x}_0$和计算$P{x}_0,P\cdot P{x}_0,P\cdot P^2{x}_0$不同，前者需要的计算更多。

Chp5 特征值和特征向量

1. 对于一般的矩阵($n\ge 5$)，没有公式或有限算法求解特征方程.

2. MatLab先求特征值${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n$，然后展开$(\lambda -{\lambda }_1)(\lambda -{\lambda }_2)\cdots (\lambda -{\lambda }_n)$的积来的得到特征多项式

3. 有相同的特征值不一定相似

4. 有些算法基于矩阵相似估计特征值，QR算法能有效估计特征值。当$A=A^T$时，可以使用雅可比方法来计算形如
   $$A_1=A,A_{k+1}=P_k^{-1}{A}_k{P}_k(k=1,2,\cdots )$$
   的矩阵序列。序列中每个矩阵都相似与$A，因此与$$A$有相同的特征值，当k增大时，$A_{k+1}$的非对角线元素趋于0，而主对角线上的元素就近似为A的特征值。

5. 幂算法迭代估计绝对值最大特征根Page317

6. 逆幂法估计A的特征值和特征向量Page320

7. 幂算法和逆幂法对于简单情况适用

QR算法求特征值

• 算法产生矩阵序列，序列中的矩阵全部相似与A，矩阵几乎是上三角的，并且主对角线上的元素近似为$A$的特征值
• 主要思想：做QR分解。把A或相似于A的矩阵分解为$A=Q_1{R}_1$,这里$Q_1^T=Q_1^{-1}$，Q为正交矩阵，而$R_1$是上三角矩阵。交换$Q_1$与$R_1$形成$A=R_1{Q}_1$，$A_1$又被分解成为$A_1=Q_2{R}_2$，然后令$A_2=R_2{Q}_2$，依次构造$A_k$
• MatLab使用QR算法计算特征值和特征向量，参见一般数值分析教材。对一元高次方程，MatLab将其转化为特征方程形式求特征值，特征值即为方程的根。(roots函数)
• 一个豪斯霍尔德矩阵或者基本矩阵镜像具有如下形式$Q=I-2u^{T}u$，此处u为单位向量，则Q为正交矩阵。基本镜像经常用于在计算机程序中产生矩阵A的一QR分解。如果A具有线性无关的列，那么一系列基本镜像的左乘可以产一个上三角矩阵

Chp6 正交性和最小二乘法

1. QR 分解得到$R$为可逆矩阵(坐标唯一性)，

2. 在使用格拉姆-施密特方法的过程中，因为每个向量都会四舍五入，引起后面得到的向量误差逐渐增大，对于较大的$j和k$($j̸=k$)，$u_j^T{u}_k$也许不会充分接近0，通过重新安排计算的阶，这类正交性的损失可以大大减少。矩阵的QR分解常常使用格拉姆-施密特正交化方法，因为它会得到更为精确的标准正交基，但是分解时间有可能加倍

3. 为了实现矩阵$A$的QR分解，程序常常会对$A$左乘一些列正交矩阵使得结果变成一个上三角矩阵，这个构造过程有点类似于$A$左乘一系列初等矩阵，最后得到$A$的LU分解。

4. 对矩阵$A$做QR分解时，使用格拉姆-施密特方法得到正交矩阵Q，那么$R=Q^{-1}A=Q^{T}A$

5. $A$的QR分解，$A=[x_1,x_2,\cdots ,x_p]$，若$n\times k$矩阵Q的列构成A的前k列子空间$W_k$的一个标准正交基，那么对于${\mathbb{R}}^n$中的向量$x$，$Q{Q}^{T}x$是$x$在$W_k$上的正交投影。如果$X_{k+1}$是A的下一列，那么可以取
   $$v_{k+1}=x_{k+1}-Q(Q^T{x}_{k+1})$$
   括号用于减少计算，取$u_{k+1}=v_{k+1}/\mid \mid v_{k+1}\mid \mid$，得到新的Q，重复上述步骤即可的得到$A$的QR分解。

6. 某些时候，最小二乘问题的法方程(正规方程)可能是病态的(A的两列高度相关)，也就是$A^{T}A$中的元素在计算中出现较小的误差，可能导致$\hat{x}$的较大误差。如果A的列线性无关，最小二乘解常常可通过A的QR分解可靠的求出。

7. 最小二乘求解
   $$\begin{array}{rl}A& =QR\\ A^{T}A\hat{x}& =A^{T}b\\ A^{T}A& =R^T{Q}^{T}QR=R^{T}R\\ R^{T}R\hat{x}& =R^T{Q}^{T}b\\ R\hat{x}& =Q^{T}b\end{array}$$
   解最后一个方程，使用行变换会更快。

8. $rankA=rank{A}^{T}A$

Chp7 对称矩阵和二次型

谱定理

一个对称的$n\times n$实矩阵具有下面性质

1. A有n个实特征值，包括重复的特征根
2. 特征空间相互正交，这种正交性是在特征向量对应不同特征值意义下成立的
3. A可以对角化
4. 对每一个特征值，对应特征空间的维数等于$\lambda$作为特征方程的更的重数

证明思路

性质1

由代数基本定理，任意一元n次方程有n个复根，设矩阵A的特征方程根为$\lambda$，对应的特征向量为$x$

$$\begin{array}{rl}Ax& =\lambda x\\ A\bar{x}& =\bar{\lambda }\bar{x}\\ {\bar{x}}^{T}A& =\bar{\lambda }{\bar{x}}^T\\ {\bar{x}}^{T}Ax& =\bar{\lambda }{\bar{x}}^{T}x\end{array}\text{\hspace{1em}(1)}$$
又因为
$$\begin{array}{rl}Ax& =\lambda x\\ {\bar{x}}^{T}Ax& =\lambda {\bar{x}}^{T}x\end{array}\text{\hspace{1em}(2)}$$
所以由(1)(2)可得
$$\bar{\lambda }{\bar{x}}^{T}x=\lambda {\bar{x}}^{T}x$$
所以
$$\bar{\lambda }=\lambda$$
所以$\lambda $是实数

性质2

设$v_1,v_2$是对应不同特征值${\lambda }_1,{\lambda }_2$的特征向量，那么
$$\begin{array}{rl}{\lambda }_1{v}_1^T{v}_2& =(A{v}_1{)}^T{v}_2\\ & =v_1^{T}A{v}_2\\ & ={\lambda }_2{v}_1^T{v}_2\end{array}$$
所以$({\lambda }_1-{\lambda }_2)v_1^T{v}_2=0$，所以$v_1^T{v}_2=0$

性质3

首先引入矩阵舒尔分解(实舒尔分解，$A$,$U$,$R$均为实矩阵)。

舒尔分解：$n\times n$矩阵A可以分解成$A=UR{U}^T$，且$U$为正交矩阵，$R$为上三角矩阵。

1. 如果矩阵$A$可以舒尔分解，那么$A$具有n个特征值，计算包含重数，证明简单从略
2. 假设$A$具有n个实特征根(${\lambda }_1,{\lambda }_2，\cdots ,{\lambda }_n$)，则A有一个实舒尔分解

设$u_1$是对应于${\lambda }_1$的单位特征向量，$u_2,\cdots ,u_n$是其余向量，且 { u 1 , u 2 , ⋯   , u n } \{u_1,u_2,\cdots,u_n\} {u1​,u2​,⋯,un​}是${\mathbb{R}}^n$的单位正交基，取$U=[u_1,u_2,\cdots ,u_n]$，那么$U^{T}AU$的第一列为
$$U^T{\lambda }_1{u}_1={\lambda }_1\left[\begin{array}{c}{u}_1^T\\ u_2^T\\ ⋮\\ u_n^T\end{array}\right]u_1={\lambda }_1\left[\begin{array}{c}1\\ 0\\ ⋮\\ 0\end{array}\right]$$
这表明$U^{T}AU$具有以下形式
$$\left[\begin{array}{ccccc}{\lambda }_1& \ast & \ast & \ast & \ast \\ 0& & & & \\ ⋮& & A_1& & \\ 0& & & & \end{array}\right]$$
其中$A_1$的特征向量为${\lambda }_2,\cdots ,{\lambda }_n$。

由舒尔分解性质2知，对称矩阵有舒尔分解$A=UR{U}^T$，又因为A是对称的，两边转置得到
$$A^T=U{R}^T{U}^T$$
所以$R^T=R$，又因为R是上三角矩阵，R必须是对角矩阵。于是矩阵A的舒尔分解$UR{U}^T$就是A的正交对角化。

性质4

由性质3，容易知道对于k重根$\lambda$，其必定对应k个线性无关特征向量。
$$\begin{array}{rl}A& =QD{Q}^T=\left[\begin{array}{c}{u}_1,\cdots ,u_n\end{array}\right]\left[\begin{array}{cccc}{\lambda }_1& & 0& \\ & \ddots & & \\ 0& & {\lambda }_n& \end{array}\right]\left[\begin{array}{c}{u}_1^T\\ ⋮\\ u_n^T\end{array}\right]\\ & =\left[\begin{array}{c}{\lambda }_1{u}_1,\cdots ,{\lambda }_n{u}_n\end{array}\right]\left[\begin{array}{c}{u}_1^T\\ ⋮\\ u_n^T\end{array}\right]\\ & ={\lambda }_1{u}_1{u}_1^T+\cdots +{\lambda }_n{u}_n{u}_n^T\end{array}$$

数值计算注解

• 对于对称矩阵，计算机可以快速精确计算特征值和特征向量，见Chp5对角化注解，利用正交矩阵可逆避免误差累计
• 确定对称矩阵A是正定的最快方式：A有n个正主元，类似于LU分解，通过对A消元即可。当A是对称矩阵时，$A=R^{T}R$，其中R是上三角矩阵。$A=QD{Q}^T$，令$D=C^{T}C$，C也是对角矩阵，$A=QC^TCQT = $，令$R=QCQ^T$，则$RTR = QC^TQTQCQ^T = QDQ^T$。
• 估计大矩阵的秩时，最可靠的方法是计算非零奇异值的个数(MatLab使用该方法[^1])，在这种情况下，特别小的非零奇异值在实际计算中常假定为零，矩阵A的有效秩是剩余非零奇异值数目。
• 奇异值计算应该避免计算$A^{T}A$，原因是任何A中元素的误差在$A^{T}A$中被平方，存在快速的迭代方法，可计算到精确很多位数的矩阵A的奇异值和左右奇异向量。

[1] MatLab 2018a rank 函数源代码

Refference :Linear Algebra and It’s Applications 4th, David C. Lay, China Machine Press