CODEWITHSUNDEEP
c++matrixeigen3
PG_LoLo
PG_LoLo

Reputation: 178

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.

Upvotes: 4

Views: 396

Answers (2)

pingul
pingul

Reputation: 3473

"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();

Upvotes: 3

chtz
chtz

Reputation: 18807

pingul's answer describes why your solution fails, I just wanted to add what you can do, if you want to avoid the temporary. Essentially, the trick is to assign the expression just to one triangular half of your matrix:

#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';
}

In fact, if you just want 'a random symmetric matrix' and not necessarily the symmetric part of a given matrix M, you can simply copy the upper part to the lower part (or the other way around):

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

Upvotes: 2

Related Questions