New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
Bug 259 - Aliasing in combound assignment + reduction
Summary: Aliasing in combound assignment + reduction
Status: ASSIGNED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Core - general (show other bugs)
Version: 3.0
Hardware: All All
: --- Unknown
Assignee: Hauke Heibel
URL:
Whiteboard:
Keywords:
Depends on:
Blocks:
 
Reported: 2011-05-03 12:07 UTC by Hauke Heibel
Modified: 2011-09-12 11:04 UTC (History)
4 users (show)



Attachments
Nest expressions by value and everything else by reference. (2.69 KB, patch)
2011-05-12 09:55 UTC, Hauke Heibel
no flags Details | Diff
Nest expressions as PlainObjects and everything else by reference. (2.15 KB, patch)
2011-05-12 13:22 UTC, Hauke Heibel
no flags Details | Diff
Fixed a typo. (2.15 KB, patch)
2011-05-12 20:55 UTC, Hauke Heibel
hauke.heibel: review+
Details | Diff

Description Hauke Heibel 2011-05-03 12:07:46 UTC
Observed: The code below does not demean the data.
Expected: Printout of approx. [0 0]


#include "Eigen/Core"

#include <iostream>

using namespace Eigen;

typedef Matrix<float,Dynamic,Dynamic> MatrixType;

int main(int argc, char* argv[])
{
  MatrixXd m = MatrixXd::Random(5,2);
  m -= m.colwise().mean().replicate(m.rows(),1);
  std::cout << m.colwise().mean() << std::endl;

  m = MatrixXd::Random(5,2);
  m.colwise() -= m.colwise().mean();
  std::cout << m.colwise().mean() << std::endl;
}

Can we detect this and automatically compute a minimal temporary for the RHS. In these cases, I think one really wants a temporary.
Comment 1 Hauke Heibel 2011-05-03 12:12:17 UTC
It should be:

m.rowwise() -= m.colwise().mean();
Comment 2 Gael Guennebaud 2011-05-08 22:30:14 UTC
The aliasing pb has to be expected here and I'm not sure we want to deal with this particular case of aliasing. On the other hand your example shows something more interesting: Replicate should really evaluate the colwise means into a temporary, not for aliasing, but for perfomrance reason to avoid multiple evaluation of it.

btw, the second version is not correct, it should be a .rowwise() on the left hand side:

