New user self-registration is disabled due to spam. Please email eigen-core-team @ if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
Bug 1543 - "Fix linear indexing in generic block evaluation" breaks Matrix*Diagonal*Vector product
Summary: "Fix linear indexing in generic block evaluation" breaks Matrix*Diagonal*Vect...
Alias: None
Product: Eigen
Classification: Unclassified
Component: Core - expression templates (show other bugs)
Version: 3.4 (development)
Hardware: x86 - general All
: Highest Wrong Result
Assignee: Nobody
Depends on:
Blocks: 3.3 3.4
  Show dependency treegraph
Reported: 2018-04-15 00:47 UTC by Linus Mårtensson
Modified: 2018-07-10 07:17 UTC (History)
3 users (show)

Minimal reproduction showing error using real life data (19.66 KB, text/x-c++src)
2018-04-15 00:47 UTC, Linus Mårtensson
no flags Details

Description Linus Mårtensson 2018-04-15 00:47:22 UTC
Created attachment 850 [details]
Minimal reproduction showing error using real life data


LDLT convergence loop fails to converge after commit

Steps to reproduce

Build and run l.cpp with:
g++ l.cpp -Ipath/to/Eigen && ./a.out

Actual results:

minimizeGroupL1 failed to converge in 100 iterations

Expected results:

minimizer successful!

Build date & Hardware:

Build 2018-04-15 on Ubuntu 17.10

Additional information:

Still present on master as of today.
Reverting the specified commit resolves the issue.
Comment 1 Christoph Hertzberg 2018-04-15 15:05:45 UTC
The problem seems to be already in the Matrix*Diagonal*Vector product. The following minimal example also exposes the problem:

int main()
    Eigen::MatrixXd A; A.setRandom(100,4);
    Eigen::VectorXd d,e; d.setRandom(100), e.setRandom(100);
    Eigen::MatrixXd D = d.asDiagonal();
    Eigen::VectorXd B1 = A.transpose()*d.asDiagonal()*e;
    Eigen::VectorXd B2 = A.transpose()*D*e;
    std::cout << B1 << "\n\n" << B2 << "\n";

The currently failing sparse_product and permutationmatrices unit tests may be related to this.
Gael: Can you have a look at this?
Comment 2 Christoph Hertzberg 2018-04-19 12:05:46 UTC
The issue seems to be that `A.transpose()*d.asDiagonal()*e;` is evaluated (slightly simplified) like so:

for(Index i=0; i<rows; ++i)
   dest.coeffRef(i) = (lhs.row(i).cwiseProduct(actual_rhs.transpose())).sum();

Where `lhs` is `A.transpose()*d.asDiagonal()`, and the corresponding `.row(i)` claims to be a block which has linear access. However, for the (dense*diagonal).row(i) this evaluates to:
  A.packet(linear_start + j) * d.packet(linear_start + j),
and the latter is obviously wrong.

A possible solution would be to make sure that `Block<Xpr,... >` only sets the `InnerPanel` flag, if `Xpr` has the LinearAccessBit set. Unfortunately, setting this seems not to be done consistently -- sometimes it is done in `traits<Xpr>::Flags`, sometimes in `traits<Xpr>::EvaluatorFlags` and sometimes in `evaluator<Xpr>::Flags`.
Comment 3 Gael Guennebaud 2018-04-23 12:23:18 UTC
An even simpler test is:

  sq_m2 = v1.asDiagonal();
  sq_m2 = sq_m1 * sq_m2;
  VERIFY_IS_APPROX( (sq_m1*v1.asDiagonal()).col(i), sq_m2.col(i) );

I added it in test/diagonalmatrices.cpp
Comment 4 Gael Guennebaud 2018-04-23 12:45:24 UTC
Fix in default and 3.3 branches:

It is simply that the condition for forwarding linear indexing was not strict enough. Hopefully, I did not overseen anything this time.
Comment 5 Gael Guennebaud 2018-04-23 12:51:35 UTC
hm, perhaps it would be better to simply define

   ForwardLinearAccess = bool(evaluator<ArgType>::Flags&LinearAccessBit)

thus ignoring the InnerPanel flag which seems irrelevant here. This way (vec1+vec2).segment(i,n) would also forward linear indices.
Comment 6 Christoph Hertzberg 2018-04-23 14:22:19 UTC
I assume the idea of the InnerPanel flag was that, e.g., A.middleCols(...) can be accessed linearly (if A itself is column-major and can be accessed linearly), but A.block(...) cannot.
Comment 7 Gael Guennebaud 2018-04-23 15:39:55 UTC
For the record, there are 3 cases:

1) Block<Sparse>: enable write access 
2) Block<Dense,DirectAccess>: linearization+alignment
3) Block<Dense,Generic>: linearization

For cases 2 and 3, this is a welcome information when dealing with small matrices. For 1), this information is vital.
Comment 8 Gael Guennebaud 2018-07-10 07:17:19 UTC
I improved the ForwardLinearAccess flag to activate it for any vector-expression:

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