Why we can't replace all functions accepting Matrix with MatrixBase. Reproducible example inside.

[Daniel Lee]

More discussion has come up on math issue #62: https://github.com/stan-dev/math/issues/62

We can't blindly use MatrixBase even though it seems natural. One example is if you pass a matrix multiplication (we use this in the math library all the time), coefficient access doesn't work. 

I haven't tried finding other operations that dont work. 

Here's a simple example:

#include <Eigen/Dense>
#include <iostream>

template <typename Derived>
void bar(const Eigen::MatrixBase<Derived>& a) {
  std::cout << "a(0): " << a(0) << std::endl;
}

int main() {
  Eigen::MatrixXd u(1,1);
  u << 1;
  Eigen::MatrixXd v(1,2);
  v << 1, 2;

  bar(u*v);
  return 0;

}

Looks simple enough. Compiling and running results in:

> ./a.out
Assertion failed: (this->rows() == 1 && this->cols() == 1), function coeff, file /Users/daniel/dev/stan/lib/eigen_3.2.2/Eigen/src/Core/ProductBase.h, line 150.
a(0): Abort trap: 6

Fun. You might think, why would anyone access anything as an element? Unfortunately, we do. There isn't a list of operations that fails with every template expression, so it's not clear if there's more unsafe behavior lurking.

The response from the Eigen developers is forwarded below. It's a design decision to not allow that particular operation. In Eigen2, there was a way to get around it, but not with Eigen3.

Daniel

Begin forwarded message:

From: Gael Guennebaud <gael.gu...@gmail.com>
Date: March 13, 2015 at 4:57:23 PM EDT
To: eigen <ei...@lists.tuxfamily.org>
Subject: Re: [eigen] help using foo(const MatrixBase<T>&) with (u * v)
Reply-To: ei...@lists.tuxfamily.org

Hi,

as Christoph said evaluating matrix product coefficient-wise is highly inefficient, therefore in Eigen 3 we decided to forbid coefficient accesses to product expression unless this a inner product. This is what the assertion exactly does.

The following:

bar((u*v).transpose())

works because (u*v) gets evaluated into a MatrixXd temporary by the Transpose expression.

If you really want to evaluate the product coefficient-wise, then call:

bar(u.lazyProduct(v)) 

lazyProduct returns a "true" expression template.

if you are going to access to all coefficients then rather call bar like this:

bar((u*v).eval()).

gael

On Fri, Mar 13, 2015 at 3:57 PM, Daniel Lee <bea...@alum.mit.edu> wrote:

Hi Christoph,

Thank you for the response. I think the access within the function masked the real problem we are having.

The issue: matrix multiplication does not behave like other expression templates, e.g. transpose, Map.

I am following instructions here: http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html

Example function (bar.hpp):
------------------------------------------------------------
#include <Eigen/Dense>

#include <iostream>

template <typename Derived>
void bar(const Eigen::MatrixBase<Derived>& a) {
  std::cout << "a(0): " << a(0) << std::endl;
}
------------------------------------------------------------

Example that works:
------------------------------------------------------------

#include "bar.hpp"

int main() {
  Eigen::MatrixXd u(1,1);
  u << 1;
  Eigen::MatrixXd v(1,2);
  v << 1, 2;

  bar((u*v).transpose());
  bar((u*v).transpose().transpose());

  return 0;

}

------------------------------------------------------------

Example that does not work, but I expect to work:
------------------------------------------------------------

#include "bar.hpp"

int main() {
  Eigen::MatrixXd u(1,1);
  u << 1;
  Eigen::MatrixXd v(1,2);
  v << 1, 2;

  bar(u*v);

  return 0;

}

------------------------------------------------------------

This results in:
> ./a.out
Assertion failed: (this->rows() == 1 && this->cols() == 1), function coeff, file /Users/daniel/dev/stan/lib/eigen_3.2.2/Eigen/src/Core/ProductBase.h, line 150.
a(0): Abort trap: 6

I am not expecting this assertion failure. With EIGEN2_SUPPORT defined, this assertion goes away.

Question: did we do something incorrectly with the definition of the function bar?

Thanks again.

Daniel

On Fri, Mar 13, 2015 at 9:33 AM, Christoph Hertzberg <ch...@informatik.uni-bremen.de> wrote:

Am 12.03.2015 um 21:44 schrieb Daniel Lee:

