漫步线性代数七——特殊矩阵和应用

本篇文章有两个目标。第一是解释实际问题中大型线性方程组Ax=bAx=b的一种解法，事实是，工程或经济学中大型和现实的问题能够引导我们更深入理解这些知识。但是有一个很重要应用却不需要大量的准备工作。

另一个目标是说明系数矩阵具有的一些特殊性质，为了方便我们用同一个应用进行讲解。大型矩阵几乎总是有一个清晰的模式-对称和很多零元素。因为一个稀疏矩阵包含的信息远小于n2n^2个，所以计算应该更快。我们将观察带状矩阵，看看集中在对角线附近是如何加快消元的，为此我们将看到一个特殊的三对角矩阵。

看方程(6)中的矩阵，它是通过将微分方程变化为矩阵方程得到的。这是对每个xx求u(x)u(x)的连续问题，很明显计算机不能解决它，所以它必须近似为一个离散的问题-我们保留更多的未知变量，结果的精度就越好，当然计算代价也就越高。作为一个简单但仍然具有代表性的连续问题，我们选择微分方程

−d2udx2=f(x)0≤x≤1(1)

\begin{equation} -\frac{d^2u}{dx^2}=f(x)\qquad 0\leq x\leq 1\tag1 \end{equation}

这是关于位置函数u(x)u(x)的线性方程，可以在解中加上任何组合C+DxC+Dx依然满足要求，因为C+DxC+Dx的二阶导为零，不影响结果。对于两个任意常数C,DC,D的不确定性，通过在区间的两端添加一个边界条件就能够移除：

u(0)=0,u(1)=0(2)

\begin{equation} u(0)=0,\quad u(1)=0\tag2 \end{equation}这个结果是一个两点边值问题，描述的不是瞬变而是稳态现象-例如一根棒的温度分布，它的一端固定为 00C0^0C并且热源为 f(x)f(x)。

记住，我们的目标是产生一个离散的问题-换句话说，一个线性代数中的问题。为此我们只可以接受f(x)f(x)有限的信息，描述在nn个相等的区间点x=h,x−2h,…,x=nhx=h,x-2h,\ldots,x=nh上的值，对于同样位置处的真是解uu我们计算近似解u1,…,unu_1,\ldots,u_n，在端点处x=0,x=1=(n+1)hx=0,x=1=(n+1)h处，边界值是u0=0,un+1=0u_0=0,u_{n+1}=0。

第一个问题是：我们如何替换导数d2u/dx2d^2u/dx^2？一阶导数可以近似表示为有限步长内停止的Δu/Δx\Delta u/\Delta x并且不允许h(orΔx)h(or\Delta x)趋近于零，Δu\Delta u可以是前面的，后面的或中间的：

ΔuΔx=u(x+h)−u(x)h or u(x)−u(x−h)h or u(x+h)−u(x−h)2h(3)

\begin{equation} \frac{\Delta u}{\Delta x}=\frac{u(x+h)-u(x)}{h}\ or\ \frac{u(x)-u(x-h)}{h}\ or\ \frac{u(x+h)-u(x-h)}{2h}\tag3 \end{equation}

最后一个关于xx对称，它是最精确的。对于二阶导数，只是利用x,x±hx,x\pm h处数值的一个组合：

d2udx2≈Δ2uΔx2=u(x+h)−2u(x)+u(x−h)h2(4)

\begin{equation} \frac{d^2u}{dx^2}\approx\frac{\Delta ^2u}{\Delta x^2}=\frac{u(x+h)-2u(x)+u(x-h)}{h^2}\tag4 \end{equation}

它也有关于xx对称的优点。重复一遍：当h→∞h\to \infty时右边接近du/dx2d^u/dx^2 的真实值，但是我们必须让hh停在某个正数上。

对每个点x=jhx=jh，方程−d2u/dx2=f(x)-d^2u/dx^2=f(x)可以用它的离散模拟(5)代替，我们通过乘以h2h^2来得出nn个方程Au=bAu=b：

−uj+1+2uj−uj−1=h2f(jh)for j=1,…,n(5)

\begin{equation} -u_{j+1}+2u_{j}-u_{j-1}=h^2f(jh)\quad for\ j=1,\ldots,n\tag5 \end{equation}

第一个和最后一个(j=1,j=nj=1,j=n)包含u0=0,un+1=0u_0=0,u_{n+1}=0，他们是已知的边界条件，如果这些值非零的话，他们就转化成右边的值。这nn个方程(5)的结构可以更矩阵形式来更好的可视化，我们选择h=16h=\frac{1}{6}，从而得到5×55\times 5 的矩阵：

