如何计算float类型的大矩阵的逆？

Question

我用高斯消元法对矩阵求逆，矩阵是12 x 12。

当求小矩阵时结果是正确的，但当矩阵变大时结果不正确。

我猜这是浮点小数精度问题。

有没有办法使用单精度浮点小数计算精确的逆矩阵？

template<typename T = float>
inline bool inverse(const MATRIX<Mtx_Rt, Mtx_Rt, T>& matrix, MATRIX<Mtx_Rt, Mtx_Rt, T>& out, void* pTmp) {

       if (matrix.GetNumOfRow() != matrix.GetNumOfCol())
        return false;

    MATRIX<Mtx_Rt, Mtx_Rt, T> augmentedMatrix = MATRIX<Mtx_Rt, Mtx_Rt, T>(matrix.GetNumOfRow(), matrix.GetNumOfCol() * 2, (T*)pTmp);
    memset(pTmp, 0, augmentedMatrix._sizeof());

    //[matrix | identity]
    for (int i = 0; i < matrix.GetNumOfRow(); i++) {
        for (int j = 0; j < matrix.GetNumOfCol(); j++)
            augmentedMatrix(i, j) = matrix(i, j);

        augmentedMatrix(i, i + matrix.GetNumOfRow()) = 1.0;
    }

    for (int i = 0; i < matrix.GetNumOfRow(); i++) {
        int pivotRow = i;
        for (int j = i + 1; j < matrix.GetNumOfRow(); j++)
            if (augmentedMatrix(j, i) > augmentedMatrix(pivotRow, i))
                pivotRow = j;

        for (int k = 0; k < 2 * matrix.GetNumOfRow(); k++) {
            T tmp = augmentedMatrix(i, k);
            augmentedMatrix(i, k) = augmentedMatrix(pivotRow, k);
            augmentedMatrix(pivotRow, k) = tmp;
        }

        for (int j = 0; j < matrix.GetNumOfRow(); j++)
            if (j != i) {
                T ratio;
                if (augmentedMatrix(i, i) != 0.0)
                    ratio = augmentedMatrix(j, i) / augmentedMatrix(i, i);
                else
                    return false;

                for (int k = 0; k < 2 * matrix.GetNumOfRow(); k++)
                    augmentedMatrix(j, k) -= ratio * augmentedMatrix(i, k);
            }
    }

    for (int i = 0; i < matrix.GetNumOfRow(); i++) {
        T diagonalElement = augmentedMatrix(i, i);
        for (int j = 0; j < 2 * matrix.GetNumOfRow(); j++)
            augmentedMatrix(i, j) /= diagonalElement;

    }

    for (int i = 0; i < matrix.GetNumOfRow(); i++)
        for (int j = 0; j < matrix.GetNumOfCol(); j++)
            out(i, j) = augmentedMatrix(i, j + matrix.GetNumOfRow());

    return true;
}

//The following is the test data
0x000001DEBFA44B70     0.00901737064    0.00851101335    0.00903018098    0.00852382369   1.71711508e-06   1.07878850e-05   1.07854212e-05  -5.71605915e-06       0.00000000       0.00000000       0.00000000       0.00000000  
0x000001DEBFA44BA0     0.00851101428    0.00921893492    0.00835073180    0.00905865245   0.000142186589   6.97691357e-05   6.96772186e-05   0.000119451797       0.00000000       0.00000000       0.00000000       0.00000000  
0x000001DEBFA44BD0     0.00903018005    0.00835073087    0.00916963257    0.00849018339  -0.000121859746  -5.83902292e-05  -5.83116416e-05  -0.000103386286       0.00000000       0.00000000       0.00000000       0.00000000  
0x000001DEBFA44C00     0.00852382462    0.00905865245    0.00849018339    0.00902501121   1.86097168e-05   5.91026037e-07   5.80158485e-07   2.17815614e-05       0.00000000       0.00000000       0.00000000       0.00000000  
0x000001DEBFA44C30    1.71710417e-06   0.000142186575  -0.000121859761   1.86097132e-05      0.141640529       0.00000000       0.00000000       0.00000000    0.00899084471    0.00772444298   -0.00682374742   -0.00809014961  
0x000001DEBFA44C60    1.07878877e-05   6.97691430e-05  -5.83902220e-05   5.91027856e-07       0.00000000      0.122075945      0.122108623     0.0703755319     0.0109803220    0.00996958558   -0.00884070992   -0.00985144638  
0x000001DEBFA44C90    1.07854303e-05   6.96772186e-05  -5.83116343e-05   5.80155756e-07       0.00000000      0.122108623      0.122141339     0.0704018399     0.0109902890    0.00998063106   -0.00885062385   -0.00986028183  
0x000001DEBFA44CC0   -5.71606142e-06   0.000119451797  -0.000103386286   2.17815686e-05       0.00000000     0.0703755319     0.0704018399      0.110758379    0.00697547942    0.00992908329   -0.00901818927   -0.00606458541  
0x000001DEBFA44CF0        0.00000000       0.00000000       0.00000000       0.00000000    0.00899084471     0.0109803220     0.0109902890    0.00697547942      0.157142207      0.148548171      0.145359069      0.136765048  
0x000001DEBFA44D20        0.00000000       0.00000000       0.00000000       0.00000000    0.00772444159    0.00996958464    0.00998063106    0.00992908329      0.148548156      0.157674491      0.136261418      0.145387754  
0x000001DEBFA44D50        0.00000000       0.00000000       0.00000000       0.00000000   -0.00682374695   -0.00884070992   -0.00885062385   -0.00901818834      0.145359069      0.136261433      0.156348825      0.147251189  
0x000001DEBFA44D80        0.00000000       0.00000000       0.00000000       0.00000000   -0.00809014961   -0.00985144544   -0.00986028090   -0.00606458541      0.136765033      0.145387754      0.147251174      0.155873895  

Answer 1

使用模板数值工具包 [TNT] https://math.nist.gov/tnt/ 只需要很少的包含（无需构建库）和行代码（如下所示）即可使用 LU 分解计算逆矩阵。 它是一个模板类，因此您可以将它与 float 一起使用。我使用它来处理大小高达 2000 的矩阵已经有 10 多年了。 TNT 还包括 QR、特征值求解器和 SVD 分解。

template <class T>
TNT::Array2D<T> inverse(const TNT::Array2D<T>& M)
{
    TNT::Array2D<T> LU(M);
    int nline = dim1();
    int ncol = dim2();
    TNT::Array2D<T> input(nline, ncol);
    for (int i = 0; i < nline; i++)
    {
        for (int j = 0; j < ncol; j++)
        {
            if (i == j) input[i][j] = 1.;
            else input[i][j] = 0.;
        }
    }
    return LU.solve(input);
}