Suppose I am receiving chunks of m doubles from an external sensor (think of this being a single observation or sample). These m doubles represent a single row.

I would like to do an svd() on the last n samples of data every time the data arrives, which is basically an n by m matrix. Some things I am having a problem with:

(1) I know the total number of samples that will arrive. I am thinking of storing them all in a giant Eigen matrix, and doing svd() on minors, to mimic a sliding window. Is there a better way of doing this?
(2) The data arrives in a form of a double pointer; is there an efficient way of copying this data into a single row (or column) of an Eigen matrix, something like memcpy()?
(3) I would like to be able to call some LAPACK routines on that data as well, I guess this means I would need access to the underlying column-major storage?

Thanks.

to trade the memory consumption/memory copies you could have a (m+k)*n matrix, at each step you would add a new row:

mat.row(i+m) = RowVectorXd::Map(ptr,n);
++i;

compute the svd of:

mat.middleRows(i,m);

call lapack routines using:

mat.middleRows(i,m).data()
mat.outerStride() (== leading dim)

and of course before inserting a new row:

if(i==k) {
mat.topRows(m) = mat.middleRows(i,m);
i = 0;
}

Thank you!

How efficient is the insertion
mat.row(i+m) = RowVectorXd::Map(ptr,n);
?

My guess is this would depend on the internal data storage in Eigen.
If the storage is column-major, the above operation would
insert n numbers in n different places of memory. In this case
I should think of incoming data as separate columns in the matrix,
so that this assignment inserts a contiguous chunk of data.
Thanks a lot for your suggestions.

you can still use a row major matrix to make this insertion efficient but 1) you won't any perf gain because the SVD is much more mostly anyway, and 2) this will make the calls to Lapack routines more difficult .

This makes sense. Thank you.

Regarding the operation

mat.topRows(m) = mat.middleRows(i,m);

is it done "in-place", or is there any memory allocation going on behind the scenes?

Thanks.

it's done inplace, so the reverse operation:

mat.middleRows(i,m) = mat.topRows(m);

would not work as expected if k<m but in your case it's ok.

I decided to represent each observation as a column in a matrix, so my matrices are now MxN where N is the number of M-dimensional samples. I will have to transform this matrix anyway before applying svd(), so the storage format won't matter.

But here is a puzzle. I wrote a code that repeatedly inserts a column into a matrix, and wraps around when needed. There are two functions, one keeps the data as an Eigen matrix, the other one mimics the same behavior with a double pointer. The function that uses Eigen::MatrixXd is considerably slower. That function pretty much follows the code that you suggested. I am pasting the code (compiled with g++ on 64 bit linux). I am a newbie with Eigen, and you'll probably realize that I have no clue how to make it fast; hoping for some suggestions.

The #define NO_EIGEN switches between the function using double and Eigen code. The Eigen header might have a different path on your system.

Thank you.

Code: Select all

#include <eigen3/Eigen/Eigen>
#include <iostream>
#include <iomanip>
#include <stdlib.h>
#include <sys/time.h>
#include <string.h>
#include <sstream>

using namespace Eigen;
using namespace std;

const long rows = 200;
const long cols = 500;
const long actual_cols = 1500;

const long ntimes = 5;
const long ninsertions = 45000;

#define NO_EIGEN

std::string toString(double *p) {
    std::stringstream ss;
    ss << "Full matrix:\n";
    for (long i = 0; i < rows; ++i) {
        for (long j = 0; j < actual_cols; ++j) {
            ss << p[j*rows + i] << " ";
        }
        ss << '\n';
    }
    return ss.str();
}

EIGEN_DONT_INLINE void  doublefunc(double *p, double *r) {
    double *ptr_p = p + rows * cols;
    long k = 0;
    for (long i = 0; i < ninsertions; ++i) {
        for (long m = 0; m < rows; ++m) {
            r[m] = ++k;
        }
        if (actual_cols == (ptr_p - p)/rows) {
            memcpy(p, ptr_p-rows*(cols-1), rows*(cols-1)*sizeof(double));
            ptr_p = p + rows * (cols-1);
        }
        memcpy(ptr_p, r, sizeof(double) * rows);
        //std::cerr << toString(p) << std::endl;
        ptr_p += rows;
    }
    return;
}

EIGEN_DONT_INLINE void  eigenfunc(double *p, double *r) {
    int current_col = 0;
    int k = 0;
    Eigen::MatrixXd mat(rows, actual_cols);
    for (long i = 0; i < ninsertions; ++i) {
        for (long m = 0; m < rows; ++m) {
            r[m] = ++k;
        }
        if (current_col == actual_cols - cols) {
            mat.leftCols(cols) = mat.rightCols(cols);
            current_col = 0;
        }
        mat.col(current_col + cols) = VectorXd::Map(r,rows);
        //std::cerr << mat << std::endl << std::endl;
        ++current_col;
    }
}

int main(void) {
    EIGEN_ALIGN16 double *p = new double[rows * actual_cols];
    EIGEN_ALIGN16 double r[rows];
    memset(p, 0, rows * actual_cols * sizeof(double));

    struct timeval tv;
    double total_time = 0.;

    for (long i = 0; i < ntimes; ++i) {
        gettimeofday(&tv, 0);
        double start_time = tv.tv_sec + tv.tv_usec * 1.0e-6;
#ifdef NO_EIGEN
        doublefunc(p, r);
#else
        eigenfunc(p, r);
#endif
        gettimeofday(&tv, 0);
        double end_time = tv.tv_sec + tv.tv_usec * 1.0e-6;
        double time_spent_microsec = end_time - start_time;
        total_time += time_spent_microsec;
    }

    std::cerr << std::setprecision(9) << "total time spent: " << total_time << std::endl;

}



Here I observe the same level of performance. Make sure you compile with optimizations.

Thanks you, indeed I was not using optimization, turning it on was a huge improvement.

Code: Select all

          DOUBLE           EIGEN   
NO OPTIM   1.50157           10.75500   (average seconds per run)
-O            1.06710   1.40339   
-O2         1.06269       1.27520   
-O3         1.05312           1.17498