New user self-registration is currently disabled. Please email eigen-core-team @ if you need an account.
Bug 1435 - Matrix subtraction gives wrong result
Matrix subtraction gives wrong result
Product: Eigen
Classification: Unclassified
Component: Core - general
3.3 (current stable)
x86 - 64-bit Linux
: Normal Wrong Result
Assigned To: Nobody
Depends on:
  Show dependency treegraph
Reported: 2017-06-07 13:29 UTC by Andrew Palmer
Modified: 2017-06-08 10:57 UTC (History)
3 users (show)

Example code (3.31 KB, application/gzip)
2017-06-07 13:29 UTC, Andrew Palmer
no flags Details

Description Andrew Palmer 2017-06-07 13:29:10 UTC
Created attachment 788 [details]
Example code

This gives the wrong results:

P = P - K * H * P;

While this gives the correct result:

P -= K * H * P;

Please see the attached code for an example. The correct result should be:

0.677419 0.322581        0
0.322581 0.877419        0
       0        0      1.3

but the incorrect result is:

0.677419 0.322581        0
0.322581  1.09594        0
       0        0      1.3

Please let me know if you need any additional information.
Comment 1 Gael Guennebaud 2017-06-07 14:51:27 UTC
This is an aliasing issue:

Basically, P = P - K * H * P; is decomposed as:

P = P;
T1.noalias() = K * H;
P.noalias() -= T1 * P; // aliasing issue

instead of:

T1.noalias() = K * H;
T2.noalias() = T1 * P;
P = P - T2; // no aliasing

thus saving one temporary. Unfortunately, this rewriting (introduced in some 3.3.x release) can introduce additional aliasing issues. I'm not sure what to do about it. Would it be good enough to detect them an trigger a runtime assertion in debug mode? similar to what we do with A = A.transpose() ??
Comment 2 Christoph Hertzberg 2017-06-07 15:41:53 UTC
So a workaround in this case would actually be to write:
  P = P - K*(H*P);

But writing something like this would actually be even worse?
  P = A - B * C * P;

Generally, I think we should make sure that as long as products are involved (and the user did not write .noalias()) we should assume that aliasing may happen.
Making noalias() instead of may_alias the default behavior would be a huge API-break, IMO.

Except for the memory allocation (of T0), do we actually save anything in comparison to

  T0 = P; T1 = K*H; T0-=T1*P;
Comment 3 Gael Guennebaud 2017-06-08 09:22:43 UTC
Ok, so this rewriting rule should be enabled only if the destination is flagged as no-alias, otherwise it is indeed useless. Let's consider a more canonical cases:

(1) R = C + A*B;
(2) R.noalias() = C + A*B;

without rewriting, there is no aliasing issue, noalias() is ignored and in both cases we have only one temporary:

T0 = A*B;
R = C + T0;

with a careful rewriting we would generate:

for case (2):

R = C;
R += A*b;

-> zero temporary

for case (1):

T0 = C;
T0 += A*B;

-> one temporary and the swap might be costly if R is a block expression.
Comment 4 Gael Guennebaud 2017-06-08 10:57:51 UTC
nevermind, I've already handled this issue with operator +, and it's just that I have forgotten to a specialization when I added the rewriting for operator -, fixed: (devel) (3.3)

Note You need to log in before you can comment on or make changes to this bug.