N-dimensional Matrices in C++

I’m working on coding up an algorithm in C++, based on a Matlab implementation I wrote, and I’ve been using the Eigen library (http://eigen.tuxfamily.org/) for my linear algebra. It’s been relatively straightforward to use so far, but unfortunately it doesn’t (as far as I can tell) have any built in support for N-dimensional matrices. In fact, it appears that N-dimensional matrices aren’t really addressed by any of the standard C++ linear algebra libraries. So, to get my 4-dimensional matrix, I have to make an array of pointers pointing to an array of pointers pointing to my Eigen MatrixXf objects. I initialize my 4-D matrices using this function:

MatrixXf*** initialize_4dMatrix(MatrixXf*** matrix, int sizes[], float value=0.0) {
    matrix = new MatrixXf** [sizes[0]];
    for (int d1 = 0; d1 < sizes[0]; d1++) {
        matrix[d1] = new MatrixXf* [sizes[1]];
        for (int d2 = 0; d2 < sizes[1]; d2++) 
            matrix[d1][d2] = new (MatrixXf)(MatrixXf::Ones(sizes[2],sizes[3]).array() * value);
    return matrix;

and call it with:

int sizes[4] = {1,2,3,4};
MatrixXf ***m = initialize_4dMatrix(m, {1,2,3,4});

To then access the matrices, it’s (*m[d1][d2])(d3,d4); which isn’t bad for accessing one matrix at a time, but what if I want the jth element of the ith row of every matrix (i.e. element (i,j))? The Python access method m[:,:,i,:] (or just m[:,:,i]) would be ideal. Unfortunately, I don’t think there’s any way around writing another function to do it:

    static VectorXf* get_elements_across_matrix_array(MatrixXf** matrix, int size, int element[]) { 
        VectorXf* vector = new VectorXf(size);
        for (int i = 0; i < size; i++) 
            (*vector)[i] = (*matrix[i])(element[0],element[1]);
        return vector;

And I would need to write another function that returns a MatrixXf pointer if I wanted to access the ith rows of each matrix.

To delete these matrices I also need to pass into my delete function the sizes of the array. It’s a pain that it can’t dynamically discover the array size, but I don’t think there’s any way around it in C++.

static void delete_4dMatrix(MatrixXf*** matrix, int sizes[]) {
    for (int d1 = 0; d1 < sizes[0]; d1++) {
        for (int d2 = 0; d2 < sizes[1]; d2++)
            delete matrix[d1][d2];

Additionally, I need to have a separate function for initializing my 3D matrices, and any other dimension that I would need the way I have it set up now. Less than ideal. My guess is that there’s a elegant way to handle this with void pointers, but at the moment I only need 3D and 4D matrices so it’s not worth the effort for me to set it up.

So there are a few obstacles to using N-dimensional matrices with Eigen, but it’s really not a problem once you get a handle on it and how to properly reference different elements and matrices through your (N-2)-dimensional array of pointers.

Tagged , ,

6 thoughts on “N-dimensional Matrices in C++

  1. Aaron Schurger says:

    Hi Travis,
    Just wanted to thank you for this posting, which has been enormously helpful to me as I am just getting to know Eigen and odeint, which I am using to set up and solve ODEs.
    Many thanks!

    • travisdewolf says:

      Hi Aaron,

      thanks for the kind words! Really glad that you found this useful! đŸ˜€

      • Aaron Schurger says:

        You’re welcome! I’ve made some progress thanks to your posting, now I am just a bit stuck trying to get Eigen and odeint to work together. I am assuming that you can use an Eigen MatrixXd as the “state_type” for odeint (as below). But that doesn’t seem to work. Have you used Eigen and ideint together? Does this look like the proper way to do it?

        typedef Eigen::Matrix state_type;

        void odefun( const state_type &u , state_type &dudt , double t )


        void write_odefun( const state_type &x , const double t )
        int i;
        cout << t << '\n';
        for (i=0; i<=9; i++) {
        cout << x(i) << '\n';

        int main(int argc, char **argv)

        state_type u(3,3); // initial conditions
        srand((unsigned int) time(0));
        u.setRandom(); // picks random numbers from -1 to 1
        u*=0.1; // make the random numbers smallish

        integrate( odefun , u , 0.0 , 100.0 , 0.1 , write_odefun );


      • travisdewolf says:

        Hmm. I’m unfortunately a bit swamped at the moment, and don’t have time to really dive in, but it seems like there may have been similar problems from others: http://stackoverflow.com/questions/22782224/using-boostodeint-with-eigenmatrix-as-state-vector
        if those don’t do it for you I can try to take a look asap but it might be a bit! Please let me know if you get it! Good luck!!

  2. Aaron Schurger says:

    Thanks for the links! The first one seems highly relevant. I’ll mess around with it a bit and see if I can get it working. If I get really desperate, then I may take you up on your offer! I’m not in a huge hurry. Thanks again!

  3. Aaron Schurger says:

    Everything compiles and runs OK, as long as I comment out the critical line…

    integrate( dudt_fun , u , 0.0 , 100.0 , 0.1 , write_dudt );

    So essentially once I try to use the ode solving capabilities of odeint, I get a compile error:

    error: no matching function for call to ‘integrate(void (&)(state_type&, state_type&, double), state_type&, double, double, double, void (&)(const state_type&, double))’

    If you have any ideas, I’d be very grateful!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: