稀疏矩阵运算- MATLAB和Simulink金宝app - 金宝app,下载188bet金宝搏,金宝搏官方网站

稀疏矩阵的操作

操作效率

计算复杂度
算法细节

计算复杂度

稀疏操作的计算复杂度与nnz，矩阵中非零元素的个数。计算复杂度也线性地依赖于行大小米和列的大小n矩阵，但与乘积无关m * n，即零元素和非零元素的总数。

一些相当复杂的操作，如稀疏线性方程的解，其复杂性涉及排序和填充等因素，这在前一节中已经讨论过。然而，一般来说，稀疏矩阵运算所需的计算机时间与非零量上的算术运算的数量成正比。

算法细节

稀疏矩阵按照以下规则进行计算:

接受矩阵并返回标量或常量大小向量的函数总是生成完整存储格式的输出。例如,大小功能始终返回完整矢量，无论是完整还是稀疏。
接受标量或向量并返回矩阵的函数，例如0，的，兰德,眼睛，始终返回完整结果。这是必要的，以避免意外引入稀疏性。稀疏的模拟0 (m, n)仅仅是稀疏(m, n)．的稀疏类比兰德和眼睛是sprand和speye,分别。这个函数没有稀疏模拟的．
接受矩阵并返回矩阵或向量的一元函数保留操作数的存储类。如果年代是稀疏矩阵吗胆固醇(S)也是一个稀疏矩阵，并且诊断接头(S)是一个稀疏向量。列函数，例如马克斯和总和也返回稀疏向量，即使这些向量可以完全非零。这条规则的重要例外是稀疏的和完整的功能。
如果两个操作数都是稀疏的，则二元操作符产生稀疏结果，如果两个操作数都是满的，则产生完整结果。对于混合操作数，除非操作保持稀疏性，否则结果是满的。如果年代是稀疏的,F满了,然后S + F，S * F,F \ S是完整的,而美国* F和S&F是稀疏的。在某些情况下，即使矩阵只有很少的零元素，结果也可能是稀疏的。
的矩阵连接猫函数或方括号对混合操作数产生稀疏结果。

稀疏矩阵的行和列的排列年代可以用两种方式表示:

一个置换矩阵P行动在一排排年代作为P * S或者列为* P '．
一个排列向量p，它是一个包含的排列的完整向量1: n，作用于一排排年代作为: S (p)，或列为S (: p)．

例如:

p = [1 3 4 2 5] I = eye(5,5);P = I(P，:) e = ones(4,1);S = diag(11:11:55) + diag(e,1) + diag(e，-1)

p = 1 2 3 4 5 p = 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 S = 22 11 1 0 0 0 1 1 0 0 0 1 33 1 0 0 0 1 44 1 0 0 0 1 55

现在你可以使用排列向量来尝试一些排列p还有置换矩阵P．例如，语句: S (p)和P * S返回相同的矩阵。

: S (p)

Ans = 11 1 0 0 0 1 33 1 0 0 0 1 44 11 22 1 0 0 0 0 1 55

P * S

Ans = 11 1 0 0 0 1 33 1 0 0 0 1 44 11 22 1 0 0 0 0 1 55

同样的,S (: p)和* P '两种生产

* P '

Ans = 11 0 0 1 0 11 0 22 0 0 33 11 0 0 0 1 44 0 1 0 0 1 0 0 1 55

如果P是一个稀疏矩阵，那么两种表示使用的存储成正比n你可以申请其中之一年代在时间上正比于nnz (S)．矢量表示稍微紧凑且有效，因此各种稀疏矩阵置换程序都返回全行向量，除了LU（三角形）分解中的枢转置换，返回与完整LU分解兼容的矩阵。

在两种排列表示之间进行转换:

n = 5;I = speye (n);公关=我(p:);电脑=我(:p);电脑= (1:n) *电脑

PC = 1 3 4 2 5

公关= (pr * (1: n)“)”

Pr = 1 3 4 2 5

的倒数P仅仅是R = P '．你可以计算的倒数p与r (p) = 1: n．

r (p) = 1:5

1 .答案为c

重新排序稀疏性

对矩阵的列重新排序通常可以使其LU或QR因子更稀疏。对行和列重新排序通常可以使其Cholesky因子变得更稀疏。最简单的这种重新排序是按非零计数对列进行排序。对于结构非常不规则的矩阵，这有时是一种很好的重新排序方法，特别是在行或列的非零计数有很大变化的情况下。

的colperm计算一种排列，该排列将矩阵的列按每列中非零的数目从最小到最大进行排序。

