matlab - 求解 *sparse* 上三角系统

Question

如果我想解决一个完整的上三角系统，我可以调用linsolve(A,b,'UT'). 但是，稀疏矩阵目前不支持此功能。我该如何克服呢？

Answer 1

UT 和 LT 系统是最容易解决的系统之一。在 wiki 上阅读有关它的信息。知道了这一点，就很容易编写自己的 UT 或 LT 求解器：

%# some example data
A = sparse( triu(rand(100)) );
b = rand(100,1);

%# solve UT system by back substitution    
x = zeros(size(b));
for n = size(A,1):-1:1    
    x(n) = (b(n) - A(n,n+1:end)*x(n+1:end) ) / A(n,n);    
end

LT 系统的过程非常相似。

话虽如此，使用 Matlab 的反斜杠运算符通常更容易和更快：

x = A\b

正如 nate 已经指出的，这也适用于备件矩阵。

请注意，此运算符还解决了具有非正方形A或如果在主对角线上A具有等于零（或 ）的某些元素的 UT 系统。< eps它以最小二乘的方式解决了这些情况，这可能对您来说是可取的，也可能不是您想要的。您可以在执行求解之前检查这些情况：

if size(A,1)==size(A,2) && all(abs(diag(A)) > eps)
    x = A\b;
else
    %# error, warning, whatever you want
end

通过键入阅读有关（反）斜杠运算符的更多信息

>> help \

或者

>> help slash

在 Matlab 命令提示符下。

Answer 2

编辑由于您需要的是三角求解过程，也称为后向/前向替换，因此您可以使用普通的 MATLAB 反斜杠\运算符：

x = U\b

如原始答案中所述，MATLAB 将识别出您的矩阵是三角形的事实。为了确定这一点，您可以将性能与SuiteSparsecs_usolve中的过程进行比较。它是一个用 C 语言实现的 mex 函数，用于计算上三角稀疏矩阵的稀疏三角求解（那里也有类似的函数：和）。cs_lsolvecs_utsolvecs_ltsolve

您可以查看原生 MATLAB 和稀疏 Cholesky 分解上下文中的性能比较。cs_l(t)solve从本质上讲，MATLAB 性能很好。唯一的陷阱是如果你想解决一个转置系统

x = U'\b

MATLAB无法识别并显式创建U. 在这种情况下，您应该cs_utsolve明确调用。

原始答案如果您的系统是对称的并且您只存储上三角矩阵部分（这就是我在您的问题中完全理解的方式），并且如果 Cholesky 分解适合您，如果您的矩阵是正定的，则​​ chol处理对称矩阵。对于不定矩阵，您可以使用ldl。两者都处理稀疏存储并处理对称矩阵部分。

较新的 matlab 版本为此使用cholmod 和 suitesparse。这是迄今为止我所知道的表现最好的 Cholesky 分解。在 matlab 中，它也使用并行 BALS 进行并行化。

您从上述函数中获得的因子是上三角矩阵 L 使得

A=LL'

您现在需要做的就是执行前向和后向替换，这既简单又便宜。在 matlab 中，这是在反斜杠运算符中自动完成的

x=L'\(L\b)

矩阵可以是稀疏的，matlab 会识别出它是上三角/下三角。您还可以将此调用与前向替换一起用于使用 cholesky 分解获得的因子。

Answer 3

您可以在稀疏矩阵上使用 MLDIVIDE( \ ) 或 MRDIVIDE( / ) 运算符...