Re: [eigen] Tiny matrix in Eigen2

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2009/9/17 WANG Xuewen <xuewen.wang@xxxxxxxxx>:
> Hi,
>
> 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>
>       TinyRMatrix;
>
> typedef Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic,
>                    Eigen::RowMajor | Eigen::DontAlign,
>                    NMaxTinyMatrixDimension, NMaxTinyMatrixDimension>
>       TinyCMatrix;
>
> typedef Eigen::Matrix<double, Eigen::Dynamic, 1,
>                     Eigen::RowMajor | Eigen::DontAlign,
>                     NMaxTinyMatrixDimension, 1>
>       TinyRVector;
>
> typedef Eigen::Matrix<std::complex<double>,
>                     Eigen::Dynamic, 1,
>                     Eigen::RowMajor | Eigen::DontAlign,
>                     NMaxTinyMatrixDimension, 1>
>       TinyCVector;
>
> 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
NMaxTinyMatrixDimension.

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.

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

>
> 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
do:

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

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


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

Benoit