重新排序以减少带宽

反向Cuthill-McKee顺序旨在减少矩阵的剖面或带宽。它不能保证找到最小的可能带宽，但它通常能做到。的symrcm函数实际上作用于对称矩阵的非零结构+“，但其结果对非对称矩阵也是有用的。这种排序对于来自一维问题或某种意义上的问题的矩阵是有用的细长．

近似最小度排序

图中一个节点的度是到该节点的连接数。这与邻接矩阵对应行中非对角非零元素的个数相同。近似最小度算法根据这些度在高斯消去或乔尔斯基分解过程中如何改变而生成一个排序。它是一个复杂而强大的算法，通常会导致比大多数其他排序(包括列计数和反向Cuthill-McKee)更稀疏的因素。由于跟踪每个节点的度非常耗时，所以近似最小度算法使用的是度的近似，而不是精确度。

这些MATLAB^®函数实现近似最小度算法:

symamd-用于对称矩阵。
colamd-使用非对称矩阵和对称矩阵的形式*“或“*．

看到稀疏矩阵的重新排序和分解举个例子symamd．

属性可以更改与算法细节相关的各种参数spparms函数。

使用的算法colamd和symamd,请参阅［5］．算法使用的近似程度是基于的[1]．

嵌套解剖订单

与近似最小次排序相似，嵌套的分解排序算法由解剖函数将矩阵视为图的邻接矩阵，从而对矩阵的行和列重新排序。该算法通过将图中的顶点对折叠在一起，将问题缩减到一个更小的规模。在对小图重新排序后，该算法然后应用投影和细化步骤将图扩展回原来的大小。

嵌套解剖算法产生高质量的重排，与其他重排技术相比，在有限元矩阵中表现尤其出色。有关嵌套解剖排序算法的更多信息，请参见［7］．

因子稀疏矩阵

LU分解
柯列斯基分解
QR分解
不完全分解

LU分解

如果年代是一个稀疏矩阵，下面的命令返回三个稀疏矩阵l，U,P这样p * s = l * u．

陆(L U P) =(年代);

陆用部分枢轴的高斯消去法得到因子。置换矩阵P只有n非零元素。和密集矩阵一样，这个表述陆[L U] = (S)返回一个排列单位的下三角矩阵和其乘积为的上三角矩阵年代．就其本身而言,陆(S)返回l和U在一个没有主元信息的矩阵中。

三个输出语法陆(L U P) = (S)选择P通过数值部分枢轴，但不枢轴来提高稀疏性陆的因素。另一方面，四输出语法[L U P, Q] =陆(S)选择P通过阈值部分枢轴，并进行选择P和问以改善稀疏性陆的因素。

你可以在稀疏矩阵中使用

陆(年代,打)

在哪里阈值是[0,1]中的主阈值。当一列中的对角线分量的大小小于时，就会发生旋转阈值乘以这一列中任何次对角线元素的大小。打= 0部队斜旋转。打= 1是默认的。(默认为阈值是0．1对于四输出语法）。

当你打电话陆有三个或更少的输出，MATLAB自动分配必要的内存来保存稀疏l和U分解过程中的因子。除了四输出语法，MATLAB没有使用任何符号逻辑单元预分解来确定内存需求和预先设置数据结构。

稀疏矩阵的重新排序和分解

打开生活的脚本

这个例子展示了重排序和分解对稀疏矩阵的影响。

如果你得到一个好的列排列p这减少了填充物，也许是从symrcm或colamd,然后计算陆(S (: p))比计算花费更少的时间和存储空间陆(S)．

使用Bucky球示例创建一个稀疏矩阵。

B =巴基;

B在每一行和每一列中恰好有三个非零元素。

创建两个排列,r和米使用symrcm和symamd分别。

r = symrcm (B);m = symamd (B);

这两种排列是对称反Cuthill-McKee排列和对称近似最小次排列。

创建间谍图来显示带有这三个不同数字的Bucky Ball图的三个邻接矩阵。当地的，以五角大楼为基础的原始编号结构在其他地方并不明显。

图subplot(1,3,1)“原始”)子情节(1,3,2)间谍(B(r,r))标题(“反向Cuthill-McKee”)子情节(1,3,3)间谍(B(m,m))标题(最小程度的）

与Cuthill-McKee的顺序相反，r，减少了带宽，并将所有非零元素集中在对角线附近。近似最小次排序，米，产生具有大块零块的形状状结构。

要查看在Bucky球的LU分解中生成的填充，请使用speye的对角线上插入-3B．

B = B - 3*speye(size(B));

因为每个行和现在是零，这个新的B实际上是奇异的，但计算其逻辑单元分解仍然具有指导意义。当只使用一个输出参数调用时，陆返回两个三角形因子，l和U，在一个稀疏矩阵中。矩阵中非零的数量是求解涉及的线性系统所需时间和存储的度量B．

以下是正在考虑的三种排列的非零计数。

陆(B)(原):1022
陆(B (r, r))（反向Cuthill-McKee）：968
陆(B (m m))(近似最低学位):636

尽管这只是一个小示例，但结果却是典型的。原始的编号方案导致最多的填充。反向Cuthill-McKee订单的填充集中在带内，但它几乎和前两个订单一样广泛。对于近似最小次排序，在消去过程中保留了相对较大的零块，且填充量明显小于其他排序生成的零块。

的间谍下面的图反映了每个重新排序的特征。

图子情节(1,3,1)间谍(陆(B))标题(“原始”)子情节(1,3,2)间谍(陆(乙(r,r)))“反向Cuthill-McKee”)子情节(1,3,3)间谍(陆(B(m,m)))最小程度的）

