Re: [eigen] Tiny matrix in Eigen2

[ Thread Index | Date Index | More Archives ]

First, thanks for all the replies, I understand now better Eigen2!

Benoit Jacob wrote:
2009/9/17 WANG Xuewen <>:

We are trying to replace our home made matrix vector library with Eigen2
which provides an unified interface between tiny matrix and medium size
matrix and sounds promising on performance. Surprisingly, it is 10%
slower for tiny matrix (with 2.0.5, not sure about the devel branch
since same code doesn't compile using it), so I wonder if there is
something wrong that I've made and anything that can help improving it.

Our tiny matrix whose dimension limits to 7 but the exact dimension is
not known at compile time. The most costly computation on the matrix is
to compute the exponential of a real or complex matrix. We use the
method similar to MatrixExponential in unsupported.

I've used something like

typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic,
                    Eigen::RowMajor | Eigen::DontAlign,
                    NMaxTinyMatrixDimension, NMaxTinyMatrixDimension>

typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic,
                   Eigen::RowMajor | Eigen::DontAlign,
                   NMaxTinyMatrixDimension, NMaxTinyMatrixDimension>

typedef Eigen::Matrix<double, Eigen::Dynamic, 1,
                    Eigen::RowMajor | Eigen::DontAlign,
                    NMaxTinyMatrixDimension, 1>

typedef Eigen::Matrix<std::complex<double>,
                    Eigen::Dynamic, 1,
                    Eigen::RowMajor | Eigen::DontAlign,
                    NMaxTinyMatrixDimension, 1>

Where NMaxTinyMatrixDimension is 7.

My questions are:

1. Is the chosen storage the right choice for my problem? Should I just
use Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> etc?

That seems indeed what you described: Here you're telling Eigen: I
don't know the exact size of my matrix, but it's not bigger than

Thus Eigen is able to create your matrices on the stack without a
dynamic memory allocation, which is good, but it won't be able to
unroll any loop, because the exact size isn't known at compile time.

If you use instead MatrixXd which is a typedef for
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>, then the matrix
arrays will be dynamically allocated (on the heap).

The only way in which MatrixXd could be faster, is that it allows
vectorization, while with a fixed max-size 7*7 you won't get
vectorization. Perhaps (don't know) for some operations it can start
being faster, but I wouldn't count on it in general.

I have no problem to change NMaxTinyMatrixDimension to a pair (for example to 6 or 8) and change to use Aligned so can possibly get the vectorization. I've done a quick test and it seems it did improve ( I'll do a bit more test to be sure).

Even if the dynamic dimension is not necessarily the same as 8 or 6, for example, if it is 3, will I be able to gain from vectorization?
2. Even if the storage is in the stack, but the dimension is not known
at compile time. Will the loop get unrolled for most operations?

No, no loop will be unrolled at all, since that requires knowing the
exact size at compile time.

Is it possible to extend Eigen2 to allow unroll the loop? We did it for the home made library with meta programming, but I'm not sure if it is possible to do it with Eigen2. I may take a look if you tell me that it is possible and give me hint on how. Or will vectorization be more interesting anyway?
3. We do something like M = N + scalar * Matrix::Identity(); as seen in
unsupported. Is this optimal? Does it really matter?

Ah, this isn't really optimal, because everytime a coefficient is read
in the Identity() expression, that causes an if(). We need to
implement a special optimized code path for that. Meanwhile, rather

M = N;
M.diagonal().cwise() += scalar;

(don't remember, perhaps this cwise()+= requires #include<Eigen/Array>)

This doesn't work with 2.0.5 (doesn't compile, I tried it from the beginning) but apparently works with devel branch so thanks for the hint! If it will get optimized automatically in the future, that would be great!
4. Our internal library stores a complex matrix by using two matrix, one
for the real part and another for the imaginary part. It seems that this
improves but is it still preferred even with Eigen2?

In Eigen, we don't do that. I haven't done benchmarks so I don't know
which approach is faster, perhaps for small sizes yours is better....

I got comparable number now but I'll do more computation.

thanks again,


Mail converted by MHonArc 2.6.19+