LAPACKE C++ 实矩阵求逆答案

【问题标题】：LAPACKE C++ Real Matrix InversionLAPACKE C++ 实矩阵求逆
【发布时间】：2020-05-11 02:15:18
【问题描述】：

我正在尝试在 C++ LAPACKE 中反转实数矩阵。对于复杂的矩阵，我有相同的功能并且它有效。但真实案例给出了错误的答案。这是我的功能：

void inv(std::vector<std::vector<double>> &ans, std::vector<std::vector<double>> MAT){

    int N = MAT.size();

    int *IPIV = new int[N];

    double * arr = new double[N*N];
    for (int i = 0; i<N; i++){
        for (int j = 0; j<N; j++){
            int idx = i*N + j;
            arr[idx] = MAT[i][j];
        }
    }

    LAPACKE_dgetrf(LAPACK_ROW_MAJOR, N, N, arr, N, IPIV);
    LAPACKE_dgetri(LAPACK_ROW_MAJOR, N, arr, N, IPIV);

     for (int i = 0; i<N; i++){
        for (int j = 0; j<N; j++){
            int idx = i*N + j;
            ans[i][j] = arr[idx];
        }
    }

    delete[] IPIV;
    delete[] arr;
}

我尝试反转一个 24 x 24 的双精度矩阵。虽然程序似乎几乎就在那里，但逆还没有完全出现，它与 python linalg inverse 给我的有很大不同（python 就在这里，因为我将矩阵乘以逆，结果非常接近缩进）。在 LAPACKE 输出中，我将矩阵乘以它的逆矩阵，我得到对角线为 1，但非对角线的值高达 0.17，与 0 相比，这是巨大的。有没有办法让 LAPACKE 程序提供更好的结果？谢谢！

【问题讨论】：

是定义好的矩阵吗？矩阵行列式的值是多少？
2.2864066779666567e+37。它的值很大，所以我认为这就是程序遇到问题的原因？
是的，这可能解释了效果。您可以尝试将输入矩阵（每个项）乘以（1/25），这将使行列式减少 25^{-24}=2.8e-34 的因子（如果我没记错的话），从而使输入矩阵行列式大约 1000。然后计算逆并乘以 1/25。我在答案中放了简单的 python 代码
查看 numpy 的 linalg.inv() 的来源，特别是 umath_linalg.c 的第 1555 - 1693 -1727 行，它解决了 AX=Id，使用 dgesv() 而不是调用 dgetrf() 和dgetri()。 Scipy 使用 dgetrf() / dgetri()。您可以尝试在 python 中使用 scipy 的 scipy.linalg.inv() 并查看它是否也会导致正确的输出？

标签： c++ matrix lapack matrix-inverse lapacke

【解决方案1】：

对于具有大行列式的矩阵，您可以重新调整输入，计算逆，然后重新调整输出。这是一个非常简单的 Python 示例，您的比例因子应该是 1/25 左右才能获得 (1/25)²⁴=2.8e-34 的总乘数，从而使输入矩阵行列式约为 1000 .

import numpy as np

scale = 0.5

i = np.array([[1,2],[3,4]])
print(i)
print(np.linalg.det(i))
print("-----------------------------------")
x = np.multiply(i, [scale]) # rescale matrix
print(x)
print(np.linalg.det(x))     # determinant should be less
print("-----------------------------------")
y = np.linalg.inv(x)
print(y)
print(np.linalg.det(y))
print("-----------------------------------")
o = np.multiply(y, [scale]) # rescale matrix
print(o)
print(np.linalg.det(o))
print(np.dot(i, o))

您可以使用scale 并验证代码的任何值是否返回正确的反转输入矩阵。

所以你的代码是

double scale = 1.0/25.0;

for (int i = 0; i<N; i++){
    for (int j = 0; j<N; j++){
        int idx = i*N + j;
        arr[idx] = scale*MAT[i][j];
    }
}

....

for (int i = 0; i<N; i++){
    for (int j = 0; j<N; j++){
        int idx = i*N + j;
        ans[i][j] = scale * arr[idx];
    }
}

【讨论】：