数值模拟导论-第四讲 线性稀疏矩阵的直接解法 Luca Daniel 感谢Deepak Ramaswamy, Michal Rewienski,Karen Veroy and Jacob White ?回顾LU分解法 ?稀疏矩阵 ?三角矩阵分解 ?稀疏矩阵数据结构 概述 SMA-HPC ?2003 MIT —珩架和节点，电阻网，3d热流 ―一般的稀疏矩阵分解 ―填充和重排列 —图表逼近 —散布 LU分解基础 SMA-HPC ?2003 MIT 分解 图片 上图便是LU分解的图形表示。第一步，用第一个方程消去第二到第四方程中 的x1。这一过程我们用除第一行外的各行分别减去第一行乘以某个比例因子，从 而使系数a21，a31，a41变为零。再用比例因子（又称之为乘子）代替这些零 位。对于第二行，乘子是a21/a11，因为第二行减去第一行乘以a21/a11，a21位 正好为零。由于在消去过程中a22，a23，a24的值也会随之改变，因此我们将他 们变成蓝色。同样在消去a31和a41的过程中，a31和a41也被他们的乘子所代 替。在这一过程中第三行其余的位置的值也会随之改变，因此也将他们变为蓝 色。 用同样的方法处理第二行。计算消去第三行和第四行中x2的乘子，并且用这 些乘子代替出现的零。并且注意在消去过程中改变的量，将他们改为绿色。最后 一步，便是用第三行消去第四行中的x3，更新第四行的各个位置，并且将a44变 为粉红色。 我们可以看到乘子在代替矩阵中的零的位置之后，在消去过程中他们并没有 改变。 LU分解基础 SMA-HPC ?2003 MIT 矩阵分解 算法 for i=1 到n-1 {每一行 for j=i+1 到n {每一要消去的目标行 ji ji ii M M M = for k=i+1 到n {对角元后的元素 jk jk jiik M MMM← ? 乘子 } } } 对角元 LU分解基础 SMA-HPC ?2003 MIT 矩阵分解 对角占优矩阵的性质 A）对一个对角占优的矩阵进行LU分解时不会 产生零对角元。 B）严格对角占优矩阵经过LU分解它的各个位 置上的值增加不会超过 (1) 2 n? 。 定理：在对严格对角占优的矩阵进行高斯消元时 不会产生零对角元。 证明：1）求出第一步消元后的矩阵。 2）考察（n－1）×（n－1）的次矩阵。 仍然是完全对角占优矩阵。 第一步 第一步消元后的第二行 由此得出 稀疏矩阵 SMA-HPC ?2003 MIT 应用 空间珩架 空间珩架节点矩阵 未知量：节点位置 方程：合力＝0 稀疏矩阵 SMA-HPC ?2003 MIT 应用 电阻网 未知：节点电压 方程：电流和＝0 电阻网是一种特殊情况，它的数学模型是偏微分方程。（我们将在以后 学习到）我们现在考察一下这个节点矩阵并注意矩阵中的非零数。从一个 4×4的例子中我们可以很容易看出矩阵的特点。一个4×4系统它的节点矩阵 如下： 三角带取决于结触结点之间沿栅格行的相互作用。与三角带相距4的非零数 是由行节点与三角带的耦合产生的。 稀疏矩阵 SMA-HPC ?2003 MIT 应用 电阻网 100×10电阻网矩阵的非零带 稀疏矩阵 SMA-HPC ?2003 MIT 应用 四方体中的温度场 电路 模型 稀疏矩阵 SMA-HPC ?2003 MIT 三角带矩阵 矩阵形式 稀疏矩阵 SMA-HPC ?2003 MIT 三角带矩阵 高斯算法 for i=1 到n-1 {每一行 for j=i+1 到n {每一要消去的目标行 ji ji ii M M M = 对角元 for k=i+1 到n {对角元后的元素 jk jk ji ik M MMM← ? 乘子 } } } 需要N步运算 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵的填充 例子 电阻例子 节点矩阵 此矩阵为对 称且对角占 优矩阵 回想第二章的内容，节点矩阵的维数可以通过计算电阻器来确定，请看下面： 这个电阻占据节点矩阵的如下区域 也可以理解为：G ii 就等于解点i处的电阻率（一个以上电阻）之和 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵的填充 例子 矩阵的非零数据结构LU分解后的矩阵形式 X都不为零 在LU分解的过程中，基准行以外的各行要减去基准行乘以一个因子。既然 在同一列中这两行没有必要有非零数，相减的目的就是要在目标行中添加非零 数。 请通过下面的简单例子，思考LU分解。 结果是 我们看到分解后的矩阵的右下角有一个非零位。而在原先的矩阵中他是 零。这种将原先的零位变为非零位的运算称之为填充。 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵填充 例2 填充的传递 第一步的填充影响第二步的填充 在例子中，这个4×4的矩阵有7个零位。在LU分解的过程中，5个零位变 为了非零。深入了解填充我们会发现：LU分解的第一步即用第二行减去第一 行乘以一个系数，填充了第三列和第四列。当用第三行和第四行减去第二行 乘以一个系数后，那么第三行和第四就产生了二级填充。 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵的填充 重排序 将结点重新排序可以减少填充 －可以保留某些性质（对称，对角占优） －等同于交换行和列 填充 未填充 在文中的节点方程式中，将结点重新编号是减小填充的简单的方法，因为结 点的编号可以任意选择。要记住，在节点方程中重新编号节点与在矩阵中交 换两行或列的位置是相对应的。 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵的填充 重排序 在什么地方可以产生填充？ 已经分解 可能产生填充 的位置 乘子 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵的填充 重排序 Markowitz重排序 For i=1 到n 找到最小的Markowitz乘积的对角元j。 交换第j行和第i行，第j列和第i列交换。 分解新的第i行，并确定填充。 结束 理想算法 为了能理解Markowitz的算法，我们要考虑算法所花费的机时。第一步 我们要确定具有最小的Markowitz乘积的对角元。这一步的耗时为K×N步 运算其中K是每一行非零位数的平均值。第二步就是在分解中交换行和 列。一个好的数据结构将会使交换很省时。 第三步就是分解重新排列的矩阵并且填充矩阵。如果矩阵本来非常稀疏 那么第三步将会很省时。 接下来又必须计算分解后的矩阵的Markowitz乘积，这有又需要K× （N-1）步运算。继续，依此类推，只计算Markowitz乘积就需要约 1/2KN2步的运算。 因此，我们要改善这种情况，方法便是计算矩阵每经过一次分解后 Markowitz乘积的微小变化。通过观看矩阵图我们可以很容易的找到这一 优化技巧。 ?与在节点公式中节点重排列相协调。 ?减少搜索耗时 ?保留矩阵原有性质 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵的填充 重排序 为什么只讨论对角元？ －对角占优 －对称 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵的填充 填充之后矩阵的样式 非常稀疏 非常稀疏稠密 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵的填充 未分解的任意矩阵 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵的填充 未分解的任意矩阵 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵图 意义 对称矩阵的结构和据阵图 ·矩阵每一行一个节点 ·每一个非对角元对应一条边线。 如果矩阵在结构上是对称的（ 当且仅当 在图中矩阵的每一行有一个节点 如果aij≠0，在节点i和j之间有一个边界。 矩阵图有两个重要的性质： 1)节点次数的平方就是Markowitz乘积 2）据阵图在完成一次矩阵分解后可以进行简单变化，就可适应矩阵的 变化。 矩阵图使得在分解对称矩阵时两步法变得非常有效。第一步用矩阵 图找出一个Markowitz乘积最小的矩阵排序。然后，在新排列的矩阵中 进行数值分解。 ），矩阵可以 与一个简单的图表相联系。 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵图 MarkowitZ乘积 Markowitz乘积＝（节点次数的平方） 第i个节点次数的平方等于Markowitz乘积，它可以由图中很容易看出。节点 的次数就是此节点与其他节点相连的边界，并且每个边界代表了本行中非对 角的非零位和本列中非对角的非零位。因此，行中非对角的非零位和列中非 对角的非零位的乘积就等于节点次数的平方。 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵图 矩阵分解 LU分解的第一步 ·删掉与对角行相对应的节点。 ·连接矩阵图的边线。 每一步LU分解需要一步浮点运算并产生一个简化矩阵，如下： 完成了第i步分解之后，未分解的矩阵就被缩小为（i－1）×（i －1）大小的矩阵，并且如果有填充这个矩阵还会变稠密。矩阵图 可以用来表示未分解矩阵部分的非零位，但是有两点必须改变： 1）在未分解矩阵部分少了一行之后，矩阵图也要相应减少一个 节点。 2）必须加入与填充对应的矩阵图边界。 在图中，我们通过一个例子显示了在经过一步LU分解之后矩阵 图是怎样变化的。我们可以明显得看到，在矩阵中如果i行被消去， 那么在矩阵图中节点i也必须除去。另外，所用与节点i相邻的节点（ 相邻节点通过边界相联）将会通过增加必要的边界来使他们彼此相 邻。这个增加的边界就是填充。 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵图 例子 据阵图 Markowitz乘积 ＝（节点次数） 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵图 例子 交换第2和第1行 下面的例子是关于不产生填充的因子。 对角矩阵 更加平行的对角矩阵的另一顺序。 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵图 电阻网的例子 未知：节电电压 方程：电流和等于零 电阻网是一种特殊情况，它的数学模型是偏微分方程。（我们将在 以后看到） 我们现在考察一下这个节点矩阵并注意矩阵中非零数。从一个4×4 的例子中我们可以很容易看出矩阵的特点。一个4×4系统它的节点矩阵 如下： 三角带取决于结触结点沿栅格行的相互作用。与三角带相距4的非 零数是由于row coupling与三角带的耦合产生的。 稀疏矩阵 SMA-HPC ?2003 MIT 矩阵图 网格的例子 分解一个M×M的矩阵需要花费多少机时？ 推想中间列最后是否被消去？ 要分解一个与类似由电阻器序列产生的M×M个删格相对应的 矩阵，一个快速的方法就是考察矩阵图。如果将中心列节点M排在矩阵 图中的最后，那么他们最终将按照图中显示进行连接。因此，产生了一个 最终与稠密矩阵相对应的矩阵图。 虽然要得到这个问题的精确解已经超出了本课程的内容但我们可以看到 ，既然分解这个稠密矩阵需要 步运算，所以这表明分解这个由M×M个 步运算。 删格组成的矩阵也需要 矩阵分解的方法 SMA-HPC ?2003 MIT 1）用非数值的元素代替原来的矩阵 对角占优或对称正定 2）用矩阵图确定矩阵的次数 要用到许多矩阵图操作的技巧 3）形成存储填充矩阵的数据结构 加入了许多额外的非零元 4）将原来的数值代入这个数据结构中并分 解。 必须认真的组织计算 稀疏矩阵 SMA-HPC ?2003 MIT 稀疏矩阵的数 据结构 行指针向量 每一行的数据排列 为了有效的储存一个稀疏矩阵，需要一种只记录矩阵非零元素的 数据结构。一种简单的方法是基于每一行至少有一个非零元素。因 此，我们对每一行建立一对数组，这个数组与矩阵非零位的值和位所 在的列相对应。请看下面的例子： 这个矩阵的数据结构是 我们可以看到这里没有储存任何的零位。 稀疏矩阵 SMA-HPC ?2003 MIT 稀疏矩阵的数据结 构 冗余数据 从目标行j中消去原始行i。 考察j行所有位，我们发现其中3个位置与i行相对应。 第i行 第j行 为了有效的储存一个稀疏矩阵，需要一种只记录矩阵非零元素的数 据结构。一种简单的方法是基于每一行至少有一个非零元素。因此，我 们对每一行建立一对数组，这个数组与矩阵非零位的值和位所在的列相 对应。请看下面的例子： 这个矩阵的数据结构是 我们可以看到这里没有储存任何的零位。 稀疏矩阵 SMA-HPC ?2003 MIT 稀疏矩阵的数据结构 冗余数据 每一个多余的数据都会占用不必要的内存空间 为了有效的储存一个稀疏矩阵，需要一种只记录矩阵非零元素的数 据结构。一种简单的方法是基于每一行至少有一个非零元素。因此，我 们对每一行建立一对数组，这个数组与矩阵非零位的值和位所在的列相 对应。请看下面的例子： 这个矩阵的数据结构是 我们可以看到这里没有储存任何的零位。 稀疏矩阵 SMA-HPC ?2003 MIT 稀疏矩阵的数据结构 避免冗余数据的扩散 第j行 1）考察j行所有的元素，将他们扩展成为长度为n的向量 2）用标定数组储存需要的元素 总结 SMA-HPC ?2003 MIT ?LU分解和对角占优 －无数据元的因子 ?稀疏矩阵 －压杆，电阻网，3D热流 ?对角矩阵的分解 －O(N)步分解 ?稀疏矩阵分解概述 －Markowitz重排步以达到最小填充 ?基于表格的方法 －分解和填充 －对降低用高斯消去解稀疏矩阵的复杂性是有帮助的