Eigen与Boost矩阵求逆结果不一致问题求助

问题背景 #

开发C++矩阵求逆代码时，发现Eigen3（3.3.9版本）与Boost uBLAS（1.73版本）对18x18矩阵求逆的结果存在显著差异，差异矩阵的无穷范数约为130，直接导致最终算法输出完全不同。已排查相关技术线程未找到有效方案，希望实现不依赖特定求逆库的鲁棒解法。

测试矩阵构造逻辑：矩阵行由幂函数及其导数构成，完整测试代码如下：

#include <Eigen/Dense>

#include <boost/numeric/ublas/lu.hpp>
using namespace boost::numeric;

#include <cmath>
#include <stdexcept>
#include <iostream>
#include <iomanip>

#define SIZE 18

Eigen::MatrixXd initMatrix_Eigen() {
    Eigen::MatrixXd matrix(SIZE, SIZE);
    for (size_t i = 0; i < SIZE / 2; i++) {
        for (size_t j = 0; j < SIZE; j++) {
            matrix(2 * i, j) = pow((double)i + 1, (double)j + 2);
            matrix(2 * i + 1, j) = (j + 2) * pow((double)i + 1, (double)j + 1);
        }
    }
    return matrix;
}

ublas::matrix<double> initMatrix_Boost() {
    ublas::matrix<double> matrix(SIZE, SIZE);
    for (size_t i = 0; i < SIZE / 2; i++) {
        for (size_t j = 0; j < SIZE; j++) {
            matrix(2 * i, j) = pow((double)i + 1, (double)j + 2);
            matrix(2 * i + 1, j) = (j + 2) * pow((double)i + 1, (double)j + 1);
        }
    }
    return matrix;
}

Eigen::MatrixXd invertMatrix_Eigen(const Eigen::MatrixXd& matrix) {
    return matrix.inverse();
}

ublas::matrix<double> invertMatrix_Boost(const ublas::matrix<double>& matrix) {
    if (matrix.size1() != matrix.size2()) throw("Input matrix must be square");
    ublas::matrix<double> tmpMatrix(matrix);
    ublas::permutation_matrix<size_t> permutationMatrix(tmpMatrix.size1());
    if (ublas::lu_factorize(tmpMatrix, permutationMatrix) != 0) throw std::invalid_argument("LU factorization failed");
    ublas::matrix<double> inverse(tmpMatrix.size1(), tmpMatrix.size2());
    inverse.assign(boost::numeric::ublas::identity_matrix<double>(tmpMatrix.size1()));
    ublas::lu_substitute(tmpMatrix, permutationMatrix, inverse);
    return inverse;
}

int main(int argc, char** argv)
{
    std::cout << std::setprecision(12);

    Eigen::MatrixXd matrix_Eigen = initMatrix_Eigen();
    Eigen::MatrixXd inverse_Eigen = invertMatrix_Eigen(matrix_Eigen);
    std::cout << "Eigen inverse matrix :" << std::endl;
    std::cout << inverse_Eigen << std::endl;

    ublas::matrix<double> matrix_Boost = initMatrix_Boost();
    ublas::matrix<double> inverse_Boost_tmp = invertMatrix_Boost(matrix_Boost);
    // Convert inverse_Boost_tmp into an Eigen matrix called inverse_Boost :
    const size_t size = inverse_Boost_tmp.size1();
    Eigen::MatrixXd inverse_Boost(size, size);
    for (size_t i = 0; i < size; i++)
        for (size_t j = 0; j < size; j++)
            inverse_Boost(i, j) = inverse_Boost_tmp(i, j);
    std::cout << "Boost inverse matrix :" << std::endl;
    std::cout << inverse_Boost << std::endl;

    Eigen::MatrixXd difference = inverse_Eigen - inverse_Boost;
    std::cout << "Difference between inverse matrices :" << std::endl;
    std::cout << difference << std::endl;
    std::cout << "NormInf(difference) = " << difference.lpNorm<Eigen::Infinity>() << std::endl;
    return 0;
}

