Eigen Matrix of complex not entirely happy yet

[orania]

Hello all,

I'm wanting to express my system in terms of Eigen matrices of complex. Apparently, this stretches the current implementation in eigen_algebra.hpp a bit. I've managed to twist my brain around the first hurdle, I think (see first hack in my code example), but I'm afraid I'm stumped on how to proceed. Eigen now responds with a "YOU_MIXED_DIFFERENT_NUMERIC_TYPES" assertion, which, I'm thinking, might actually be telling me that I'm violating Eigen's lazy evaluation semantics somehow (I really don't understand expression templates yet). In any event, if anyone could point me in the right direction, I would be very, very grateful.

#include <iostream>
#include <complex>
#include <Eigen/Core>
#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/external/eigen/eigen.hpp>

// hack #1: not happy with the lpNorm return type given a Matrix of complex<T> (is there a cleaner solution?)
namespace boost { namespace numeric { namespace odeint {
template<typename T,int S1,int S2,int O, int M1, int M2>
struct vector_space_norm_inf< Eigen::Matrix<std::complex<T>,S1,S2,O,M1,M2> > {
    typedef T result_type;
    result_type operator()( const Eigen::Matrix<std::complex<T>,S1,S2,O,M1,M2> &m ) const {
        return m.template lpNorm<Eigen::Infinity>();
    }
};
} } } // end boost::numeric::odeint namespace

//#define _COMPILES_    // comment in/out to switch use cases

using namespace std;
using namespace Eigen;
using namespace boost::numeric::odeint;

#ifdef _COMPILES_
typedef double value_type;
#else // DOESN'T COMPILE
typedef complex<double> value_type;
#endif //

typedef Matrix<value_type, 3, 2> state_type;

void solve(const state_type &x, state_type &dxdt, const double t) {
    dxdt = x;
}

int main() {
    state_type x;
    x <<    1., 2., 3.,
            4., 5., 6.;

    typedef runge_kutta_dopri5<state_type> stepper_type;
    integrate_adaptive(make_dense_output<stepper_type>(1.0e-6, 1.0e-6), solve, x, 0., 10., .01);
    std::cout << x << std::endl;
}

I'm using current git ODEINT, Eigen 3.2.0, GCC 4.8.2. And thank you!

[orania]

Please ignore my chatter if it is not helpful. Without my hack above, default_error_checker::error() complains about returning a complex, rather than a double (since it's presumably ended up with double as value_type, though I'm still trying to understand from where, exactly). While that seems reasonable, I'm beginning to think that it's precisely here that I'm taking the wrong path - at least in terms of wanting to retain the efficiency of Eigen's expression templates. Now I have removed my hack and tried:

typedef runga_kutta_dopri5<state_type, value_type, state_type, double, vector_space_algebra> stepper_type

... but that just gets very ugly. Pity, I was hoping to see value_type find its way into default_error_checker and make everyone play nice, but no such luck as yet. A luta continua!

[mariomulansky]

Hi, an Eigen-Matrix of complex values has indeed not been tested yet. Thanks fro bringing this up here. I think your first approach is the right way to fix it and it looks reasonable. Can you provide a full minimal example so I can test it for myself and maybe see where things go wrong? Using complex as odeint's value_type is not what you want. The value type represents the type of the internal numerical constants of the algorithms and should remain double in your case.

[orania]

Thanks sincerely for your reply. I understand the code I gave in my first post to be a full minimal example (I've included a #define _COMPILES_ switch in order to compare compiling against non-compiling cases). Would you like me to fill out the system as defined in solve()? The system I'm employing is essentially an extension of Stuart-Landau (you have an example of your own of that system here.), but I'm sure a simple dxdt = x.cwiseProduct(alpha) would suffice, with alpha a matrix of complex<double> , same dimensions as x, all elements having non-zero values. Sorry, if that's not clear I can edit my post above. Or, if it's something else that's missing, please let me know. If you'd prefer the full system I'm employing, I'm happy to post that too, but it's big.

[mariomulansky]

ah right. sorry in the email from github i saw only the first part. I will have a look.

[headmyshoulder]

It seems that eigen cannoo construct a complex matrix from a matrix of doubles. The following code fails:

Matrix< complex< value_type > , 3, 2> m1;
Matrix< double , 3 , 2 > m2;
m1 = m2;

And something similar also happens in odeint. I don't know if this can easily be fixed.

[mariomulansky]

hrm ok, i was afraid of something like that... anyways I will try to see where exactly things fail,

[headmyshoulder]

I know where it fails in default_operations::rel_error. I also think I found a solution, wait a second :)

[mariomulansky]

I've tried with the latest version of Eigen from trunk and it compiles with your fix. So it seems Eigen has fixed the problem just recently... Anyhow, karsten is currently looking into a fix in odeint for the Eigen release version, so either you can switch to Eigen trunk or wait a bit more until karsten comes up with a fix in odeint.

[orania]

Thank you Mario, thank you Karsten. What I was suggesting earlier regarding the choice of value_type relates to Karsten's observation above. Eigen is built on lazy evaluation of expression templates (please excuse me using big words that I don't fully understand yet). My first versions (pre-Eigen) iterated through std::vector<complex<double> >, but even after squeezing everything I could out of OpenMP, my execution time was disappointing (my target is real-time processing of digital audio). My first single-process attempts with Eigen matched those I'd achieved with vector<> and OpenMP, so I've been hoping to squeeze even more out of my system that way, and then perhaps a little more with the reintroduction of OpenMP. This is why I'm being pedantic about trying to preserve all possible opportunities for Eigen to do its vectorisation magic (again, I don't understand, but I see it's fast). Having said all of that... I can confirm that I am able to compile my minimal example above with Eigen trunk (sorry for not thinking to try that earlier!), and I am now fervently hacking my way through my full system in order to construct a way to draw fair comparisons between these two implementations. I will report back, but in the meantime: thank you sincerely for your help. It is truly appreciated!

[mariomulansky]

If you really need speed you should also consider to drop complex and formulate the ODE for the real and imaginary component separately, i.e. arriving at 2N-dimensional system. complex does not support expression templates and this might give you severe performance drops. Eigen might have a way around this, but I wouldnt know how to be honest - well except introducing expression templates for complex... Anyways, also make sure to compile with -DNDEBUG to get optimal performance out of eigen (if you use bjam this is automatically set in release mode, though)

[orania]

Thanks Mario, all advice is much appreciated. I been quite far down this road (separating real and imaginary components), and initially the performance gain was astonishing. However, at some point (I'm still not precisely sure why or where), I started seeing negligible difference between my complex and separated implementations, so I switched back. I will heed your advice and be sure to check again. I know nothing of bjam, so I'll look into that. I don't honestly know that I'll be able to attain real-time, but at least if I give it my best shot, then time spent waiting for experimental results might come down from weeks or months to hours or days. Right now I'm doing 1 sec of audio in anything between 5 and 300 sec (depending on various options I'm toying with). Once again, thank you sincerely.

[orania]

Just a quick report-back: I've replaced my earlier "iteration over vector<> implementation (with OpenMP pragmas) with Eigen matrices - in both cases wrapping complex<double>s, and I'm happy that 1) it compiles, and 2) that the performance appears very similar, even though I'm invoking no explicit OpenMP pragmas in the second case. I will be tinkering to squeeze even more out of this approach, and post my conclusions as they emerge. My sincere thanks to the odeint team! I consider this issue closed.