Search code examples
c++eigeneigen3

Vectorize a Symmetric Matrix


I would like to write a function with the following signature

VectorXd vectorize (const MatrixXd&);

which returns the contents of a symmetric matrix in VectorXd form, without repeated elements. For example,

int n = 3;    // n may be much larger in practice.
MatrixXd sym(n, n);
sym << 9, 2, 3,
       2, 8, 4, 
       3, 4, 7;

std::cout << vectorize(sym) << std::endl;

should return:

9
2
3
8
4
7

The order of elements within vec is not important, provided it is systematic. What is important for my purposes is to return the data of sym without the repeated elements, because sym is always assumed to be symmetric. That is, I want to return the elements of the upper or lower triangular "view" of sym in VectorXd form.

I have naively implemented vectorize with nested for loops, but this function may be called very often within my program (over 1 million times). My question is thus: what is the most computationally efficient way to write vectorize? I was hoping to use Eigen's triangularView, but I do not see how.

Thank you in advance.


Solution

  • Regarding efficiency, you could write a single for loop with column-wise (and thus vectorized) copies:

    VectorXd res(mat.rows()*(mat.cols()+1)/2);
    Index size = mat.rows();
    Index offset = 0;
    for(Index j=0; j<mat.cols(); ++j) {
        res.segment(offset,size) = mat.col(j).tail(size);
        offset += size;
        size--;
    }
    

    In practice, I expect that the compiler already fully vectorized your nested loop, and thus speed should be roughly the same.