MATLAB 学习指南

MATLAB 学习指南
第四章. 高级矩阵运算
4.1. 稀疏矩阵
稀疏矩阵
为了展现使用稀疏矩阵带来的有效性，我们将两次使用有限差分法来解一个PDE。

首先，我们将会使用矩阵指令。

矩阵指令会使用到目前为止我们已经学过的全矩阵。

其次，我们将使用新的指令。

新的指令充分利用大多数元素都是零这一事实来极大地减少内存需求量和解PDE所需的浮点运算的次数。

使用全矩阵格式进行计算
通过使用中心有限差分法来离散1维空间中的拉普拉斯算子可得到下列矩阵。

设定离散PDE所得的矩阵，使其1维方格包含点间距为1的num_pts点。

在x=1和x=num_pts处的Dirichlet边界条件已经被设定好了。

对于内部点，我们有一个源项。

现在，我们使用标准MATLAB解算程序来获得格点处PDE的解。

现在我们进一步观察Afull。

指令spy(A)能够绘制矩阵A的图形，具体做法是：无论何处，只要矩阵A的一个元素非零，就编写一个点。

底部的数字nz是非零元素的个数。

我们可以看到只有一小部分矩阵元素非零。

由于我们按常规的方式为格点编号，每个格点的邻近格点都存储在邻近区域，所以此矩阵中的非零元素位于主对角线上以及紧邻两条对角线的上方和下方。

即使我们不按常规给格点编号，我们得到的非零点仍然是这么少。

这种情况是我们经常遇到的，即，在PDE的数值模拟中，我们常常遇到的矩阵都是稀疏的，那就是说，只有一小部分点是非零的。

对于这个矩阵，元素的总数目是：
这就意味着Afull的大部分空间是空的，我们浪费了很多内存来存储我们已知的值为零的元素。

除了Afull，清除内存中的所有变量。

稀疏矩阵
使用指令“sparse”，我们可将一个矩阵转化成为稀疏的格式。

MATLAB把一个稀疏矩阵当作一个NZ×3型的数组来存储，其中，NZ是非零元素的个数。

第一列是行数，第二列是非零元素的列数，第三列是非零元素的实际值。

总内存使用量远小于用全矩阵格式时的总内存使用量。

现在，我们将使用稀疏矩阵格式来解决问题
接下来，我们设定格点值。

现在我们声明：从一开始矩阵A就有稀疏矩阵的结构。

首先，我们计算非零元素的个数（或此数的上限）。

我们可以看到，对于相应于每个内部点的每一行，我们都有3个值，但是，对于第一行和最后一行，我们仅有一个值。

因此，非零元素的个数是：
我们现在使用“spalloc(m,n,nz)”指令来为m×n型的，仅有nz个非零元素的稀疏矩阵分配内存。

我们现在设定A矩阵的值。

在x=1和x=num_pts处的Dirichlet边界条件已经被设定好了。

对于内部点，我们有一个源项。

现在，当我们调用MATLAB标准解算程序时，它能自动辨别出A是一个稀疏矩阵并利用这一事实使用解算程序的运算法则。

从A和Afull的几行代码中可以看出稀疏矩阵格式所需的内存量远小于全矩阵格式所需的内存量。

同样，如果N是格点数，我们看到全矩阵的大小是N^2；但是，稀疏矩阵所需的内存大小仅约为3*N。

因此，随着N的增加，稀疏矩阵格式将比全矩阵格式有效得多。

对于具有数千格点的复杂模拟来说，我们只能寄希望于利用稀疏性来解决这些问题。

使用稀疏矩阵还可以提高执行速度，为看到执行速度的提高，仔细检查下面两种矩阵乘法的算法。

矩阵乘法的全矩阵算法
矩阵乘法的稀疏矩阵算法
最后，我们注意到求一个稀疏矩阵的逆阵通常会极大地破坏其稀疏性。

因此，如果我们有了矩阵A和矩阵C=A*B的值，想要计算矩阵B，我们不要使用inv(A)*C，而要使用“left matrix division”算符A\C。

它会返回一个与inv(A)*C等价的矩阵，但使用MATLAB解算程序可以充分利用稀疏性。

我们看到消去法产生的误差已经把非常小的非零值插入到远离中心三对角线的元素中了。

