templatesc++14eigencomplex-numbersceres-solver

Eigen Matrix Assignment Operator with User-Defined Data Types


I am using the Eigen matrix library to deal with matrices of std::complex<T> data types, where T is either of type double or of type ceres::Jet<double,...>. Eigen documentation indicates that << is the correct operator to use for assignment, but it seems that << is not overloaded for matrices of user-defined data types. Is there a different method I can use to initialize Eigen matrices that works for both data types?


Solution

  • The problem with your code is not that Eigen is not overloaded for matrices of user-defined type (it is) but that you have a nested template argument of the type Eigen::Matrix<std::complex<T>, Eigen::Dynamic, Eigen::Dynamic> (where T is something like ceres::Jet<double> resulting in Eigen::Matrix<std::complex<ceres::Jet<double>>, Eigen::Dynamic, Eigen::Dynamic>). You will have to explicitly construct the elements of the nested type T before using the streaming operator <<.

    Eigen::Matrix<std::complex<T>, 3, 1> mat;
    mat << T{1.0}, T{2.0}, T{3.0};
    

    More detailed explanation

    To further explain this let's start with a simple class we will use for the element type

    class SomeClass {
      public:
        constexpr SomeClass(double const& some_data = 0.0) noexcept
          : some_data{some_data} {
          return;
        }
    
      private:
        double some_data;
    };
    

    A matrix of this user-defined type can be successfully assigned with:

    Eigen::Matrix<SomeClass, 3, 1> mat;
    mat << 1.0, 2.0, 3.0;
    

    The problem only emerges if you apply it to nested element types such as std::complex<SomeClass>. We would like to use the following overload:

    CommaInitializer<Derived> Eigen::DenseBase<Derived>::operator<<(const Scalar& s)
    

    where Scalar is the type of the matrix coefficients (element type, e.g. Derived = Eigen::Matrix<std::complex<T>, 3, 1> and Scalar = std::complex<T>):

    Eigen::DenseBase<Derived>::Scalar
    

    Therefore when compiling GCC will tell you that it can't match the template:

    note: candidate: Eigen::CommaInitializer<Derived> Eigen::DenseBase<Derived>::operator<<(const Scalar&) [with Derived = Eigen::Matrix<std::complex<SomeClass>, 3, 1>; Eigen::DenseBase<Derived>::Scalar = std::complex<SomeClass>
    inline CommaInitializer<Derived> DenseBase<Derived>::operator<< (const Scalar& s)
    
    note: no known conversion for argument 1 from ‘double’ to ‘const Scalar& {aka const std::complex<SomeClass>&}’
    

    What you will have to do to make it compile is to construct the elements with a type which can be converted to Scalar and can then call the corresponding streaming operator as follows:

    mat << std::complex<SomeClass>{SomeClass{1.0}}, std::complex<SomeClass>{SomeClass{2.0}}, std::complex<SomeClass>{SomeClass{3.0}};
    
    mat << std::complex<SomeClass>{1.0}, std::complex<SomeClass>{2.0}, std::complex<SomeClass>{3.0};
    
    mat << SomeClass{1.0}, SomeClass{2.0}, SomeClass{3.0};
    

    This should be avoidable by changing the Eigen library to something like

    CommaInitializer<Derived> Eigen::DenseBase<Derived>::operator<<(Elem const& s)
    

    and require with std::enable_if_t, static_assert or C++20 concepts that Scalar can be constructed by an element of type Elem:

    std::is_constructible_v<Scalar, Elem>
    

    If you really need you could try to write an overload for it yourself that makes sure that std::is_constructible_v<Scalar, Elem> but !std::is_same_v<Scalar, Elem> but I personally think it is never a good idea to add such a functionality to an existing library yourself. Somebody else copying fragments of your code and expecting them to work might end up with not working code.


    Alternative options for assignments

    As alternatives to assignment with << you could

    For the first two options you will have to use the same logic as discussed in the previous paragraph or use double braces

    Eigen::Matrix<std::complex<SomeClass>, 3, 1> mat = { {1.0, 2.0}, {2.0, 3.0}, {3.0, 4.0} };
    

    or

    std::vector<std::complex<SomeClass>> vec = { {1.0, 2.0}, {2.0, 3.0}, {3.0, 4.0} };
    Eigen::Matrix<std::complex<SomeClass>, 3, 1> mat{vec.data()};
    

    while something like

    mat(0,0) = 1.0;
    
    mat(0,0) = {1.0, 2.0};
    

    will create a complex element with a real part 1.0 and imaginary part 0.0 or 2.0 of type SomeClass without having to call the constructor SomeClass{} explicitly.