柯列斯基分解

如果年代是对称(或厄米)，正定，稀疏矩阵，下面的陈述返回一个稀疏，上三角矩阵R这R * R =年代．

R =胆固醇(S)

胆固醇不自动枢轴为稀疏性，但您可以计算近似最小度和轮廓限制排列使用胆固醇(S (p, p))．

由于Cholesky算法不使用旋转来表示稀疏性，也不需要旋转来表示数值稳定性，胆固醇快速计算所需的内存数量，并在分解开始时分配所有内存。您可以使用symbfact，它使用的算法与胆固醇，以计算分配了多少内存。

QR分解

MATLAB计算稀疏矩阵的完整QR分解年代与

(Q, R) = qr (S)

或

[Q, R E] = qr (S)

但这通常是不切实际的。酉矩阵问零元素的比例往往不高。还有一种更实用的替代方法，有时被称为“无q QR分解”。

有一个稀疏输入参数和一个输出参数

R = qr (S)

只返回QR分解的上三角形部分。矩阵R提供了与正规方程相关的矩阵的Cholesky分解:

‘* R = S *

然而，在计算中固有的数值信息的损失*年代年代是可以避免的。

使用两个具有相同行数的输入参数和两个输出参数，语句

[C R] = qr(年代,B)

应用正交变换到B、生产C =问' * B没有计算问．

Q-Lsol QR分解允许稀疏最小二乘问题的解决方案

$最小化 {为一个 x - b 为}_{2}$

有两个步骤:

[c R] = qr (A, b);x = R \ c

如果一个是稀疏的，但不是正方形的，MATLAB使用以下步骤来求解线性方程的反斜杠运算符:

x = A \ b

或者，你可以自己分解并检查R为等级不足。

也有可能解出一系列右边不同的最小二乘线性方程组，b这并不一定是在什么时候R = qr (A)计算。该方法解决了“半正规方程”‘* R * x =“* b与

x = R \ (R \(“* b))

然后采用一步迭代细化来减少舍入误差:

r = b - A*x;e = R \ (R \ (* R));X = X + e

不完全分解

的ilu和ichol函数提供近似,不完整的分解，它作为稀疏迭代方法的先决条件是有用的。