⎡⎣⎢⎢⎢⎢⎢⎢2−1−12−1−12−1−12−1−12⎤⎦⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢u1u2u3u4u5⎤⎦⎥⎥⎥⎥⎥⎥=h2⎡⎣⎢⎢⎢⎢⎢⎢f(h)f(2h)f(3h)f(4h)f(5h)⎤⎦⎥⎥⎥⎥⎥⎥(6)

\begin{equation} \begin{bmatrix} 2&-1&&&\\-1&2&-1&&\\&-1&2&-1&\\&&-1&2&-1\\&&&-1&2 \end{bmatrix} \begin{bmatrix} u_1\\u_2\\u_3\\u_4\\u_5 \end{bmatrix} =h^2\begin{bmatrix} f(h)\\f(2h)\\f(3h)\\f(4h)\\f(5h) \end{bmatrix}\tag6 \end{equation}

现在，我们将求解方程(6)，它的系数矩阵非常有规律，有许多特殊的性质，其中有三个是非常基本的：

矩阵AA是三对角的。所有非零元素位于主对角线以及附近的两条对角线上，这条窄带以外的aij=0a_{ij}=0，这些零大大简化了高斯消元过程。
矩阵是对称的。每个元素aija_{ij}等于它的镜像ajia_{ji}，使得AT=AA^{T}=A。上三角矩阵UU将是下三角矩阵LL的转置，A=LDLTA=LDL^{T}。AA的对称性反映了d2u/dx2d^2u/dx^2的对称性，奇导数像du/dx,d3u/dx3du/dx,d^3u/dx^3将破坏对称性。
矩阵是正定的。这个额外的性质说明主元是正的，在理论和实践中不需要行变换。这和本文末尾要将的矩阵BB(非正定的)正好相反，在没有行变换的情况下，它将对舍入非常敏感。另外关于正定的概念我会在以后的文章中详细介绍！

我们返回到AA是三对角矩阵这个事实，它对消元会有什么影响呢？消元过程的第一步是在第一个主元下面产生零：

⎡⎣⎢⎢⎢⎢⎢⎢2−1−12−1−12−1−12−1−12⎤⎦⎥⎥⎥⎥⎥⎥→⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢20−132−1−12−1−12−1−12⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥

\begin{bmatrix} 2&-1&&&\\-1&2&-1&&\\&-1&2&-1&\\&&-1&2&-1\\&&&-1&2 \end{bmatrix} \to \begin{bmatrix} 2&-1&&&\\0&\frac{3}{2}&-1&&\\&-1&2&-1&\\&&-1&2&-1\\&&&-1&2 \end{bmatrix}

跟一般的5×55\times 5矩阵相比，这一步主要有两个简化：

在主元下面只有一个非零元素
主元所在的行非常短

乘数因子是ℓ=−12\ell=-\frac{1}{2}，新的主元是32\frac{3}{2}。更进一步，三对角模式保留着：每步消元都允许这两个简化。

最终结果是LDU=LDLTLDU=LDL^{T}，注意主元！

A=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢1−121−231−341−451⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢2132435465⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢1−121−231−341−451⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

A= \begin{bmatrix} 1&&&&\\-\frac{1}{2}&1&&&\\&-\frac{2}{3}&1&&\\&&-\frac{3}{4}&1&\\&&&-\frac{4}{5}&1 \end{bmatrix} \begin{bmatrix} \frac{2}{1}&&&&\\&\frac{3}{2}&&&\\&&\frac{4}{3}&&\\&&&\frac{5}{4}&\\&&&&\frac{6}{5} \end{bmatrix} \begin{bmatrix} 1&-\frac{1}{2}&&&\\&1&-\frac{2}{3}&&\\&&1&-\frac{3}{4}&\\&&&1&-\frac{4}{5}\\&&&&1 \end{bmatrix}

三对角矩阵的L,UL,U分解因子是二对角矩阵，三个因子和矩阵AA一样对角线有同样的带状结构，注意L,UL,U互相之间存在转置关系，我们从对称就预期到会如此，主元2/1,3/2,4/3,5/4,6/52/1,3/2,4/3,5/4,6/5都是正的，他们的乘积就是AA的行列式detA=6detA=6。当nn不断变大时，这些主元明显收敛到1，这样的矩阵计算起来非常方便。

稀疏因子L,UL,U完全改变了通常的操作次数，每列的消元只需要两步，对于nn列来说，次数从n3/3n^3/3降到了2n2n，三对角方程组几乎能很快解决，求解它的代价和nn成正比。

带状矩阵就是在|i−j|<w|i-j|内aij=0a_{ij}=0(图1)，对角矩阵的半个带宽w=1w=1，三对角矩阵的为w=2w=2,，当w=nw=n时就是就是一般的矩阵了。对每一列，消元法需要w(w−1)w(w-1)次操作：长为ww的行对下面的w−1w-1行进行操作，对于nn列的带状矩阵大约需要w2nw^2n次操作。

