This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
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...
Status: RESOLVED FIXED
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
URL:
Whiteboard:
Keywords:
Depends on:
Blocks: 3.3 3.4
  Show dependency treegraph
 
Reported: 2018-04-15 00:47 UTC by Linus Mårtensson
Modified: 2019-12-04 17:38 UTC (History)
3 users (show)



Attachments
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

Overview

LDLT convergence loop fails to converge after commit
https://bitbucket.org/eigen/eigen/commits/4aa7f03ab605a48e9af22fcbf5040ffd51d1f442

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_m1.setRandom();
  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:

https://bitbucket.org/eigen/eigen/commits/728dfff683d6/
https://bitbucket.org/eigen/eigen/commits/5967bc3c2cdb/

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:
https://bitbucket.org/eigen/eigen/commits/a8600da73e11/
Comment 9 Nobody 2019-12-04 17:38:10 UTC
-- GitLab Migration Automatic Message --

This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity.

You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/1543.

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