11

As a newcomer to Eigen there is something I am battling to come to terms with.

With matrix multiplication Eigen creates a temporary by default to avoid aliasing issues:

matA = matA * matA; // works fine (Eigen creates a temporary before assigning)

If it is safe to assume no aliasing, we can use .noalias() to allow for optimization:

matB = matA * matA; // sub-optimal (due to unnecessary temporary)
matB.noalias() = matA * matA; // more efficient

So Eigen here by default avoids making the unsafe assumption unless it is explicitly told that it is safe to assume no aliasing. So far so good.

However, for many other expressions, such as a = a.transpose() Eigen makes the unsafe assumption (that there are no aliasing issues) by default, and needs explicit intervention to avoid that unsafe assumption:

a = a.transpose().eval(); // or alternatively: a.transposeInPlace()

So if we are concerned about efficiency, it is not sufficient to be careful when there is potential aliasing, and it is not sufficient to be careful when there is no potential aliasing. We have to be careful both when there is potential aliasing, and when there is no potential aliasing, and decide whether a special intervention is warranted, based on whether the expression involves matrix multiplication or not. Is there some design rationale in Eigen for this "mixture of default expectations" in Eigen's interface?

Museful
  • 6,711
  • 5
  • 42
  • 68
  • 1
    I just checked on Eigen page. They mentioned effectively that *you should not do that* `a = a.transpose()`, but they also mentioned that *Eigen uses a run-time assertion to detect this and exits with a message*. I guess that the problem is that as they state, *In general, aliasing cannot be detected at compile time*, and they favor efficiency. Performing `b = a.transpose()` while creating a temporary would be quite inefficient – Damien May 16 '19 at 14:26
  • 4
    Essentially it is a design choice. As @Damien said, in general Eigen prefers efficiency, but matrix-products are a common source of accidental aliasing (and the extra copy usually is far less expensive than doing the product). For transpose better use `a.transposeInPlace();`. Most other expressions barely have aliasing issues, e.g., when adding two matrices you can only get problems if you access overlapping (but non-equal) memory regions. – chtz May 16 '19 at 14:42
  • Another aspect is that performing multiplication costs about O(n^3), while creating a temporary costs O(n^2). For the transpose operation, both transposition and copying costs O(n^2) – Damien May 16 '19 at 14:43
  • @Damien Then the question is why they no longer favor efficiency (over safety by default) when it comes to matrix multiplication. (O(n^3) is brute-force matrix multiplication; surely they re-use some sub-expressions for better complexity in general.) – Museful May 16 '19 at 14:47
  • @Damien Eigen's page says: ""Eigen does detect ***some*** instances of aliasing, albeit at run time" [emphasis added]. It sounds like the onus is still on the user to ensure it doesn't happen. – Museful May 16 '19 at 14:53
  • Even if multiplication is not really O(n^3), it should be much more costly than copying. Another point is that `a.transposeInPlace()` exists as @chtz said. Therefore, `a = a.transpose()` could be discouraged anyway – Damien May 16 '19 at 14:53
  • My understanding is that an aliasing like `a = a.transpose()` should be detected easily. But all that are hypothesis. Only *Eigen* developers could provide a real answer – Damien May 16 '19 at 14:57
  • @Damien I know that `a.transposeInPlace()` exists, and that `a = a.transpose()` is incorrect. I know it is possible to write correct code in Eigen. My question is whether there is rationale behind the apparent inconsistency in their interface design. – Museful May 16 '19 at 14:58
  • As you noticed yourself, it seems there is not a strict rationale. I assume a case per case decision, weighting the possible efficiency loss in one hand, and the security risk in the other hand. – Damien May 16 '19 at 15:15
  • 3
    As chtz and Damien said, the rationale is simply that products are much more expensive than temporary creation + O(n^2) copy. And yes, Eigen, just like any other BLAS library, use a O(n^3) matrix product algorithm. Theoretically more efficient algorithms payoff for huge products (n>3000 maybe even more), and more importantly, they yield to huge and unacceptable roundoff errors with floating point arithmetic (float/double). – ggael May 16 '19 at 15:16
  • 4
    btw, this decision was made at the very early stage of Eigen developments, 10 years later I'd rather assume no aliasing unconditionally for consistency. – ggael May 16 '19 at 15:19
  • @chtz I was going to say I need more convincing but now we've heard it from the horse's mouth. – Museful May 16 '19 at 15:27
  • @Damien Regarding detection of `a = a.transpose()`: This is in fact done at runtime (unless assertions are disabled), but detecting this at compile-time is generally not possible (just consider the case where you pass two matrices by reference, which could but don't have to be the same). Also, what is hard to detect (even at runtime) is aliasing between sub-blocks of matrices. – chtz May 16 '19 at 16:03

0 Answers0