Why does Eigen make inconsistent default assumptions about aliasing?

Question

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?

Humming

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 – Manrope

Essentially it is a design choice. As @Manrope 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. – Fixed

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) – Manrope

@Manrope 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.) – Humming

@Manrope 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. – Humming

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 @Fixed said. Therefore, a = a.transpose() could be discouraged anyway – Manrope

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 – Manrope

@Manrope 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. – Humming

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. – Manrope

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). – Campuzano

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. – Campuzano

@Fixed I was going to say I need more convincing but now we've heard it from the horse's mouth. – Humming

@Manrope 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. – Fixed

Answer 1

The comments present an answer to the question and are summarized here:

Eigen documentation 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. 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.

Essentially it is a design choice. 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.

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).

This decision was made at the very early stage of Eigen developments, 10 years later I'd rather assume no aliasing unconditionally for consistency.

Nodical