标签: matrix-inverse

我在 Matlab 中有一个大小为 10800x10800 的矩阵，我直接使用以下函数计算其逆矩阵：

inv(A)

仅计算一次就需要 3 至 4 分钟。这是迭代算法的一部分，需要超过 20 次迭代，所以总体来说速度会非常慢。有一个更好的方法吗？也许是一些数学公式或者更好的 Matlab 函数？

编辑：矩阵是对角的。每次迭代，对角线元素都会根据拟合因子分析器的公式进行更新。但这并不重要，重要的是它是一个对角矩阵，并且每次迭代都会改变。

谢谢

我正在尝试计算矩阵的逆P,但如果我乘以inv(P)*P,MATLAB不会返回单位矩阵.它几乎是身份(非对角线值的顺序10^(-12)).但是,在我的应用程序中我需要更高的精度.

在这种情况下我该怎么办？

所以我试图反转一个大的（449x449）协方差矩阵，因此它是对称且正定的。（我想做的是将这个矩阵反转，作为适合 Mauna Loa CO2 数据集的高斯过程的一部分。）

这个反演相当长，所以我想使用 chol2inv 而不是求解。但是 chol2inv 方法给了我一个非常奇怪的结果：一个非常接近 0 的矩阵（它的总和等于 10^(-13)）。

为什么 chol2inv 会给我这个？

我是matlab的新手,我试图找出具有非常小值的矩阵的逆.当我试图找到逆时,我得到一个错误,说矩阵是单数.建议的解决方案之一是尝试向对角线元素添加一些元素.我知道我必须使用眼睛和诊断方法,但我无法找到正确的解决方案.

任何评论都会有所帮助.

我有一个矩阵M,它有576行和15列.第一列是全部1.第二列从1开始,逐一增加到576.第三列是第二列的平方.第四列是第二列的第三列.第五列是第二列的第四列.它是这样的.最后第十五列是第二列的第十四次幂.操作是将矩阵M的转置与其自身相乘并取结果的倒数.然而,Matlab会发出这样的警告; 警告:矩阵接近单一或严重缩放.结果可能不准确.RCOND = 1.889839e-79.这是一种反过来的方式吗？我也试过这个操作符\但是发生了同样的警告.

我是使用 LAPACK/BLAS 的新手。我想计算一个方程的解：

AU=F

我想知道这部分代码的逻辑错误是什么。我使用求解器输入大小为 ((xdiv-1) (ydiv-1),(xdiv-1) (ydiv-1)) 的矩阵 A。然后，随后求解方程。U=逆（A）* f。

其中 U 和 f 大小相同。(u((xdiv-1) (ydiv-1),1),f((xdiv-1) (ydiv-1),1))。执行矩阵求逆时出现分段错误错误。

这是我的代码：

program main
double precision, allocatable :: A(:,:)
double precision, allocatable :: u(:,:), f(:,:)
double precision mesh(2), dx, dy
integer  xdiv, ydiv
 xdiv=55
 ydiv=55
 mesh(1)=.001
 mesh(2)=.001
 dx=mesh(1)
 dy=mesh(2)

allocate (A((xdiv-1)*(ydiv-1),(xdiv-1)*(ydiv-1))) 
allocate (Ainv((xdiv-1)*(ydiv-1),(xdiv-1)*(ydiv-1)))
allocate (u((xdiv-1)*(ydiv-1),1),f((xdiv-1)*(ydiv-1),1))

         do i =1,(xdiv-1)*(ydiv-1)
         A(i,i)=-2.d0*(1.d0/(dx**2)+1.d0/(dy**2))
         enddo
         do i=1,(xdiv-2)
         do j=1,(ydiv-1)
         A(i+(j-1)*(xdiv-1),i+(j-1)*(xdiv-1)+1)=1.d0/(dx**2)
         A(i+(j-1)*(xdiv-1)+1,i+(j-1)*(xdiv-1))=1.d0/(dx**2)
         enddo
         enddo
         do i=1,(xdiv-1)
         do j=1,(ydiv-2)
         A(i+(j-1)*(xdiv-1),i+(j)*(xdiv-1))=1.d0/(dy**2)
         A(i+(j)*(xdiv-1),i+(j-1)*(xdiv))=1.d0/(dy**2)
         enddo
         enddo

        do i=1,(xdiv-1)
         do j=1,(ydiv-1)
       xcoord = (i-1)*mesh(1) …

我目前正在 Python 上使用最小二乘算法，涉及一些大地测量计算。

我选择了Python（它不是最快的）并且它运行得很好。然而，在我的代码中，我有大型稀疏对称（非正定，因此不能使用 Cholesky）矩阵的逆来执行（下图）。我目前使用 np.linalg.inv() ，它使用 LU 分解方法。

我很确定在速度方面需要进行一些优化。

我想到了 Cuthill-McKee 算法来重新排列矩阵并取其逆。您有什么想法或建议吗？

非常感谢您的回答！

我的代码执行4x4矩阵求反（128、256、512）数时，遇到了精度问题。当我使用原始版本（即numpy函数np.linalg.inv或）时np.linalg.pinv，一切正常。

不幸的是，使用下面的CUDA代码，我将nan和inf值转换为倒置矩阵。

更明确地说，我将此矩阵求反：

2.120771107884677649e+09 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
0.000000000000000000e+00 3.557266600921528288e+27 3.557266600921528041e+07 3.557266600921528320e+17
0.000000000000000000e+00 3.557266600921528041e+07 3.557266600921528288e+27 3.557266600921528041e+07
0.000000000000000000e+00 3.557266600921528320e+17 3.557266600921528041e+07 1.778633300460764144e+27

如果使用经典的numpy“ inv”，则会得到以下3x3倒置矩阵：

4.715266047722758306e-10 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
0.000000000000000000e+00 2.811147187396482366e-28 -2.811147186834252285e-48 -5.622294374792964645e-38
0.000000000000000000e+00 -2.811147186834252285e-48 2.811147187396482366e-28 -5.622294374230735768e-48
0.000000000000000000e+00 -5.622294374792964645e-38 -5.622294374230735768e-48 5.622294374792964732e-28

为了检查该逆矩阵的有效性，我将其乘以原始矩阵，结果是单位矩阵。

但是使用CUDA GPU反转后，我得到了以下矩阵：

0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
-inf -inf -9.373764907941219970e-01 -inf
inf nan -inf nan

所以，我想竟被以提高精度到我的CUDA内核或Python代码，以避免这些nan和inf值。

这是CUDA内核代码，并调用了我的主要代码的一部分（我已经用numpy inv函数注释了经典方法：

    # Create arrayFullCross_vec array
    arrayFullCross_vec = np.zeros((dimBlocks,dimBlocks,integ_prec,integ_prec)) …

我编写了两个程序来使用高斯消元计算矩阵的逆，第一个程序是在 C# 中，第二个是在 CUDA C++ 中。这两个程序遵循完全相同的程序并给出相同的最终结果。但是，当我检查中间步骤中的值时，我发现值略有不同，小于 1e-5 的相对误差。

这是两个程序的每个代码的一部分。

C＃

int i, j, i1, n, y, z;
double[,] M = new double[n, n];
double[,] inv = new double[n, n];
for (i = 0; i < n; i++)
    inv[i, i] = 1;
for (i = 0; i < n; i++)
{
    for (j = i + 1; j < n; j++)
        M[i, j] /= M[i, i];
    for (j = 0; j < n; j++)
        inv[i, j] /= M[i, i]; …