It should be:

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

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();

(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

yes, that's what I meant. This is just a matter of using nested<>.

(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

Created attachment 167 [details]
Nest expressions by value and everything else by reference.

(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?

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

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).

Created attachment 169 [details]
Fixed a typo.

Changed MatrixType to OtherDerived.

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;
}

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...

-- 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/259.