使用 Matlab 的反斜杠运算符来反转稀疏矩阵会导致某些条目向下舍入为零

问题描述

我正在使用 Matlab 的 backslash 运算符来求解写为两个矩阵 M1 和 M2 的方程组。这两个矩阵是方阵和 tridiagonal，因此我将它们定义为稀疏矩阵。例如，每个条目的尺寸为 5x5，它们的定义如下，每个条目中的值取决于某个常量 a：

N = 5;
a =  1e10;
M1 = spdiags([-a*ones(N,1)...        % Sub diagonal
              (1 + 2*a)*ones(N,1)... % Main Diagonal
              -a*ones(N,1)],...      % Super diagonal
              -1:1,N,N);
M2 = spdiags([+a*ones(N,1)...
              (1 - 2*a)*ones(N,1)...
              +a*ones(N,...
              -1:1,N);
M_out = M1\M2;

例如，M1 的完整形式如下所示：

>> full(M1)

ans =

   1.0e+10 *

    2.0000   -1.0000         0         0         0
   -1.0000    2.0000   -1.0000         0         0
         0   -1.0000    2.0000   -1.0000         0
         0         0   -1.0000    2.0000   -1.0000
         0         0         0   -1.0000    2.0000

现在，如果我检查结果 M_out 中非零条目的数量，那么我可以看到它们都是非零的，这很好：

>> nnz(M_out)

ans =

    25

问题是我还需要为常量 a 的较大值执行此操作。但是，例如，如果改为 a=1e16，那么 M_out 的非对角线条目会自动设置为零，大概是因为它们变得太小了：

>> nnz(M_out)

ans =

    5

在 Matlab 中是否有更好的方法来解决稀疏矩阵求逆的问题？还是我以错误的方式使用了反斜杠运算符？

解决方法

如果矩阵的大小不会增长太多，我建议进行完整的符号计算：

N = 5;
syms a
M1 = diag(-a*ones(N-1,1),-1) + diag((1 + 2*a)*ones(N,0) + diag(-a*ones(N-1,+1);
M2 = diag(+a*ones(N-1,-1) + diag((1 - 2*a)*ones(N,0) + diag(+a*ones(N-1,+1);
M_out = M1\M2;

M_num_1e10 = subs(M_out,a,1e10);
M_num_1e16 = subs(M_out,1e16);

vpa(M_num_1e10)
vpa(M_num_1e16)

在这种情况下，您将需要 Symbolic Math Toolbox。如果您没有它，我认为您应该考虑迁移到 Python 并使用 SymPy。

编辑：

考虑到您定义问题的方式，您的计算需要更高的精度。双精度是不够的。例如，在双精度中 (1e16+1) 必须四舍五入为 (1e16)，也就是说 (1e16+1)-(1e16) 等于零。所以你的问题始于矩阵的主对角线。 MATLAB 仅通过其符号工具箱提供扩展精度。

如果您想坚持使用双精度，您可以依靠所谓的双双算术自行扩展双精度。我说你必须自己做，因为我认为MATLAB没有开源的double-double库。