Re: [eigen] Eigen is slow with complex matrices

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

Hi,

2009/10/19 Robert Lupton the Good <rhl@xxxxxxxxxxxxxxxxxxx>:
> We're using eigen on a large astronomical project;  it's great -- no need to
> find the proper combination du jour of optimised blas and lapack.
>
> Unfortunately, in the context of complex matrices we are having to
> reconsider the use of BLAS. The reason appears to be in the complex class
> itself (see below)  Is there any way that Eigen can work around this?

There are several ways in which Eigen is currently not optimally fast
with complex numbers, some are Eigen's fault, some are the GNU STL's
fault as you mention. Perhaps in the particular case of your project,
it reduces to one precise reason;

>> A better implementation would be to have op* directly calculate both
>> components of the result, and return it as a new complex object:
>> {
>>   return complex<_Tp>(
>>       _x.real() * _y.real() - _x.imag() * _y.imag(),
>>       _x.real() * _y.imag() + _x.imag() * _y.real());
>> }

Have you tried it? Does it fix your performance issues? If yes, you
have a pretty clear path from here: just copy this standard complex
implementation as a new class, say "ourcomplex" and make that change
there; Then add support for "ourcomplex" as explained here,
http://eigen.tuxfamily.org/dox/CustomizingEigen.html#CustomScalarType

Otherwise, here are some thoughts on the performance issues with complex.

A. Not our fault
0) What you mention --- I didn't know about it.
1) The basic operations are not vectorized by GCC as they should. For
example, an addition of two complex<double> should be trivial to
vectorize.
2) Because of that, we haven't added vectorization in Eigen (though
for complex<float> we should be able to do something)
3) So we've been considering adding our own complex class, see the
TODO. Needs someone with time to do it.

B. Our fault
4) Eigen should be modified to avoid taking abs() of complex numbers
too often, because that means a sqrt(). The idea would be to introduce
ei_abs_or_abs2 that does abs for real and abs-squared for complex, so
that no sqrt() is needed.

To check if 4) is affecting you, edit Core/MathFunctions.h and put
EIGEN_DONT_INLINE in front of the ei_abs() specializations for
complex. Then profile your program, is ei_abs() accounting for a large
part of the execution time?

Cheers,
Benoit