m.rowwise() -= m.colwise().mean().eval();
Comment 3 Hauke Heibel 2011-05-11 09:28:42 UTC
(In reply to comment #2)
> The aliasing pb has to be expected here and I'm not sure we want to deal with
> this particular case of aliasing. On the other hand your example shows
> something more interesting: Replicate should really evaluate the colwise means
> into a temporary, not for aliasing, but for perfomrance reason to avoid
> multiple evaluation of it.

> m.rowwise() -= m.colwise().mean().eval();

That was a typo, which I fixed in the second comment.

I also thought about the speed implications. Fixing them would automatically resolve the aliasing issue.

If I am not wrong, a column-wise -= is implemented like this (VectorwiseOp.h, ll.438ff):

    template<typename OtherDerived>
    ExpressionType& operator-=(const DenseBase<OtherDerived>& other)
    {
      EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
      for(Index j=0; j<subVectors(); ++j)
        subVector(j) -= other.derived();
      return const_cast<ExpressionType&>(m_matrix);
    }

Here, I think that it might also be a good idea to evaluate "other.derived()" into a temporary. Do you agree?

Evaluating column-wise and replicate arguments into temporaries will resolve both, aliasing as well as speed issues.

- Hauke
Comment 4 Gael Guennebaud 2011-05-11 22:46:58 UTC
yes, that's what I meant. This is just a matter of using nested<>.
Comment 5 Hauke Heibel 2011-05-12 09:54:33 UTC
(In reply to comment #4)
> yes, that's what I meant. This is just a matter of using nested<>.

I have prepared a patch. It contains some code duplication but I wanted to check back before cleaning up.

I think it is not just a matter of using nested<>. In fact, our replicate expression is using nested. I think we need exactly the opposite of nested<>. If I am not completely wrong, we need to nest every expression by value and every plain object by reference.

So right now, I added something like this

typedef typename internal::eval<MatrixType>::type PlainObjectType;
typedef typename internal::conditional<
  internal::is_same<PlainObjectType,MatrixType>::value,
  PlainObjectType&,
  PlainObjectType
>::type NestedType;

const NestedType m_matrix;

to the replicate expressions and something similar to += and -= of VectorwiseOp.

My patch is attached and some feedback would be great.

- Hauke
Comment 6 Hauke Heibel 2011-05-12 09:55:30 UTC
Created attachment 167 [details]
Nest expressions by value and everything else by reference.
Comment 7 Jitse Niesen 2011-05-12 12:42:39 UTC
(In reply to comment #5)
> I think it is not just a matter of using nested<>. In fact, our replicate
> expression is using nested. I think we need exactly the opposite of nested<>.
> If I am not completely wrong, we need to nest every expression by value and
> every plain object by reference.

I'm confused how this is the opposite of nested<>. I thought that nested (with a high enough value for n) does precisely that: expressions are evaluated and nested by value and plain objects are nested by reference.

As you say, Replicate already uses nested<> but with the default value of n = 1. Perhaps simply changing this will solve the issue? That is, change Replicate.h:51 to 

  /* untested code */
  enum { ColTimesRowFactor = (ColFactor == Dynamic || RowFactor == Dynamic)
                             ? Dynamic : ColFactor * RowFactor };
  typedef typename nested<MatrixType,ColTimesRowFactor>::type MatrixTypeNested;


> So right now, I added something like this
> 
> typedef typename internal::eval<MatrixType>::type PlainObjectType;
> typedef typename internal::conditional<
>   internal::is_same<PlainObjectType,MatrixType>::value,
>   PlainObjectType&,
>   PlainObjectType
> >::type NestedType;

If T is Matrix or Array, then eval<T>::type already is a reference to the Matrix or Array. So why do you need the conditional?
Comment 8 Hauke Heibel 2011-05-12 13:16:32 UTC
Thanks for the feedback Jitse. 

Its probably been too long ago since I worked on the Eigen internals. You are of course right; eval<T>::type alone is enough and the conditional is not required.

Regarding nested<> I wrote "the opposite" since I thought I remembered that in general nested<> does nest expressions by value and not by PlainObjectType. I want all expressions to be nested as PlainObjectType - no matter how small they are. So in case I am not wrong regarding nested, nesting with eval<T>::type seems to be enough.

Regards,
- Hauke
Comment 9 Hauke Heibel 2011-05-12 13:22:34 UTC
Created attachment 168 [details]
Nest expressions as PlainObjects and everything else by reference.

I updated the patch and the description. To be more clear, expressions shall be nested as PlainObjects (Matrix or Array, by value).
Comment 10 Hauke Heibel 2011-05-12 20:55:45 UTC
Created attachment 169 [details]
Fixed a typo.

Changed MatrixType to OtherDerived.
Comment 11 Bram de Jong 2011-09-12 11:02:25 UTC
I found something related to this, but with /=... I was asked to post the example in here:


// Observed: [2.00342  3.70924  1.16536  3.35758 -1.38856 -2.90473]
// Expected: [1 1 1 1 1 1]

#include "Eigen/Core"

#include <iostream>

using namespace Eigen;

int main(int argc, char* argv[])
{
    MatrixXd m = MatrixXd::Random(5,6);
    // normalize the columns:
    m.array() /= m.colwise().sum().replicate(m.rows(), 1).array();
    // after this, all sum(over column) should evaluate to 1
    std::cout << m.colwise().sum() << std::endl;
}
Comment 12 Bram de Jong 2011-09-12 11:04:53 UTC
Proof of the last comment:

    MatrixXd m = MatrixXd::Random(5,6);
    // normalize the columns:
    MatrixXd tmp = m.colwise().sum();
    m.array() /= tmp.replicate(m.rows(), 1).array();
    // after this, all sum(over column) should evaluate to 1
    std::cout << m.colwise().sum() << std::endl;

Does do the right thing...

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