当ww趋近nn时，矩阵变成一般矩阵，次数变成n3n^3。产生L,D,UL,D,U的除法和乘-减(不考虑AA为对称这个假设)精确次数是P=13w(w−1)(3n−2w+1)P=\frac{1}{3}w(w-1)(3n-2w+1)，对于一般矩阵即w=nw=n，P=13n(n−1)(n+1)P=\frac{1}{3}n(n-1)(n+1)，这是一个整数，因为n−1,n,n+1n-1,n,n+1是连贯的，所以他们之中有一个能被3整除。

图1：带状矩阵AA和它的因子L,UL,U

这是我们最终得到的操作次数，我们强调一点，像AA那样的有限差分矩阵有逆，在求解Ax=bAx=b时，知道A−1A^{-1}比知道L,UL,U要糟糕，因为A−1A^{-1}乘bb需要n2n^2步，而前向消元和回代(产生x=U−1c=U−1L−1b=A−1bx=U^{-1}c=U^{-1}L^{-1}b=A^{-1}b)4n4n 步就足够了。

我们希望本例增强读者对消元的理解，它是实践中遇到的一个大型线性方程实例，接下来对于nn个未知量的mm个方程，我们将转向讨论xx的存在性和唯一性。

舍入误差

理论上非奇异情况有完整的主元(考虑行变换)，实践中，需要更多的行交换否则计算的解可能变得毫无价值。接下来我们重点来讲如何使消元更加稳定-为什么需要它以及如何做。

对于中等大小的系统，比如说100×100100\times 100，消元可能涉及一百万的三分之一次(13n3\frac{1}{3}n^3)，对于每步操作，我们必须预期舍入误差。通常情况下，我们固定有效位数，然后将两个大小不同的数相加给出一个误差：

.456+.00123→.457丢掉数字2和3

.456+.00123\to .457\quad \text{丢掉数字2和3}

那么所有这些误差是如何影响Ax=bAx=b的最终误差的呢？

这不是一个简单的问题。而这个问题被约翰⋅\cdot冯⋅\cdot 诺依曼碰到了，在计算机使得100万步操作成为可能的时候，是他引领了数学。事实上高斯和冯⋅\cdot诺依曼的结合给出了简单的消元算法，虽然冯诺依曼高估了最后的舍入误差。威尔金森(Wilkinson)找到这个问题的正确方法，并且他的书到现在依然是经典。

我们将给出两个简单的例子来说明舍入误差的重要性：

A=[1.1.1.1.0001]B=[.0011.1.1.]

A= \begin{bmatrix} 1.&1.\\1.&1.0001 \end{bmatrix} \qquad B= \begin{bmatrix} .001&1.\\1.&1. \end{bmatrix}

AA几乎是奇异的而BB离奇异很远，如果我们稍微将AA中的元素改一下a22=1a_{22}=1，它既是奇异的。考虑两个非常接近的向量bb：

ill−conditioned uu++v1.0001v==22andwell−conditioned uu++v1.0001v==22.0001

ill-conditioned\ \begin{matrix} u&+&v&=&2\\u&+&1.0001v&=&2 \end{matrix} \quad and\quad well-conditioned\ \begin{matrix} u&+&v&=&2\\u&+&1.0001v&=&2.0001 \end{matrix}

第一个解是u=2,v=0u=2,v=0，而第二个解是u=v=1u=v=1。bb中第五位数的改变放大到解中第一位数的改变，没有任何数值方法可以避免它对小扰动这种敏感，病态条件能够从一个地方转到另一个地方，但是无法移除。真实解都如此敏感更别说计算解了。

第二点如下：

15、即便是一个well−conditionedwell-conditioned矩阵BB也可能别差的算法个毁掉

我们很遗憾的说：对于矩阵BB，直接使用高斯消元就是一种差的算法。假设.0001作为第一个主元，那么第二行需要减去第一行的10000倍，右下角就变成了-9999，但是第三位四舍五入将变为-10000，为1的元素将会消失：

.0001u+vu+v==22→.0001u+v−9999v==1−9998

\begin{matrix} .0001u+v&=&2\\u+v&=&2 \end{matrix} \to \begin{matrix} .0001u+v&=&1\\-9999v&=&-9998 \end{matrix}

四舍五入将得到−10000v=−10000-10000v=-10000或者v=1v=1。首先我们用正确值v=.999v=.999(保留三位小数)回代将得到u=1u=1：

.0001u+.9999=1oru=1

.0001u+.9999=1\quad or\quad u=1

但是如果接受v=1v=1，我们将得到u=0u=0：

.0001u+1=1oru=0

.0001u+1=1\quad or\quad u=0

计算结果完全不对，BB是well-conditioned但是消元明显不稳定，L,D,UL,D,U也明显是比BB大许多：

