Eigen::MatrixXd 高性能实战

1，删除Eigen::MatrixXd的某一行或某一列

A是n*n的矩阵，同时删除i行i列
参考 stackoverflow

核心代码实现片段v1
 

        auto removeRow = [](auto& matrix, const Eigen::Index rowTo) {
            Eigen::Index rows = matrix.rows() - 1;
            Eigen::Index cols = matrix.cols();
            if (rowTo < rows) {
                matrix.block(rowTo, 0, rows - rowTo, cols) = matrix.bottomRows(rows - rowTo);
            }
            matrix.conservativeResize(rows, cols);
        };

        auto removeCol = [](Eigen::MatrixXd& matrix, const Eigen::Index colTo) {
            Eigen::Index rows = matrix.rows();
            Eigen::Index cols = matrix.cols() - 1;
            if (colTo < cols) {
                matrix.block(0, colTo, rows, cols - colTo) = matrix.rightCols(cols - colTo);
            }
            matrix.conservativeResize(rows, cols);
        };
        
        Eigen::MatrixXd A1 = A;
        for (const auto& i : zeros) {
            removeRow(A1, i);
            removeClo(A1, i);
        }

removeRow中使用 auto& matrix 是因为 Eigen::VectorXd也可以调用这个lambda

使用block函数，效率还可以，不够好。

参考 Slicing and Indexing

        std::vector<int> noZeros;
        noZeros.reserve(matSize);
        for (int col = 0; col < matSize; ++col) {
            if (!A.row(col).isZero(0)) {
                noZeros.emplace_back(col);
            }
        }
        Eigen::MatrixXd A1 = A(noZeros, noZeros);

这个是最佳实践。

2，求解Ax = b;

参考官方文档 Benchmark of dense decompositions

auto x = A.inverse() * b;//1
Eigen::VectorXd x = A.inverse() * b;//2
Eigen::VectorXd x = A.partialPivLu().solve(b); //3
Eigen::VectorXd x = A.householderQr().solve(b);//4

inverse()的说明

  * For small fixed sizes up to 4x4, this method uses cofactors.

  * In the general case, this method uses class PartialPivLU.

但效率不如直接调用partialPivLu()

最佳实践是 可逆矩阵 partialPivLu()最高效，不可逆 householderQr()最高效。

如果要访问x，使用auto非常低效
参考 C++11 and the auto keyword

《Effective Modern C++》Item6