C++ - 计算矩阵的逆

Question

我试图编写一个能够计算矩阵的逆的程序：

以下是目前为止我已经写好的代码：

#include <iostream>
#include <vector>
#include <math.h>
#include <iomanip>
#include <stdexcept>

double getDeterminant(const std::vector<std::vector<double>> vect) {
    if(vect.size() != vect[0].size()) {
        throw std::runtime_error("Matrix is not quadratic");
    } 
    int dimension = vect.size();

    if(dimension == 0) {
        return 1;
    }

    if(dimension == 1) {
        return vect[0][0];
    }

    //Formula for 2x2-matrix
    if(dimension == 2) {
        return vect[0][0] * vect[1][1] - vect[0][1] * vect[1][0];
    }

    double result = 0;
    int sign = 1;
    for(int i = 0; i < dimension; i++) {

        //Submatrix
        std::vector<std::vector<double>> subVect(dimension - 1, std::vector<double> (dimension - 1));
        for(int m = 1; m < dimension; m++) {
            int z = 0;
            for(int n = 0; n < dimension; n++) {
                if(n != i) {
                    subVect[m-1][z] = vect[m][n];
                    z++;
                }
            }
        }

        //recursive call
        result = result + sign * vect[0][i] * getDeterminant(subVect);
        sign = -sign;
    }

    return result;
}

std::vector<std::vector<double>> getTranspose(const std::vector<std::vector<double>> matrix1) {

    //Transpose-matrix: height = width(matrix), width = height(matrix)
    std::vector<std::vector<double>> solution(matrix1[0].size(), std::vector<double> (matrix1.size()));

    //Filling solution-matrix
    for(size_t i = 0; i < matrix1.size(); i++) {
        for(size_t j = 0; j < matrix1[0].size(); j++) {
            solution[j][i] = matrix1[i][j];
        }
    }
    return solution;
}

std::vector<std::vector<double>> getCofactor(const std::vector<std::vector<double>> vect) {
    if(vect.size() != vect[0].size()) {
        throw std::runtime_error("Matrix is not quadratic");
    } 

    std::vector<std::vector<double>> solution(vect.size(), std::vector<double> (vect.size()));
    std::vector<std::vector<double>> subVect(vect.size() - 1, std::vector<double> (vect.size() - 1));

    for(std::size_t i = 0; i < vect.size(); i++) {
        for(std::size_t j = 0; j < vect[0].size(); j++) {

            int p = 0;
            for(size_t x = 0; x < vect.size(); x++) {
                if(x == i) {
                    continue;
                }
                int q = 0;

                for(size_t y = 0; y < vect.size(); y++) {
                    if(y == j) {
                        continue;
                    }

                    subVect[p][q] = vect[x][y];
                    q++;
                }
                p++;
            }
            solution[i][j] = pow(-1, i + j) * getDeterminant(subVect);
        }
    }
    return solution;
}

std::vector<std::vector<double>> getInverse(const std::vector<std::vector<double>> vect) {
    if(getDeterminant(vect) == 0) {
        throw std::runtime_error("Determinant is 0");
    } 
    double d = 1.0/getDeterminant(vect);
    std::vector<std::vector<double>> solution(vect.size(), std::vector<double> (vect.size()));

    for(size_t i = 0; i < vect.size(); i++) {
        for(size_t j = 0; j < vect.size(); j++) {
            solution[i][j] = vect[i][j] * d; 
        }
    }

    return getTranspose(getCofactor(solution));
}

void printMatrix(const std::vector<std::vector<double>> vect) {
    for(std::size_t i = 0; i < vect.size(); i++) {
        for(std::size_t j = 0; j < vect[0].size(); j++) {
            std::cout << std::setw(8) << vect[i][j] << " ";
        }
        std::cout << "\n";
    }
}

int main() {
    std::vector<std::vector<double>> matrix(3, std::vector<double> (3));
    matrix = {
        {1,2,3},
        {4,5,6},
        {7,8,8}
    };

    printMatrix(getInverse(matrix));
    return 0;
}

计算行列式、转置矩阵和余子式矩阵的函数应该没有问题（就我所看到的），但计算逆矩阵的函数不正确。 我在互联网上搜索并找到了这个，它使用相同的函数来计算逆矩阵。

这个公式是否不正确，或者您是否有任何其他想法，为什么它不起作用？

我正在使用的矩阵是


它的逆矩阵应该是


- user12493956

---

你尝试过在代码执行的时候使用调试器逐步执行吗？这通常是缩小问题范围的好方法。在这种情况下，传递一个小矩阵，手动（用笔和纸）求逆并查看代码是否跟踪同样的计算。 - Cory Kramer

这个公式似乎部分正确，因为它对一些示例矩阵（如[[1,2],[3,4]]）有效。你的代码中可能存在一些更微妙的错误。你需要逐步检查并找出第一次出现不正确数字的位置。 - Bartek Banachewicz

我不记得如何计算逆矩阵（所有类型的通用公式）。也许在https://math.stackexchange.com/上，你可以找到一些专家来帮助你解决算法问题。 - Leonardo Alves Machado

另外：如果您要使用std::vector<std::vector<double>>来表示矩阵，您应该检查每个内部向量的大小是否与外部向量相同。 if(std::any_of(vect.begin(), vect.end(),[size = vect.size()](auto & inner){ return inner.size() != size; })) { throw ... } 这也将消除dimension == 0情况下的未定义行为。 - Caleth

您可能不会回来阅读评论，但大多数系统不会以这种方式计算矩阵的逆。更常见的方法是使用带有枢轴的LU分解。有一本名为“C语言数字计算法”的好书提供了很好的介绍和代码。 - duffymo

请注意，您实现的算法非常慢，因为您必须计算所有这些子行列式。我认为它的规模是O(n!)。更好的方法是使用高斯消元算法https://en.wikipedia.org/wiki/Gaussian_elimination，其时间复杂度为O(n^3)。不过，代码还是很好的。 - tomtom1-4

Answer 1

首先，感谢您的评论。

问题在于执行顺序。 正确的解决方案是：

std::vector<std::vector<double>> getInverse(const std::vector<std::vector<double>> vect) {
    if(getDeterminant(vect) == 0) {
        throw std::runtime_error("Determinant is 0");
    } 

    double d = 1.0/getDeterminant(vect);
    std::vector<std::vector<double>> solution(vect.size(), std::vector<double> (vect.size()));

    for(size_t i = 0; i < vect.size(); i++) {
        for(size_t j = 0; j < vect.size(); j++) {
            solution[i][j] = vect[i][j]; 
        }
    }

    solution = getTranspose(getCofactor(solution));

    for(size_t i = 0; i < vect.size(); i++) {
        for(size_t j = 0; j < vect.size(); j++) {
            solution[i][j] *= d;
        }
    }

    return solution;
}

- user12493956

Answer 2

我知道这是一个旧问题，但当输入矩阵的维度为1时，您的代码无法工作。 这是我使用的解决方法：

vector<vector<double>> inverse(const vector<vector<double>> A) {
    double d = 1.0/det(A);
    vector<vector<double>> solution(A.size(), vector<double> (A.size()));

    if(A.size() == 1){
        vector<double> ans = {0};
        ans[0] = 1.0/det(A);
        solution[0] = ans;
        return solution;
    }

    for(size_t i = 0; i < A.size(); i++) {
        for(size_t j = 0; j < A.size(); j++) {
            solution[i][j] = A[i][j] * d; 
        }
    }

    return transpose(cofactor(solution));
}

- Jaden Park