了解 Eigen 中的 solveInPlace 操作

Question

我在 Eigen3.3.7 中使用 LLT 时尝试探索 "solveInPlace()" 函数的选项以加速我应用程序中的矩阵逆计算。我用下面的代码来测试它。

    int main()
    {
        const int M=3;

        Eigen::Matrix<MyType,Eigen::Dynamic,Eigen::Dynamic> R = Eigen::Matrix<MyType,Eigen::Dynamic,Eigen::Dynamic>::Zero(M,M);
        // to make sure full rank
        for(int i=0; i<M*2; i++)
        {
            const Eigen::Matrix<MyType, Eigen::Dynamic,1> tmp = Eigen::Matrix<MyType,Eigen::Dynamic,1>::Random(M);
            R += tmp*tmp.transpose();
        }

        std::cout<<"R \n";
        std::cout<<R<<std::endl;
        decltype (R) R0 =  R; // saving for later comparison



        Eigen::LLT<Eigen::Ref<Eigen::Matrix<MyType,Eigen::Dynamic,Eigen::Dynamic> > > myllt(R);
        const Eigen::Matrix<MyType,Eigen::Dynamic,Eigen::Dynamic> I = Eigen::Matrix<MyType,Eigen::Dynamic,Eigen::Dynamic>::Identity(R.rows(), R.cols());

        myllt.solveInPlace(I);

        std::cout<<"I: "<<I<<std::endl;
        std::cout<<"Prod InPlace: \n"<<R0*I<<std::endl;


        return 0;
}

阅读 Eigen 文档后，我认为输入矩阵（这里 "R"）会在计算变换时被修改。令我惊讶的是，我发现结果存储在 "I" 中。这不是预期的，因为我将 "I" 定义为常量。请为此行为提供解释。

Answer 1

简单的非编译器答案是您要求 LLT 解决就地（即在传递的参数中）那么您期望结果如何成为？显然，您会认为这是一个编译器错误，因为 "in-place" 意味着更改参数，但您传递的是一个常量对象。

因此，如果我们在 Eigen 文档中搜索 solveInPlace，我们会发现唯一一个采用 const 引用的项目具有 following note:

"in-place" version of TriangularView::solve() where the result is written in other

Warning
The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here. This function will const_cast it, so constness isn't honored here.

非就地选项为：

R = myllt.solve(I);

但这并不能真正加快计算速度。在任何情况下，在您决定需要就地选项之前进行基准测试。

你的问题已经到位，因为 const_cast 的目的是剥离 references/pointers 它们的 const-ness iff 基础变量不是 const 限定*（cppref).如果你要写一些例子

const int i = 4;
int& iRef = const_cast<int&>(i); // UB, i is actually const
std::cout << i; // Prints "I want coffee", or it can as we like UB
int j = 4;
const int& jRef = j;
const_cast<int&>(jRef)++; // Legal. Underlying variable is not const.
std::cout << j; // Prints 5

i 的情况可能会按预期工作或不正常，我们依赖于每个 implementation/compiler。它可能适用于 gcc，但不适用于 clang 或 MSVC。没有任何保证。当您在示例中间接调用 UB 时，编译器可以选择执行您期望的操作或完全执行其他操作。

*从技术上讲，它是 UB 的修改，而不是 const_cast 本身。

了解 Eigen 中的 solveInPlace 操作

Understanding solveInPlace operation in Eigen

matrix-inverse

eigen

eigen3