009稀疏矩阵

下载 第9章 稀 疏 矩 阵 在许多问中提到了含有大量 0元素的矩阵，这样的矩阵称为稀疏矩阵。比如求解普通或 者部分微分方程的数值解。为了节省存储空间和计算时间， M AT L A B考虑到矩阵的稀疏性， 在对它运算时有特殊的命令。 9.1 矩阵为什么稀疏 一个稀疏矩阵中有许多元素等于零，这便于矩阵的计算和保存。如果 M AT L A B把一个矩 阵当作稀疏矩阵，那么只需在 m×3的矩阵中存储 m个非零项。第1列是行下标，第 2列是列下 标，第3列是非零元素值，不必保存零元素。如果存储一个浮点数要 8个字节，存储每个下标 要4个字节，那么整个矩阵在内存中存储需要 1 6×m个字节。 ■ 例9 . 1 A = e y e ( 1 0 0 0 ) ; 得到一个1 0 0 0×1 0 0 0的单位矩阵，存储它需要8 Mb空间。如果使用命令： B = s p e y e ( 1 0 0 0 ) ; 用一个1 0 0 0×3的矩阵来代表，每行包含有一个行下标、列下标和元素本身。现在只需 1 6 K b的空间就可以存储1 0 0 0×1 0 0 0的单位矩阵，它只需要满单位矩阵的 0 . 2 %存储空间。对于 许多的广义矩阵也可这样来作。 稀疏矩阵的计算速度更快，因为 M AT L A B只对非零元素进行操作，这是稀疏矩阵的第二 个突出的优点。 ■ 例9 . 2 假设矩阵A、B和例9 . 1中的矩阵一样。计算 2*A需要一百万次的浮点运算，而计算 2*B只 需要2 0 0 0次浮点运算。 因为M AT L A B不能自动创建稀疏矩阵，所以要用特殊的命令来得到稀疏矩阵，在下一节 中将给出这些命令。前面章节中的算术和逻辑运算都适用于稀疏矩阵。 9.2 创建和转换稀疏矩阵 在M AT L A B中，用命令s p a r s e 来创建一个稀疏矩阵。 命令集8 7 创建稀疏矩阵 s p a r s e ( A ) 由非零元素和下标建立稀疏矩阵 A。如果A已是一个稀疏 矩阵，则返回A本身。 s p a r s e ( m , n ) 生成一个m×n的所有元素都是0的稀疏矩阵。 ■ ■ s p a r s e ( u , v , a ) 生成一个由长度相同的向量 u，v和a定义的稀疏矩阵。其 中u和v是整数向量， a是一个实数或者复数向量。 (ui, vi) 对应值 ai，如果 a中有零元素，则将这个元素排除在外。 稀疏矩阵的大小为m a x (u)×m a x (v)。 s p a r s e ( u , v , a , m , n ) 生成一个m×n的稀疏矩阵，(ui, vi)对应值ai。向量u，v和a 必须长度相同。 s p a r s e ( u , v , a , m , n , 生成一个m×n的含有n z m a x个非零元素的稀疏矩阵。 (ui, n z m a x ) vi)对应值ai。n z m a x的值必须大于或者等于向量u和v的长度。 f i n d ( x ) 返回向量x中非零元素的下标。如果 x=X是一个矩阵，那 么X的向量就作为一个长向量来考虑。 [ u , v ] = f i n d ( A ) 返回矩阵A中非零元素的下标。 [ u , v , s ] = f i n d ( A ) 返回矩阵A中非零元素的下标。用向量s中元素的值及u和v中 相应的下标，实际上就是向量u、v和s作为命令s p a r s e的参数。 s p c o n v e r t ( D ) 将一个有三列的矩阵转换成一个稀疏矩阵。 D中的第1列作 为行的下标，第2列作为列的下标，最后一列作为元素值。 而且可以使用命令 f u l l将稀疏矩阵转换成一个满矩阵。 命令集8 8 转换成满矩阵 f u l l ( S ) 将稀疏矩阵S转换成一个满矩阵。 ■ 例9 . 3 (a) 创建一个5×5的单位矩阵： A = e y e ( 5 ) 将矩阵A转换成稀疏矩阵B： (b) 假设M AT L A B中给出如下的向量： 这样就有了行向量，但是也可使用列向量。运行命令S m a t r i x = s p a r s e ( i n d 1 , i n d 2 , n u m b e r ) ， 结果为： 1 2 4M ATLAB 5 手册 下载 其中有去掉了两个零元素。将这个矩阵转换成满矩阵，输入： F u l l m a t r i x = f u l l ( S m a t r i x ) 得到的结果为： 注意，稀疏矩阵和得到的满矩阵的大小是分别是由 i n d 1和i n d 2中最大元素值确定的，即 使相应的值是零，并且在列出的稀疏矩阵中去掉这个值。 输入命令w h o s 可得到： 可以看出虽然两个矩阵的大小相同，但是其中稀疏矩阵需要的存储空间更小些。 (c) 在处理稀疏矩阵时f i n d 命令很有用。命令对于稀疏矩阵或者满矩阵都返回相同的结果。 返回得到的三个向量直接用来重新创建一个稀疏矩阵。令 S m a t r i x定义在( b )中，运行命令： 得到的结果为： 用下面命令得到的矩阵和 ( b )中得到的矩阵是不一样的： 第9章 稀 疏 矩 阵 1 2 5 下载 ■ 9.3 稀疏矩阵运算 M AT L A B中对满矩阵的运算和函数同样可用在稀疏矩阵中。结果是稀疏矩阵还是满矩阵， 这取决于运算符或者函数及下列的操作数： • 当函数用一个矩阵作为输入参数，输出参数为一个标量或者一个给定大小的向量时，输 出参数的格式总是返回一个满阵形式，如命令 s i z e 。 • 当函数用一个标量或者一个向量作为输入参数，输出参数为一个矩阵时，输出参数的格式 也总是返回一个满矩阵，如命令e y e 。还有一些特殊的命令可以得到稀疏矩阵，如命令s p e y e 。 • 对于单参数的其他函数来说，通常返回的结果和参数的形式是一样的，如 d i a g 。 • 对于双参数的运算或者函数来说，如果两个参数的形式一样，那么也返回同样形式的结 果。在两个参数形式不一样的情况下，除非运算的需要，均以满矩阵的形式给出结果。 • 两个矩阵的组和[A B]，如果A或B中至少有一个是满矩阵，则得到的结果就是满矩阵。 • 表达式右边的冒号是要求一个参数的运算符，遵守这些运算规则。 • 表达式左边的冒号不改变矩阵的形式。 ■ 例9 . 4 假设有： 这是一个5×5的单位满矩阵和相应的稀疏矩阵。 (a) C = 5 *B，结果为： 这是一个稀疏矩阵。 (b) D = A + B ，给出的结果为： 这是一个满矩阵。 (c) x = B \ h ，结果为： 这是一个满向量。 1 2 6M ATLAB 5 手册 下载 ■ 有许多命令可以对非零元素进行操作。 命令集8 9 矩阵的非零元素 n n z ( A ) 求矩阵A中非零元素的个数。它既可求满矩阵也可求稀疏矩阵。 s p y ( A ) 画出稀疏矩阵 A中非零元素的分布。也可用在满矩阵中，在 这种情况下，只给出非零元素的分布。 s p y ( A , c s t r , s i z e ) 用指定的颜色 c s t r(见表1 3 - 1 )和在s i z e规定的范围内画出稀疏 矩阵A中非零元素的分布。 n o n z e r o s ( A ) 按照列的顺序找出矩阵A中非零的元素。 s p o n e s ( A ) 把矩阵A中的非零元素全换为1。 s p a l l o c ( m , n , 产生一个m×n阶只有n z m a x个非零元素的稀疏矩阵。这样可以 n z m a x ) 有效地减少存储空间和提高运算速度。 n z m a x ( A ) 给出为矩阵A中非零元素分配的内存数。不一定和 n n z ( A ) 得 到的数相同；参见s p a r s e 或者s p a l l o c 。 i s s p a r s e ( A ) 如果矩阵A是稀疏矩阵，则返回1；否则返回0。 s p f u n ( f c n , A ) 用A中所有非零元素对函数 f c n求值，如果函数不是对稀疏矩 阵定义的，同样也可以求值。 s p f u n ( A ) 求稀疏矩阵 A的结构秩。对于所有的矩阵来说，都有 s p r a n k ( A ≥r a n k ( A ) 。 ■ 例9 . 5 用下面的命令定义稀疏矩阵： 现在创建一个大矩阵： Big=kron(A, A) 这个矩阵B i g是什么样子呢？K r o n e c k e r张量积给出一个大矩阵，它的元素是矩阵 A的元素 之间可能的乘积。因为参量都是稀疏矩阵，所以得到的矩阵也是一个稀疏矩阵。可以用命令 w h o s 和i s s p a r s e 来确认一下。 查看矩阵B i g的结构图，可输入s p y ( B i g ) ，结构如图9 - 1所示。 图9-1 用s p y命令显示的矩阵结构图 第9章 稀 疏 矩 阵 1 2 7 下载 可以看出B i g是一个块双对角矩阵。 9.4 稀疏矩阵的特例 MATLAB中有四个基本稀疏矩阵，它们是单位矩阵、随机矩阵、对称随机矩阵和对角矩阵。 命令集9 0 单位稀疏矩阵 s p e y e ( n ) 生成n×n的单位稀疏矩阵。 s p e y e ( m , n ) 生成m×n的单位稀疏矩阵。 命令speye(A) 得到的结果和s p a r s e ( e y e ( A ) ) 是一样的，但是没有涉及到满阵的存储。 命令集9 1 随机稀疏矩阵 s p r a n d ( A ) 生成与A有相同结构的随机稀疏矩阵，且元素服从均匀分布。 s p r a n d ( m , n , d e n s ) 生成一个m×n的服从均匀分布的随机稀疏矩阵，有d e n s×m× n个非零元素，0≤d e n s≤1。参数d e n s是非零元素的分布密度。 s p r a n d ( m , n , d e n s , 生成一个近似的条件数为1 /rc、大小为m×n的随机稀疏矩 r c ) 阵。如果rc=rc是一个长度为 l≤l ( m i n (m, n) )的向量，那么 矩阵将rci作为它l个奇异值的第一个，其他的奇异值为 0。 s p r a n d n ( A ) 生成与A有相同结构的随机稀疏矩阵，且元素服从正态分布。 s p r a n d n ( m , n , d e n s , 生成一个m×n的服从正态分布的随机稀疏矩阵，和sprand r c ) 一样。 s p r a n d s y m ( S ) 生成一个随机对称稀疏矩阵。它的下三角及主对角线部分 与S的结构相同，矩阵元素服从正态分布。 s p r a n d s y m ( n , d e n s ) 生成一个 m×n的随机对称稀疏矩阵。矩阵元素服从正态 分布，分布密度为d e n s。 s p r a n d s y m ( n , d e n s , 生成一个近似条件数为1 /rc的随机对称稀疏矩阵。元素以0 r c ) 对称分布，但不是正态分布。如果 rc=rc是一个向量，则矩 阵有特征值rci。也就是说，如果 rc是一个正向量，则矩阵 是正定矩阵。 s p r a n d s y m ( n , d e n s , 生成一个正定矩阵。如果k= 1，则矩阵是由一正定对称矩 r c , k ) 阵经随机 J a c o b i旋转得到的，其条件数正好等于 1 /rc；如 果k= 2，则矩阵为外积的换位和，其条件数近似等于 1 /rc。 s p r a n d s y m ( S , d e n s , 生成一个与矩阵S结构相同的稀疏矩阵，近似条件数为1 /rc。 r c , 3 ) 参数d e n s被忽略，但是这个参数在这个位置以便函数能确认 最后两个参数的正确与否。 ■ 例9 . 6 (a) 假设有矩阵A： 1 2 8M ATLAB 5 手册 下载 ■ 输入R a n d o m = s p r a n d n ( A ) ，可得到随机稀疏矩阵： 矩阵中随机数的位置和矩阵A中非零元素的位置相同。 (b) 对于( a )中的矩阵A，输入： B = s p r a n d s y m ( A ) 结果为： 这是一个用矩阵A的下三角及主对角线部分创建的对称矩阵，在非零元素的位置用随机数 作为元素值。 用命令s p d i a g s可以取出对角线元素，并创建带状对角矩阵。假设矩阵 A的大小为m×n， 在p个对角线上有非零元素。 B的大小为m i n (m×n)×p，它的列是矩阵A的对角线。向量d的长 度为p，其整型分量给定了A的对角元： di<0 主对角线下的对角线。举例如 di=－1，用第一个下对角线 di=0 用主对角线 di>0 主对角线上的对角线 命令集9 2 对角稀疏矩阵 [ B , d ] = s p d i a g s ( A ) 求出A中所有的对角元，对角元保存在矩阵 B中，它们的 下标保存在向量d中。 s p d i a g s ( A , d ) 生成一个矩阵，这个矩阵包含有矩阵A中向量d规定的对角元。 s p d i a g s ( B , d , A ) 生成矩阵A，用矩阵B中的列替换d定义的对角元。 A = s p d i a g s ( B , d , m , n ) 用保存在由d定义的B中的对角元创建稀疏矩阵A。 例11 . 4给出了如何使用s p d i a g s 命令来解普通微分方程组。 9.5 系数阵为稀疏矩阵的线性方程组 在许多实际应用中要保留稀疏矩阵的结构，但是在计算过程中的中间结果会减弱它的稀 疏性，如 L U分解。这就会导致增加浮点运算次数和存储空间。为了避免这种情况发生，在 第9章 稀 疏 矩 阵 1 2 9 下载 ■ M AT L A B中用命令对矩阵进行重新安排。这些命令都列在下面的命令集 9 3中。通过h e l p 命令 可以得到每个命令更多的帮助信息，也可见 h e l p d e s k 。 命令集9 3 矩阵变换 c o l m m d ( A ) 返回一个变换向量，使得矩阵 A列的秩为最小。 s y m m m d ( A ) 返回使对称矩阵秩为最小的变换。 s y m r c m ( A ) 矩阵A的C u t h i l l - M c K e e逆变换。矩阵A的非零元素在主对角线附近。 c o l p e r m ( A ) 返回一个矩阵A的列变换的向量。列按非零元素升序排列。有时这 是L U因式分解前有用的变换：lu(A(:, j)) 。如果A是一个对称矩 阵，对行和列进行排序，这有利于C h o l e s k y分解：chol(A(j, j)) 。 r a n d p e r m ( n ) 给出正数1，2，. . .，n的随机排列，可以用来创建随机变换矩阵。 d m p e r m ( A ) 对矩阵A进行D u l m a g e - M e n d e l s o h n分解，输入help dmperm可 得更多信息。 ■ 例9 . 7 创建一个秩为4的变换矩阵，可输入： 一旦运行p e r m = r a n d p e r m ( 4 )，就会得到： 给出的变换矩阵为： 如果矩阵A为： 输入命令： 运行结果为： 1 3 0M ATLAB 5 手册 下载 有两个不完全因式分解命令，它们是用来在解大线性方程组前进行预处理的。用 h e l p d e s k 命令可得更多信息。 命令集9 4 不完全因式分解 c h o l i n c ( A , o p t ) 进行不完全C h o l e s k y分解，变量o p t取下列值之一： d r o p t o l 指定不完全分解的舍入误差， 0给出完全分解。 m i c h o l 如果m i c h o l = 1 ，则从对角线上抽取出被去掉的元素。 r d i a g 用s q r t ( d r o p t o l *n o r m ( X ( : , j ) ) ) 代替上三角分 解因子中的零元素， j为零元素所在的列。 [ L , U , P ] = 返回矩阵X的不完全分解得到的三个矩阵 L、U和P，变量o p t取 l u i n c ( X , o p t ) 下列值之一： d r o p t o l 指定分解的舍入误差。 m i l u 改变分解以便从上三角角分解因子中抽取被去掉的列元素。 u d i a g 用d r o p t o l 值代替上三角角分解因子中的对角线上的 零元素。 t h r e s h 中心极限。 解稀疏线性方程组既可用左除运算符解，也可用一些特殊命令来解。 命令集9 5 稀疏矩阵和线性方程组 s p p a r m s ( k e y s t r , o p ) 设置稀疏矩阵算法的参数，用 help spparms可得详细信息。 s p a u g m e n t ( A , c ) 根据[ c* l A; A’ 0 ]创建稀疏矩阵，这是二次线性方程组的最 小二乘问题。参见7 . 7节。 s y m b f a c t ( A ) 给出稀疏矩阵的 C h o l e s k y和L U因式分解的符号分解因子。 用help symbfact可得详细信息。 稀疏矩阵的范数计算和普通满矩阵的范数计算有一个重要的区别。稀疏矩阵的欧几里德 范数不能直接求得。如果稀疏矩阵是一个小矩阵，则用 n o r m ( f u l l ( A ) ) 来计算它的范数； 但是对于大矩阵来说，这样计算是不可能的。然而 M AT L A B可以计算出欧几里德范数的近似 值，在计算条件数时也是一样。 命令集9 6 稀疏矩阵的近似欧几里德范数和条件数 n o r m e s t ( A ) 计算A的近似欧几里德范数，相对误差为 1 0－6。 n o r m e s t ( A , t o l ) 计算A的近似欧几里德范数，设置相对误差 t o l，而不用缺省时 的1 0－6。 第9章 稀 疏 矩 阵 1 3 1 下载 ■ [ n r m , n i t ] = 计算近似n r m范数，还给出计算范数迭代的次数 n i t。 n o r m e s t ( A ) c o n d e s t ( A ) 求矩阵A条件数的1 -范数中的下界估计值。 [ c , v ] = 求矩阵A的1 -范数中条件数的下界估计值 c和向量v，使得 c o n d e s t ( A , t r ) | |Av| | = ( | |A| |·| |v| | ) / c。如果给定 t r，则给出计算的过程。 t r= 1， 给出每步过程；t r=－1，给出商c / r c o n d ( A ) 。 ■ 例9 . 8 假设给出： 用n o r m A p p r o x = n o r m e s t ( S p r s )计算出： 用t h e N o r m = n o r m ( f u l l ( S p r s ) )得： 为了找到它们之间的差别，计算 d i f f e r e n c e = t h e N o r m - n o r m A p p r o x ，得： 在许多应用中，n o r m e s t 计算得到的近似值是一个很好的近似欧几里德范数，它的计算 步数要比n o r m 要少得多；可参见7 . 6节。 用e t r e e 命令来找到稀疏对称矩阵的消元树，用向量 f来描述消元树，还可用 e t r e e p l o t 命令画出来。元素 fi是矩阵的上三角C h o l e s k y分解因子中 i行上第1非零元素的列下 标。如果有非零元素，则 fi= 0。消元树可以这样来建立： 节点i是fi的孩子，或者如果 fi= 0，则节点 i是树的根节点。 命令集9 7 矩阵的消元树 e t r e e ( A ) 求A的消元树向量 f，这个命令有可选参数；输入 h e l p e t r e e 获取帮助。 e t r e e p l o t ( A ) 画出向量f定义的消元树图形。 t r e e p l o t ( p , c , d ) 画出指针向量 p的树图形，参数 c和d分别指定节点的颜色 和分支数。e t r e e p l o t 可以调用这个命令。 t r e e l a y o u t 显示树的结构，t r e e p l o t 可以调用这个命令。 ■ 例9 . 9 假设有对称稀疏矩阵B： 1 3 2M ATLAB 5 手册 下载 ■ 运行命令b t r e e = e t r e e ( B )，结果为： 开始的数字2不难理解，它是矩阵的第 1列上第1个非零元素的行数，它决定了在 C h o l e s k y 分解因子的第1行第2列处有一个非零元素。当缩减第 1列的元素时就得到第 2列的数字5。B在 缩 减 后 ， 在 ( 5 , 2 )位 置 的 元 素 是 非 零 的 ， 这 样 消 元 树 向 量 中 第 2个 元 素 的 值 为 5。 s p y ( c h o l ( B ) ) 给出了C h o l e s k y分解因子的结构图，如图9 - 2所示： 图9-2 Cholesky分解结构图 图9-3 矩阵B的消元树 这个向量消元树可以这样来建立：上三角中只有一行有非零元素，节点 8，因此这就是树 唯一的根。节点 1是节点2的孩子，节点2和3是节点5的孩子，而节点5是节点6的孩子。节点4 第9章 稀 疏 矩 阵 1 3 3 下载 和6是节点7的孩子，而节点7又是节点8的孩子，即根的孩子。 命令e t r e e p l o t ( B ) 给出了树的结构图，如图9 - 3所示。 消元树的形状取决于列和行序，它可以用来分析消元过程。 用g p l o t 命令可以画出坐标和矩阵元素间的联系图形。必须在n×2的矩阵中给出n个坐标， 矩阵的每一行作为一个点。这样就创建出点点之间连接的 n×n矩阵，如果点 4连接到点8，则 (4, 8)的值为1。由于是一个大矩阵，而且非零元素较少，所以它应该被建成稀疏矩阵。 这个图可以说明网络问题，如传递问题。它还包含有线性方程组中未知量之间的相关信息。 命令集9 8 网络图形 g p l o t ( A , K ) 如果矩阵A的a(i, j)不为0，则将点ki连接到点kj。K是一个n× 2的坐标矩阵，A是一个n×n的关联矩阵。 g p l o t ( A , K , s t r ) 用字符串s t r给定的颜色和线型画出的同上图形。字符串 s t r的 取值参见表1 3 - 1。 [ X , A ] = u n m e s h ( E ) 求网格边界矩阵E的L a p l a c e矩阵A和网格点的坐标矩阵X。 ■ 例9 . 1 0 假设有下面的坐标矩阵K和关联矩阵A： 矩阵A在稀疏化后，用命令g p l o t ( A , K ) 画出图9 - 4，给出了点 (0, 1)和点(4, 1)之间所有可 能的路径。 图9-4 使用g p l o t的一个例子 1 3 4M ATLAB 5 手册 下载 ■ ■ 第9 章稀疏矩阵 9.1 矩阵为什么稀疏 9.2 创建和转换稀疏矩阵 9.3 稀疏矩阵运算 9.4 稀疏矩阵的特例 9.5 系数阵为稀疏矩阵的线性方程组