This seems like it must have been asked before (apologies, but I did spend quite a bit of time searching).

I have a matrix which I'm really using just a collection of column vectors. I'd like to have fast access to the coefficients of these columns. With pointers I'd do something like this:

// xp is a float pointer to the array of values, stored in column-major order
columnIndex = 27;
float* colp = xp + columnIndex * nRows; // points to beginning of desired column
float value1 = colp[7] + colp[3]*colp[5];
float value2 = colp[4] + colp[8];

This seems superior to a naive way of doing it in Eigen,
float value = X(7,columnIndex) + X(3,columnIndex)*X(5,columnIndex);
because (if my debugger is to be believed) each coefficient access involves a multiplication by the equivalent of nRows.

In Eigen I feel like I should be able to do this:
typedef Eigen::Matrix<float,10,Eigen::Dynamic> MyMatrix;
MyMatrix X(10,50);
MyMatrix::ColXpr c; // here's the problematic line
c = X.col(columnIndex);
float value1 = c(7) + c(3)*c(5);
float value2 = c(4) + c(8);
but I get compiler errors ("no matching function call...").

I know I can copy the column to a temporary, but I'd rather avoid copies if at all possible.

I'm sure I'm missing something obvious.

you have to use the ctor, e.g.:

MyMatrix::ColXpr c = X.xol(columnIndex);

or:

MyMatrix::ColXpr c(X.xol(columnIndex));

or:

MyMatrix::ColXpr c(X,columnIndex);

Also, note that the line:

c = X.col(columnIndex);

means copy each element of X.col(columnIndex) into the vector referenced by c, that is clearly not what you want here!

Many thanks! For the record (i.e., others who read this), "xol" should be "col", and here is a fully-working example of the kind of thing I was aiming for (uses placement new to reassign the column "pointer" each time):

typedef Eigen::Matrix<int,2,Eigen::Dynamic> MyMatrix;

MyMatrix m(2,3);
m <<
1,2,3,
4,5,6;

MyMatrix::ColXpr c = m.col(0);
for (int i = 0; i < 3; i++) {
new (&c) MyMatrix::ColXpr(m.col(i));
cout << "Column " << i << ":\n" << c << endl;
}

One remaining issue (low priority) is having a ColXpr as a member of a class, like this:
template <typename T,int n_dims>
class MatrixWithFastColumnAccess {
typedef Eigen::Matrix<T,n_dims,Eigen::Dynamic> MyMatrix;

private:
const MyMatrix m;
typename MyMatrix::ColXpr curCol;

public:
MatrixWithFastColumnAccess(const MyMatrix& _m) : m(_m) {;}
void set_column(int i) {
new (&curCol) typename MyMatrix::ColXpr(m.col(i));
}
T operator[](int index) const {return curCol(index);}
};

At compilation, the construction of the class causes an error (again, "no matching function call..."). I've tried changing the constructor to the following,

MatrixWithFastColumnAccess(const MyMatrix& _m) : m(_m), curCol(_m.col(0)) {;}

but that doesn't fix the problem.

this is a constness issue, use XXX::ConstColXpr

Ah. Thanks a lot!

OK, it seemed like I was almost there, and by tackling other problems I have learned a lot about Eigen/C++ in the last couple of days, but this particular one is still bugging me. Sorry if I'm being a total idiot.

I'm trying to get this program to compile:

Code: Select all

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

using namespace std;

template <typename T,int n_dims>
class MatrixWithFastColumnAccess {
public:
  typedef Eigen::Matrix<T,n_dims,Eigen::Dynamic> MyMatrix;
  typedef typename MyMatrix::ConstColXpr ConstCol;

  MatrixWithFastColumnAccess(const MyMatrix& X) : mrX(X), curCol(X.col(0)) {;}

  void set_column(int i) {
    new (&curCol) ConstCol(mrX.col(i));
  }
  T operator[](int index) const {return curCol(index);}

private:
  const MyMatrix&  mrX;
  ConstCol         curCol;
};

int main()
{
  const int n_dims = 2;
  const int n_points = 100;
  typedef Eigen::Matrix<float,n_dims,Eigen::Dynamic> MatrixType;

  MatrixType X = MatrixType::Random(n_dims,n_points);
  MatrixWithFastColumnAccess<float,n_dims> XC(X);

  // Now exercise these class elements lots of times so we can use callgrind to measure performance
  float val = 0;
  for (int i = 0; i < 1000000; i++)
    for (int j = 0; j < n_points; j++) {
      XC.set_column(j);
      int rowIndex = i % n_dims;
      val += XC[rowIndex];
    }
  cout << val << endl;
} 


However, I get this error:

> g++ test07.cpp -o test07
test07.cpp: In member function ‘void MatrixWithFastColumnAccess<T, n_dims>::set_column(int) [with T = float, int n_dims = 2]’:
test07.cpp:36: instantiated from here
test07.cpp:15: error: invalid conversion from ‘const void*’ to ‘void*’
test07.cpp:15: error: initializing argument 2 of ‘void* operator new(size_t, void*)’

(Line 15 is the instantiation of set_column.) I am having a hard time figuring out why it is trying to convert a const void* to a void* when everything is defined as const. I also don't see where I could plausibly add a const qualifier to fix things.

Now, the following alternative version of the class definition works fine:

Code: Select all

template <typename T,int n_dims>
class MatrixWithFastColumnAccess {
public:
  typedef Eigen::Matrix<T,n_dims,Eigen::Dynamic>        MyMatrix;
  typedef Eigen::Map<const Eigen::Matrix<T,n_dims,1> >  MyColumn;

  MatrixWithFastColumnAccess(const MyMatrix& X) : pX(&X(0,0)), d(X.rows()), curCol(pX,d) {;}

  void set_column(int i) {
    new (&curCol) MyColumn(pX+i*d,d);
  }
  T operator[](int index) const {return curCol(index);}

private:
  const T* pX;
  int d;
  MyColumn curCol;
}


Are there any disadvantages here? It seems that possible advantages, like the memory-alignment of columns, might be "thrown away" compared with the first version, but maybe this is all an in-principle problem rather than a real problem.

oh I see, this is because ConstCol is const at two levels, it would be like:

const Eigen::Map<const Eigen::Matrix<T,n_dims,1> >

in your second version. The solution is to remove the outside const qualifier:

typename remove_const<ConstCol>::type curCol;

Thanks, that works!

For anyone else reading this, I implemented remove_const like this:

Code: Select all

template <typename T>
struct remove_const { typedef T type; };

template <typename T>
struct remove_const<const T> { typedef T type; };


Much to my suprise, in profiling the results, the Map version was actually faster by almost a factor of 2. The issue seems to come down to the large number of calls to

Code: Select all

float* Eigen::internal::const_cast_ptr<float>(float const*)


Map was also faster if I used Unaligned rather than Aligned.

I am puzzled about why "discarding" so much information about where the column comes from is a disadvantage, but at any rate I now have a much better appreciation of the various options. Thanks again for your help!

hm, it seems you are benchmarking in debug mode, i.e., without optimizations because float* Eigen::internal::const_cast_ptr<float>(float const*) is no-op. Benchmarking Eigen code without optimizations enabled does not make any sense because of the large overhead due to the abstraction.