Search code examples
c++matrixeigen3

c++ Eigen3 matrix strange behaviour


I'm trying to get a random symmetric matrix with linear algebra c++ library Eigen3. I'm doing it like this:

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

But the reasult is totally wrong. But if I won't rewrite the m variable and simply output it to console like this:

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

All seems to work properly. I can not understand where is the problem. Is it because methods like transpose and operations like * and + do not instantly create a new matrix as instead doing it in a lazy manner and holding a reference to matrix m? But how then should I know it from official documentation? And isn't behaviour like this overwhelmingly error-prone?

UPDATE: Yes, I think my guess about lazy calculations is correct. It is mentioned in the docummentation of the transpose method:

/** \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() */

So now I am wondering what patterns should I use when I perform long chains of calculations? Everywhere writing .eval()? To be honest, it's pretty ugly and still error-prone.


Solution

  • "I think my guess about lazy calculations is correct."

    Yes, you are right. The rules for the lazy evaluation is described here. I've extracted some of the points below:

    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.

    "So now I am wondering what patterns should I use when I perform long chains of calculations?"

    In general, the expression templates should solve the issue of whether to evaluate immediately/lazily, but as you note sometimes they don't find all edge cases, and matrix.transpose() seems to be one. You can add .eval() and .noalias() to force lazy or immediate evaluation, where the other would otherwise have been chosen:

    • Force lazy evaluation: matrix1.noalias() = matrix2 * matrix2; Lazy evaluation is OK -- no aliasing between matrix1 and matrix2

    • Force immediate evaluation: matrix1 = matrix1.transpose().eval() Lazy evaluation is not OK

    For your transpose case, just add .eval():

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