C++ 近奇异矩阵求逆

在C++中使用Eigen库，想实现C=A/B的一个矩阵除法，其中B是近奇异矩阵，值很小，用Eigen中的inverse()函数求的逆是inf，请问该怎么去得到这个矩阵的逆的实现呢？？

使用矩阵分解，如QR分解或SVD分解

#include <Eigen/Dense>
Eigen::MatrixXd matrixDivision(const Eigen::MatrixXd& A, const Eigen::MatrixXd& B)
{
    Eigen::MatrixXd result;
    
    // SVD decomposition
    Eigen::JacobiSVD<Eigen::MatrixXd> svd(B, Eigen::ComputeThinU | Eigen::ComputeThinV);
    const Eigen::MatrixXd& U = svd.matrixU();
    const Eigen::MatrixXd& V = svd.matrixV();
    const Eigen::VectorXd& singularValues = svd.singularValues();
    
    // Pseudo-inverse calculation
    Eigen::MatrixXd B_pinv = V * singularValues.asDiagonal().inverse() * U.transpose();
    
    // Matrix division
    result = A * B_pinv;
    
    return result;
}

上面用SVD分解对矩阵B进行分解，并根据奇异值计算其伪逆矩阵B_pinv。然后，将矩阵A与B的伪逆相乘，得到矩阵除法的结果。
请注意，伪逆矩阵的计算可能会引入数值误差，因此结果可能不会完全精确。

如果矩阵B是近奇异矩阵，那么它的逆矩阵会非常大，这就导致了用Eigen库的inverse()函数计算出来的结果是inf。

一种常见的解决方法是使用奇异值分解（SVD）来求B的伪逆矩阵，然后再用C++中的乘法运算实现矩阵除法。

在Eigen库中，可以使用JacobiSVD类进行奇异值分解。下面是一个简单的示例代码：

#include <Eigen/Dense>
using namespace Eigen;

// ...

MatrixXd A, B, C;
// 初始化A、B、C矩阵

JacobiSVD svd(B, ComputeThinU | ComputeThinV);

double tolerance = 1e-8; // 可以根据实际情况调整容差
MatrixXd B_pinv = svd.matrixV() * DiagonalMatrix(svd.singularValues().array().abs().max(tolerance).inverse()) * svd.matrixU().transpose();

C = A * B_pinv;
其中，svd是使用JacobiSVD类对矩阵B进行奇异值分解得到的对象， B_pinv是B的伪逆矩阵。

这里的tolerance表示在奇异值分解过程中，当某个奇异值小于等于该值时，认为该奇异值为0。在上面的代码中，使用了max(tolerance)函数将小于该容差的奇异值设为0，然后再取其倒数。这样可以避免除以0的错误。

B是近奇异矩阵，导致其逆矩阵无法准确求出。这种情况下可以使用比较稳定的伪逆矩阵（pseudo-inverse）来代替准确的逆矩阵。
在Eigen库中，可以使用 BDCSVD 类来求取伪逆矩阵。具体实现方法如下：

// 假设 A 和 B 分别为 Eigen::Matrix<double, n, m> 类型的两个矩阵
Eigen::Matrix<double, m, m> BTB = B.transpose() * B;
Eigen::Matrix<double, m, m> invBTB = BTB.inverse();

// 计算伪逆矩阵
// 这里使用了 BDCSVD 类来计算伪逆矩阵
Eigen::BDCSVD<Eigen::Matrix<double, m, m>> svd(BTB);
double tol = 1e-8;  // 阈值可以根据实际情况进行调整
Eigen::Matrix<double, m, m> pinvBTB = svd.matrixV() * 
    invBTB.diagonal().array().abs().sqrt().matrix().asDiagonal() *
    svd.matrixU().transpose();

// 计算 C 矩阵
Eigen::Matrix<double, n, m> C = A * pinvBTB * B.transpose();

先计算了 B 转置乘以 B 的结果 BTB，然后对 BTB 求正常的逆矩阵 invBTB。接下来使用 BDCSVD 类计算 BTB 的奇异值分解，然后基于奇异值分解结果求取伪逆矩阵 pinvBTB。最后就可以使用 C = A * pinvBTB * B.transpose() 计算 C 矩阵了。