了解 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 就地解决（即在传递的参数中），那么您期望结果是什么？显然，您会认为这是编译器错误，因为“就地”意味着更改参数，但您传递的是 const 对象。

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

“就地”版本的 TriangularView::solve() 结果写入 other

警告
该参数仅标记为“const”以使 C++ 编译器在此处接受临时表达式。此函数将对其进行 const_cast，因此此处不支持 const 性。

非就地选项是：

R = myllt.solve(I);

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

您的问题已经到位，因为 const_cast 的意思是如果基础变量不是 const 限定的 * (cppref)，则剥离引用/指针的 const-ness。如果你要写一些例子

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 的情况可能会像预期的那样工作，我们依赖于每个实现/编译器。它可能适用于 gcc，但不适用于 clang 或 MSVC。没有任何保证。由于您在示例中间接调用了 UB，因此编译器可以选择执行您期望的操作或完全执行其他操作。

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