c++ Eigen3 矩阵奇怪的行为

Question

我正在尝试使用线性代数 C++ 库 Eigen3 获取随机对称矩阵。我是这样做的：

Eigen::MatrixXd m(3, 3);
m.setRandom();
m = 0.5 * (m + m.transpose());

但是结果完全错了。但是，如果我不重写 m 变量，而是像这样简单地将它输出到控制台：

Eigen::MatrixXd m(3, 3);
m.setRandom();
cout << 0.5 * (m + m.transpose()) << endl;

似乎一切正常。我不明白问题出在哪里。是不是因为像转置这样的方法和像 * 和 + 这样的操作不会立即创建一个新矩阵，而是以一种惰性的方式创建新矩阵并持有对矩阵 m 的引用？但是我应该如何从官方文档中知道呢？像这样的行为不是极易出错吗？

更新：是的，我认为我对惰性计算的猜测是正确的。 transpose 方法的文档中提到了它：

/** \returns an expression of the transpose of *this.
  *
  * Example: \include MatrixBase_transpose.cpp
  * Output: \verbinclude MatrixBase_transpose.out
  *
  * \warning If you want to replace a matrix by its own transpose, do \b NOT do this:
  * \code
  * m = m.transpose(); // bug!!! caused by aliasing effect
  * \endcode
  * Instead, use the transposeInPlace() method:
  * \code
  * m.transposeInPlace();
  * \endcode
  * which gives Eigen good opportunities for optimization, or alternatively you can also do:
  * \code
  * m = m.transpose().eval();
  * \endcode
  *
  * \sa transposeInPlace(), adjoint() */

所以现在我想知道在执行长链计算时应该使用什么模式？到处写.eval()？老实说，它非常丑陋，而且仍然容易出错。

Answer 1

“我认为我对惰性计算的猜测是正确的。”

是的，你是对的。 here 描述了惰性求值的规则。我提取了以下几点：

Eigen determines automatically, for each sub-expression, whether to evaluate it into a temporary variable. [...]

Expression-templates-based libraries can avoid evaluating sub-expressions into temporaries, which in many cases results in large speed improvements. This is called lazy evaluation as an expression is getting evaluated as late as possible, instead of immediately. However, most other expression-templates-based libraries always choose lazy evaluation. There are two problems with that: first, lazy evaluation is not always a good choice for performance; second, lazy evaluation can be very dangerous, for example with matrix products: doing matrix = matrix*matrix gives a wrong result if the matrix product is lazy-evaluated, because of the way matrix product works.

"所以现在我想知道在执行长链计算时应该使用什么模式？"

一般来说，表达式模板应该解决是否评估 immediately/lazily 的问题，但正如您所注意到的，有时它们并不能找到所有边缘情况，而 matrix.transpose() 似乎是一个。您可以添加 .eval() 和 .noalias() 以强制延迟或立即评估，否则会选择另一个：

强制惰性求值：matrix1.noalias() = matrix2 * matrix2;惰性求值是可以的——matrix1 和 matrix2 之间没有别名
强制立即求值：matrix1 = matrix1.transpose().eval()惰性求值不行

对于您的转置案例，只需添加 .eval():

m = 0.5 * (m + m.transpose()).eval();

Answer 2

pingul 的回答描述了您的解决方案失败的原因，我只是想添加您可以做的事情，如果您想避免暂时的。本质上，诀窍是将表达式仅分配给矩阵的一个三角形一半：

#include <iostream>
#include <Eigen/Core>

int main() {
    Eigen::Matrix4f M;
    M.setRandom();
    std::cout << "M:\n" << M << '\n';
    // for comparison, evaluate into other matrix (similar to using .eval()):
    Eigen::Matrix4f M2 = 0.5f*(M+M.transpose());
    // in place operation, only assing to lower half, then to upper half:
    M.triangularView<Eigen::StrictlyLower>() = 0.5f*(M + M.transpose());
    M.triangularView<Eigen::StrictlyUpper>() = M.transpose();
    std::cout << "M2:\n" << M2 << "\nM:\n" << M << '\n';
}

事实上，如果你只是想要'a random symmetric matrix'而不一定是给定矩阵的对称部分M，你可以简单地将上半部分复制到下半部分（或者反过来):

    M.triangularView<Eigen::StrictlyLower>() = M.transpose();

c++ Eigen3 矩阵奇怪的行为

c++ Eigen3 matrix strange behaviour

c++

matrix

eigen3