我们要想去除这些元素只能保留大于容限值的元素。

由于我们不想经历中间步骤，因为在中间步骤我们不得不存储具有许多非零元素的矩阵，所以通常我们不以这种方式来计算矩阵。

相反，我们仅限于解A*x=b这种形式的线性系统，其中，x和b是向量，A是稀疏矩阵，我们可以直接输入矩阵A的值。

因此，我们就避免了与从舍去误差中产生许多非零元素这一问题有关的内存问题。

4.2. 普通矩阵运算/本征值
使用“det”可用来计算方阵的行列式。

矩阵的其它普通函数：
使用指令“eig”可计算矩阵的本征值。

如果也要求计算本征向量，句法是：
这里V是包含本征向量（作为列向量）的矩阵，D是包含本征值的对角阵。

使用指令“eig(A,k)”可以计算A的k首项本征值；那就是说，计算具有最大模量的k本征值。

相似地，首项本征值的本征向量也可以用eigs来计算。

就稀疏矩阵来说，只有指令“eigs”可以被使用。

4.3. LU分解
使用多个b向量并使用LU分解，可解决线性系统Ax=b。

这里，我们来分解P*A=L*U，其中，P是一个置换矩阵（因此inv(P)=P’），L是一个下三角阵，U是一个上三角阵。

如果在因式分解的过程中没有进行枢轴选择（从本质上说，这是高斯消去），P就是一个恒等矩阵。

一旦完成LU因式分解，使用下列线性代数步骤就可解决Ax=b的问题。

这就得到了下列两个线性问题，包括可以用替换法解的三角矩阵。

进行LU因式分解的MATLAB指令是“lu”。

我们使用一个随机的非奇异阵来证明算法。

通过增加一个恒等矩阵来确保非奇异性。

进行LU因式分解。

比较所包含的矩阵的结构
调用LU因式分解和调用稀疏矩阵的方法相同；但是，一般而言，因式化的矩
阵L和U不如矩阵A稀疏，所以使用LU因式分解时有效性会下降一些。

矩阵的带宽越大，即，所找到的非零值离主对角线越远，这个问题就会越明显。

有时，我们仅想要一个近似的因式分解B=L*U，其中，B与A足够接近以至于C=inv(B)*A与恒等矩阵差别不大，即，C的最大和最小本征值之比小于A的最大和最小本征值之比。

在这里，B被称为预处理量，用于最优化法以及解决某类线性系统的问题。

当我们进行一个不完全LU因式分解时，我们仅计算L和U 中对应于A中非零元素的元素，或者我们另有选择，忽略绝对值小于给定值的元素。

下列代码说明了不完全LU因式分解的用法。

除了设定某些元素等于零之外，使B=A。

令远离对角线的一些元素等于零。

但是，看一下B的和近似因式分解的本征值。

4.4. QR分解
P是置换矩阵，Q是正交矩阵，R是上三角阵，通过调用指令“qr”来进行因式分解A*P=Q*R。

如果想得到A=QR分解（即P=1），指令是：
4.5. Cholesky分解
如果A是Hermetian矩阵（即A=A’）那么我们知道所有本征值都是实数。

另外，如果所有本征值都大于零，对于所有向量x都有x’*A*x>0，那么我们就说A是正定的。

在这里，进行Cholesky分解（即A=R’*R）是可能的，其中R 是上三角阵。

这就相当于写出A=L*L’，其中L是下三角阵。

首先，我们使用下面的正定矩阵。

对于稀疏矩阵，我们可以进行不完全Cholesky分解，它能给出用做预处理量的近似的因式分解，且不失稀疏性。

在此特殊的情况下，对于高度结构化的矩阵来说，不完全因式分解与完全因式分解是一样的。

4.6. 奇异值分解
只有在方阵的情况下定义本征值和本征向量才有意义。

非方阵本征值的概念的普遍化通常也很有用。

(m×n)型的矩阵A的奇异值分解（SVD）被定义为A=U*D*V’，其中，D是一个包含奇异值的(m×n)型对角阵，U是一个包含右本征向量的(m×m)型矩阵，V’是一个包含左本征向量的(n×n)型伴随（转置与共轭）阵。

在MATLAB中，使用指令“svd”可以进行奇异值分解。