的ilu函数产生三不完整的右下方(ILU)分解:补零ILU（ILU (0))，是ILU的克劳特版本(ILUC(τ))， ILU阈值下降并发生旋转(ILUTP(τ)）.的ILU (0)永远不要枢轴，结果因子只有在输入矩阵非零的位置上才有非零。这两个ILUC(τ)和ILUTP(τ)，但是，使用用户定义的丢弃容忍度进行基于阈值的丢弃τ．

例如:

一个=画廊(“纽曼”， 1600) + speye(1600);

ans = 7840

陆nnz (())

ANS = 126478.

显示,一个有7840个非零，它的完全LU分解有126478个非零。另一方面，下面的代码显示了不同的ILU输出:

[L U] = ilu(一个);nnz (L) + nnz (U)造(A, 1)

ans = 7840

规范(A - (L * U)。* spones (A),“摇来摇去”)。/规范(,“摇来摇去”）

ans = 4.8874 e-17

选择。类型=“ilutp”；opts.dropol = 1e-4;[L,U,P] = ilu(A, opts);nnz (L) + nnz (U)造(A, 1)

ans = 31147

常模(P * A - L * U,“摇来摇去”)。/规范(,“摇来摇去”）

ans = 9.9224 e-05

选择。类型=“克劳特”；[L,U,P] = ilu(A, opts);nnz (L) + nnz (U)造(A, 1)

ans = 31083

常模(P * l * U,“摇来摇去”)。/规范(,“摇来摇去”）

ans = 9.7344 e-05

这些计算表明，零填充因子有7840个非零ILUTP(1)的军医因子有31147个非零ILUC(1)的军医因子有31083个非零。此外，零填充因子的乘积的相对误差在模式上基本上为零一个．最后，由于阈值下降而产生的因子分解中的相对误差与下降公差相同，尽管这并不能保证会发生。看到ilu更多选项和细节参考页。

的ichol函数提供了零填充不完整的尖头浓缩症状（集成电路(0))以及基于阈值的不完全choolesky分解(ICT(τ)对称的，正定的稀疏矩阵。这些分解是上面不完全逻辑单元分解的类似物，并且有许多相同的特征。例如:

一个= delsq (numgrid (“年代”, 200));nnz (A)

ans = 195228

nnz（辣椒（a，“低”））

ans = 7762589

显示,一个有195228个非零，而它的完全乔尔斯基分解没有重新排序有7762589个非零。相比之下:

L = ichol(一个);nnz(左)

ans = 117216

规范(A - (L * L ')。* spones (A),“摇来摇去”)。/规范(,“摇来摇去”）

ans = 3.5805 e-17

选择。类型=“信息与通信技术”；opts.dropol = 1e-4;l = iChol（A，OPTS）;nnz(左)

ans = 1166754

规范(L * L ',“摇来摇去”)。/规范(,“摇来摇去”）

ans = 2.3997 e-04

集成电路(0)非零只存在于下三角形的模式中一个的模式一个，各因子的乘积相匹配。此外,ICT(1)的军医因子比完整的Cholesky因子要稀疏得多，且两者之间的相对误差一个和L * L '在滴差公差上是相同的。需要注意的是，不同于胆固醇，提供的默认因素ichol下三角。看到ichol参考页获取更多信息。

线性方程组

解联立线性方程组有两种不同的方法:

直接的方法通常是高斯消去法的变体。这些方法通过矩阵运算(如LU或Cholesky因数分解)直接涉及单个矩阵元素。MATLAB通过矩阵除法运算符/和\实现了直接的方法，你可以用它来解决线性系统。
迭代的方法经过有限的步骤后只能得到一个近似解。这些方法只是通过矩阵向量乘积或抽象线性算子间接地涉及系数矩阵。迭代方法通常只应用于稀疏矩阵。

直接的方法

如果有足够的存储空间，直接方法通常比间接方法更快，更普遍地适用。迭代方法通常适用于方程的限制情况，并依赖于对角占优或基础微分算子的存在等性质。直接方法在MATLAB软件的核心中实现，并使一般矩阵类尽可能有效。迭代方法通常在matlab语言文件中实现，可以直接求解子问题或预处理。

使用不同的预定顺序。如果一个不是对角线，带状，三角形或三角形矩阵的置换，反斜杠（\）重新排序一个为了减少填充的数量，即加入稀疏分解矩阵的非零项的数量。新的排序，叫做预订，在因子分解之前执行一个．在某些情况下，您可能能够提供比反斜杠算法所使用的更好的预定顺序。

要使用不同的预排序，首先关闭两个自动预排序，默认情况下，使用函数反斜杠可能执行spparms如下：

defaultParms = spparms;spparms (“autoamd”, 0);spparms (“autommd”, 0);

现在，假设你已经创建了一个置换向量p的索引的排序一个，对矩阵应用反斜杠(: p)，它的列是的列一个，根据向量排列p．

x = A(:，p) \ b;x (p) = x;spparms (currentParms);

命令spparms (defaultParms)将控件恢复到以前的状态，以防您使用一个\ b，而不指定适当的预定顺序。

迭代的方法

有11个函数可以实现线性方程组的稀疏方程组的迭代方法。

稀疏系统迭代方法的函数

函数	方法
`bicg`	双共轭梯度
`bicgstab`	双共轭梯度稳定
`bicgstabl`	双共轭梯度稳定(l)
`研究生院理事会`	共轭梯度平方
`巨磁电阻`	广义最小剩余
`lsqr`	最小二乘
`minres`	最小剩余
`pcg`	预处理共轭梯度
`qmr`	Quasiminimal残留
`symmlq`	对称的LQ
`tfqmr`	Transpose-free quasiminimal残留

这些方法旨在解决一个x＝b或者最小化范数b- - - - - -一个x．对于预条件共轭梯度法，pcg，一个必须是对称的正定矩阵。minres和symmlq可用于对称不定矩阵。为lsqr，矩阵不需要是方阵。其他七种方法可以处理非对称方阵，每种方法都有独特的优点。

这11种方法都可以使用前置条件。线性系统

$一个 x ＝ b$

是否被等效系统所取代

$米^{- 1} 一个 x ＝米^{- 1} b$

的预调节器米，以加速迭代法的收敛。在许多情况下，预处理条件在数学模型中很自然地出现。例如，变系数偏微分方程可以用常系数近似。不完全矩阵分解可以在没有自然的预处理条件的情况下使用。

拉普拉斯方程在二维平面上的五点有限差分近似提供了一个例子。下面的陈述使用了预条件共轭梯度法的预条件米＝l＊l”,即l的零填充不完全Cholesky因子一个．

一个= delsq (numgrid (“年代”，50））;b =θ（大小（a，1），1）;tol = 1e-3;maxit = 100;L = ichol(一个);[x,国旗,犯错,iter, res] = pcg (A, b,托尔,麦克斯特,L, L ');

需要21次迭代才能达到规定的精度。另一方面，使用不同的预处理程序可能会产生更好的结果。例如,使用ichol为了构造一个修正的不完全Cholesky，只经过15次迭代就满足了规定的精度:

L = ichol(结构体(“类型”，“nofill”，“michol”，“上”））;[x,国旗,犯错,iter, res] = pcg (A, b,托尔,麦克斯特,L, L ');

有关这些迭代方法和不完全因子分解的背景信息可在［2］和［6］．

特征值和奇异值

两个函数可用来计算几个指定的特征值或奇异值。圣言会是基于eigs．

用于计算几个特征值或奇异值的函数

函数	描述
`eigs`	一些特征值
`圣言会`	一些奇异值

这些函数最常用于稀疏矩阵，但它们也可以用于完整矩阵，甚至是MATLAB代码中定义的线性算子。

该声明

[V,λ]= eigs (k,σ)

找到了k矩阵的特征值和相应的特征向量一个离" shift "最近的σ．如果σ，则求得幅度最大的特征值。如果σ为零时，特征值的大小最小。第二个矩阵,B，可以包含为广义特征值问题:Aυ = λBυ。

该声明

[U, V] =圣言(A、k)

找到了k最大的奇异值一个和

(U, V) =圣言(k,“最小”）

找到了k最小奇异值。

中使用的数值技术eigs和圣言会介绍了［6］．

稀疏矩阵的最小特征值

打开生活的脚本

这个例子说明了如何求稀疏矩阵的最小特征值和特征向量。

在一个65 × 65的网格上建立5点拉普拉斯差分算子l- 复杂的二维域。

L = numgrid (“L”, 65);一个= delsq(左);

确定非零元素的顺序和数量。

尺寸（a）

ans =1×22945 2945

nnz (A)

ans = 14473

一个是一个2945阶矩阵，包含14473个非零元素。

计算最小的特征值和特征向量。

[v d] = eigs (1,“smallestabs”）;

将特征向量的分量分配到适当的网格点上，并生成结果的等值线图。

L (L > 0) =全(v (L (L > 0)));x = 1:1/32:1;轮廓(x, x, L, 15)轴广场

参考

Amestoy, P. R.， T. A. Davis, and I. S. Duff，《近似最小度排序算法》，矩阵分析与应用，第17卷第4期，1996年10月，第886-905页。

[2] Barrett, R.， M. Berry, t.f. Chan, et al.，线性系统解的模板:迭代方法的构建块， SIAM，费城，1994。

[3] Davis, T.A.， Gilbert, J. R.， Larimore, S.I.， Ng, E.， Peyton, B.，“列近似最小次排序算法”，应用线性代数SIAM会议1997年10月，第29页。

[4] Gilbert, John R.， Cleve Moler，和Robert Schreiber，“MATLAB中的稀疏矩阵:设计和实现”，SIAM J. Matrix Anal。:.，第13卷第1期，1992年1月，第333-356页。

Larimore, S. I.，近似最小度列排序算法佛罗里达大学计算机与信息科学与工程系硕士论文，盖恩斯维尔，佛罗里达，1998。

[6]萨阿德,尤瑟夫稀疏线性方程的迭代方法．出版公司，1996。

Karypis, George和Vipin Kumar。一种快速、高质量的不规则图划分多级方案。SIAM科学计算杂志．卷。20，第1,199号，第359-392页。