It seems to work fine for transpose, for Map, for MatrixXd, but
here's a minimal example of the problem we're running into with
products:

Accessing products coefficient-wise is generally not a good idea. If you really want to do that, consider calling foo(u.lazyProduct(v),0,0);
If you want to access multiple coefficients from a inside foo, passing via (const Ref<const MatrixXd>&) is probably a better idea.

If I change the #ifdef to an #ifndef, everything works as expected.
Are we somehow supposed to be setting the EIGEN2_SUPPORT property
ourselves?  If I do that in the compilation, it works as expected, but
I get a deprecation warning:

Yes, you can do that, and you can disable the deprecation warning by also defining EIGEN_NO_EIGEN2_DEPRECATED_WARNING.
But as the warning says, this will not be possible with Eigen 3.3 anymore.
See this and the pages linked from there:
http://eigen.tuxfamily.org/dox/Eigen2ToEigen3.html


Christoph


--
----------------------------------------------
Dipl.-Inf., Dipl.-Math. Christoph Hertzberg
Cartesium 0.049
Universität Bremen
Enrique-Schmidt-Straße 5
28359 Bremen

Tel: +49 (421) 218-64252
----------------------------------------------

[Ben Goodrich]

On Thursday, August 18, 2016 at 7:31:46 AM UTC-4, Daniel Lee wrote:

Here's a simple example:

#include <Eigen/Dense>
#include <iostream>

template <typename Derived>
void bar(const Eigen::MatrixBase<Derived>& a) {
  std::cout << "a(0): " << a(0) << std::endl;
}

This is not the recommended way to access elements of an Eigen object anyway. If you change it to

std::cout << "a(0): " << a.coeff(0) << std::endl;

then it does not segfault.

Ben

[Daniel Lee]

On Thu, Aug 18, 2016 at 10:29 AM, Ben Goodrich <goodri...@gmail.com> wrote:

On Thursday, August 18, 2016 at 7:31:46 AM UTC-4, Daniel Lee wrote:

Here's a simple example:

#include <Eigen/Dense>
#include <iostream>

template <typename Derived>
void bar(const Eigen::MatrixBase<Derived>& a) {
  std::cout << "a(0): " << a(0) << std::endl;
}

This is not the recommended way to access elements of an Eigen object anyway.

Really? Where do you find that sort of doc? I can't find it on Eigen's page.

 

If you change it to

std::cout << "a(0): " << a.coeff(0) << std::endl;

then it does not segfault.

Ben

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Ben Goodrich]

On Thursday, August 18, 2016 at 10:38:38 AM UTC-4, Daniel Lee wrote:

On Thu, Aug 18, 2016 at 10:29 AM, Ben Goodrich <goodri...@gmail.com> wrote:

On Thursday, August 18, 2016 at 7:31:46 AM UTC-4, Daniel Lee wrote:

Here's a simple example:

#include <Eigen/Dense>
#include <iostream>

template <typename Derived>
void bar(const Eigen::MatrixBase<Derived>& a) {
  std::cout << "a(0): " << a(0) << std::endl;
}

This is not the recommended way to access elements of an Eigen object anyway.

Really? Where do you find that sort of doc? I can't find it on Eigen's page.

I don't know. It is used in the Eigen code a lot. They mention it here

http://eigen.tuxfamily.org/dox/group__QuickRefPage.html

The way you originally had it, it checks the range. But since we are doing range checking in the Stan Math functions anyway, it is not worth triggering the segfault.

Ben

[Daniel Lee]

I didn't see anywhere where they recommended coeff() for element access. Did I miss it somewhere?

[Ben Goodrich]

On Thursday, August 18, 2016 at 11:45:08 AM UTC-4, Daniel Lee wrote:

I didn't see anywhere where they recommended coeff() for element access. Did I miss it somewhere?

Well, they actually recommend () to get the range check in user code. But they do .coeff( 1403 times in Eigen itself because they don't want to do range checks. If we range check beforehand, we can get away with using .coeff() for read and .coeffRef() for write. There may be other issues that we are not anticipating, but it doesn't seem as if this one thing should force us to wontfix the issue.

Ben

[Bob Carpenter]

Right, and same issue with all the use of std::vector::at(size_t) I've
been seeing in our code. We should be checking sizes ourselves, not
relying on std::vector to throw exceptions.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.