B=[11000001][.000100−9999][10100001]

B= \begin{bmatrix} 1&0\\10000&1 \end{bmatrix} \begin{bmatrix} .0001&0\\0&-9999 \end{bmatrix} \begin{bmatrix} 1&10000\\0&1 \end{bmatrix}

小主元.0001带来了不稳定，补救方法很明显-换行

16、一个小主元迫使消元中发生变化。通常我们比较同一列所有可能的主元，将最大主元换到当前的位置

对于BB来讲，.0001将和下面的1比较，立马进行行交换，从矩阵的角度来讲就是乘以一个置换矩阵PP，新矩阵C=PBC=PB就有好的因子：

C=[1.000111]=[1.000101][100.9999][1011]P=[0110]

C=\begin{bmatrix} 1&1\\.0001&1 \end{bmatrix}= \begin{bmatrix} 1&0\\.0001&1 \end{bmatrix} \begin{bmatrix} 1&0\\0&.9999 \end{bmatrix} \begin{bmatrix} 1&1\\0&1 \end{bmatrix} \qquad P= \begin{bmatrix} 0&1\\1&0 \end{bmatrix}

CC的主元是1和.9999，比BB的.0001和-9999要好。

还有一种策略就是与其余所有列中最大主元进行交换，这时候可能不仅是行，也会有列交换。(这时候置换矩阵乘在右边)这么做的代价太高，上面的方法其实已经够了。

我们说完了数值线性代数的基本算法：带有行变换的消元法。一些进一步的完善，比如看整行或列，都是有可能的，但本质上读者现在知道了一台电脑如何解线性方程组，相比“理论”描述-找到A−1A^{-1}并相乘A−1bA^{-1}b-我们的描述已花费了大量时间和耐心，我希望能够用更简单的方法来解释x<script type="math/tex" id="MathJax-Element-7260">x</script>是如何发现的，但是我认为我还没有找到。(博主心声：期待大家能有更好的方法提出来)

漫步线性代数七——特殊矩阵和应用相关推荐

漫步线性代数二十七——矩阵对角化
现在我们开始实质性的计算,它非常简单并且在随后的几篇文章里都会用到.特征向量对角化一个矩阵: 3.假设n×nn\times n矩阵有nn个线性无关的特征向量,如果这些向量是矩阵SS的列,那么S−1AS ...
漫步线性代数二——线性方程的几何形状
漫步线性代数二--线性方程的几何形状 2016年08月15日 23:10:10 会敲键盘的猩猩阅读数:1818 几何形状理解这个主题的方法是举例说明.我们以两个极其简单的方程开始,可以说大家在没有 ...
三阶矩阵的lu分解详细步骤_数学 - 线性代数导论 - #4 矩阵分解之LU分解的意义、步骤和成立条件...
线性代数导论 - #4 矩阵分解之LU分解的意义.步骤和成立条件目前我们用于解线性方程组的方法依然是Gauss消元法.在Gauss消元法中,我们将右侧向量b与A写在一起作为一个增广矩阵进行同步的操作 ...
线性代数之行列式矩阵术语中英对照
线性代数之行列式矩阵术语中英对照本文基于Steven J. Leon著的第九版<Linear Algebra with Applications>整理,后续会输出其它章节,敬请期待. d ...
线性代数导论2——矩阵消元
线性代数导论2--矩阵消元本文是Gilbert Strang的线性代数导论课程笔记.课程地址:http://v.163.com/special/opencourse/daishu.html 第二课时 ...
线性代数分块矩阵求逆矩阵_单位矩阵属性（AI = A）| 使用Python的线性代数
线性代数分块矩阵求逆矩阵 Prerequisites: 先决条件: Defining Matrix 定义矩阵 Identity matrix 身份矩阵 numpy.matmul( ) matrix m ...
干货来袭！！！3天0基础Python实战项目快速学会人工智能必学数学基础全套（含源码）（第1天）线性代数篇：矩阵、向量及python实战
第1天:线性代数篇:矩阵.向量.实战编程第2天:微积分篇:极限与导数.梯度下降.积分.实战编程第3天:概率分析篇:条件概率与全概率.贝叶斯公式.实战项目目录前言一.矩阵在AI中的应用二.矩 ...
MIT线性代数笔记三矩阵的乘法和逆矩阵
文章目录 1. 矩阵乘法 Matrix multiplication 1.1 标准方法(行乘以列) 1.2 列向量的线性组合 1.3 行向量的线性组合 1.4 分块乘法 2. 逆矩阵 2.1 逆矩阵的 ...
MIT线性代数1806(8) 矩阵秩特解通解
MIT线性代数1806(8) 矩阵秩特解通解

漫步线性代数七——特殊矩阵和应用

舍入误差

漫步线性代数七——特殊矩阵和应用相关推荐

最新文章

热门文章