核心原因分析 #

1. 矩阵条件数过大：你构造的矩阵属于Vandermonde类矩阵（行由幂函数构成），这类矩阵的条件数随维度增长呈指数级上升。18x18规模下，矩阵已严重病态，微小的浮点计算误差会在求逆过程中被急剧放大，导致不同库的结果出现显著差异。

2. 数值实现细节差异：

   • Eigen的inverse()默认使用部分选主元LU分解，内部针对病态矩阵做了额外优化；
   • Boost uBLAS的lu_factorize默认选主元策略、浮点精度处理逻辑与Eigen不同，在病态矩阵场景下误差累积更明显。

解决方案 #

1. 先验证矩阵数值稳定性 #

计算原矩阵的条件数（无穷范数条件数），确认是否为病态矩阵：

// Eigen中计算条件数
double cond = matrix_Eigen.lpNorm<Eigen::Infinity>() * matrix_Eigen.inverse().lpNorm<Eigen::Infinity>();
std::cout << "Matrix condition number: " << cond << std::endl;

若条件数远大于1e10，说明矩阵严重病态，求逆本身不可靠，需优先从算法层面优化。

2. 优化矩阵构造或算法逻辑 #

避免直接对病态矩阵求逆：

• 改用最小二乘求解（比如求解A*x = b时，直接解线性方程组而非求逆）；
• 对矩阵做预处理（如缩放行/列，降低条件数）；
• 更换问题建模方式，避免生成Vandermonde类矩阵。

3. 提升浮点精度 #

将double替换为long double，用更高精度的浮点类型减少误差累积：

// Eigen版本修改
Eigen::Matrix<long double, Eigen::Dynamic, Eigen::Dynamic> initMatrix_Eigen() {
    Eigen::Matrix<long double, Eigen::Dynamic, Eigen::Dynamic> matrix(SIZE, SIZE);
    for (size_t i = 0; i < SIZE / 2; i++) {
        for (size_t j = 0; j < SIZE; j++) {
            matrix(2 * i, j) = powl((long double)i + 1, (long double)j + 2);
            matrix(2 * i + 1, j) = (j + 2) * powl((long double)i + 1, (long double)j + 1);
        }
    }
    return matrix;
}
// Boost版本同理修改为ublas::matrix<long double>

4. 不依赖特定库的鲁棒求逆实现 #

基于LAPACK标准接口实现跨库兼容的求逆逻辑，LAPACK是数值计算的工业标准，Eigen和Boost底层也依赖它，直接调用可保证计算一致性：

#include <lapacke.h>

bool invertMatrix_LAPACK(double* mat, int size) {
    int* ipiv = new int[size];
    int info = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, size, size, mat, size, ipiv);
    if (info != 0) {
        delete[] ipiv;
        return false;
    }
    info = LAPACKE_dgetri(LAPACK_ROW_MAJOR, size, mat, size, ipiv);
    delete[] ipiv;
    return info == 0;
}

使用时将Eigen或Boost矩阵转换为连续内存数组即可。

5. 验证求逆结果正确性 #

无论用哪个库，都要验证求逆结果的准确性：计算A*A^{-1}与单位矩阵的差异，判断误差是否在可接受范围内：

// Eigen验证
Eigen::MatrixXd identity_check = matrix_Eigen * inverse_Eigen;
std::cout << "Identity check error norm: " << (identity_check - Eigen::MatrixXd::Identity(SIZE, SIZE)).lpNorm<Eigen::Infinity>() << std::endl;

库可靠性结论 #

Eigen和Boost uBLAS都是成熟的数值库，可靠性无本质区别，但Eigen在病态矩阵处理上的优化更友好。核心问题不在库本身，而在于矩阵的病态性——对严重病态矩阵，任何浮点求逆都会产生较大误差，优先从算法层面优化才是根本解决方案。

内容的提问来源于stack exchange，提问作者ConanLord