Eigen efficient type for dense symmetric matrix

问题:

Does Eigen have efficient type for store dense, fixed-size, symmetric matrix? (hey, they are ubiquitous!)

I.e. for N=9, it should store only (1+9)*9/2==45 elements and it has appropriate operations. For instance there should be efficient addition of two symmetric matrices, which returns simmilar symmetric matrix.

If there is no such thing, which actions (looks like this) I should make to introduce such type to Eigen? Does it has concepts of "Views"? Can I write something like "matrix view" for my own type, which would make it Eigen-friednly?

P.S. Probably I can treat plain array as 1xN matrix using map, and do operations on it. But it is not the cleanest solution.

回答1:

Yes, eigen3 has the concept of views. It doesn't do anything to the storage though. Just as an idea though, you might be able to share a larger block for two symmetric matrices of the same type:

Matrix<float,4,4> A1, A2; // assume A1 and A2 to be symmetric
Matrix<float,5,4> A;
A.topRightCorner<4,4>().triangularView<Upper>() = A1;
A.bottomLeftCorner<4,4>().triangularView<Lower>() = A2;

Its pretty cumbersome though, and I would only use it if your memory is really precious.

回答2:

Packed storage of symmetric matrices is a big enemy of vectorized code, i.e. of speed. Standard practice is to store the relevant N*(N+1)/2 coefficients in the upper or lower triangular part of a full dense NxN matrix and leave the remaining (N-1)*N/2 unreferenced. All operations on the symmetric matrix are then defined by taking into account this peculiar storage. In eigen you have the concept of triangular and self-adjoint views for obtaining this.

From the eigen reference: (for real matrices selfadjoint==symmetric).

Just as for triangular matrix, you can reference any triangular part of a square matrix to see it as a selfadjoint matrix and perform special and optimized operations. Again the opposite triangular part is never referenced and can be used to store other information.

Unless memory is a big problem, I would suggest to leave the unreferenced part of the matrix empty. (More readable code, no performance problems.)

回答3:

Efficient type for symmetric matrix

You just assign values to the lower/upper triangular parts of the matrix and use the Eigen triangular and selfadjoint views. However, I have tested both on small fixed-size matrices. I noticed that performance-wise, using views is not always the best choice. Consider the following code:

Eigen::Matrix2d m;
m(0,0) = 2.0;
m(1,0) = 1.0;
// m(0,1) = 1.0;
m(1,1) = 2.0;
Eigen::Vector2d v;
v << 1.0,1.0;
auto result = m.selfadjointView<Eigen::Lower>()*v;

The product in the last line is quite slow compared with the alternative solutions presented below (about 20% slower for double 2x2 matrices in my case). (The product using the full matrix, by uncommenting m(0,1) = 1.0;, and using auto result = m*v, is even faster for double 2x2 matrices).

EDIT: I forgot about aliasing. auto result.noalise() speeds things up (the alternatives below are still faster though).

Some alternatives.

1) Store symmetric matrix in vector

You could store your matrix in a vector of size 45. Summing 2 matrices in vector format is straightforward (just sum the vectors). But you have to write your own implementation for products.

Here is the implementation of such a matrix * vector product (dense, fixed-size) where the lower part of the matrix is stored column-wise in a vector:

template <typename T, size_t S>
Eigen::Matrix<T,S,1> matrixVectorTimesVector(const Eigen::Matrix<T,S*(S+1)/2,1>& m, const Eigen::Matrix<T,S,1>& v)
{
    Eigen::Matrix<T,S,1> ret(Eigen::Matrix<T,S,1>::Zero());
    int counter(0);
    for (int i=0; i<S; ++i)
    {
        ret[i] += m(counter++)*v(i);
        for (int j=i+1; j<S; ++j)
        {
            ret[i] += m(counter)*v(j);
            ret[j] += m(counter++)*v(i);
        }
    }
    return ret;
}

2) Store only the triangular part and implement your own operations

You could of course also implement your own product matrix * vector, where the matrix only stores the 45 elements (let's assume we store the lower triangular part). This would maybe be the most elegant solution, as it keeps the format of a matrix (instead of using a vector which represents a matrix). You can then also use Eigen functions like in the example below:

template <typename T, size_t S>
Eigen::Matrix<T,S,S> symmMatrixPlusSymmMatrix( Eigen::Matrix<T,S,S>& m1, const Eigen::Matrix<T,S,S>& m2)
{
    Eigen::Matrix<T,S,S> ret;
    ret.template triangularView<Eigen::Lower>() = m1 + m2; // no performance gap here!
    return ret;
}

In the function above (sum of 2 symmetric matrices), only the lower triangular parts of m1 and m2 are visited. Note that the triangularView does not present a performance gap in this case (I affirm this based on my benchmarks).

As for the matrix * vector product, see the example below (same performance as the product in alternative 1)). The algorithm only visits the lower triangular part of the matrix.

template <typename T, size_t S>
Eigen::Matrix<T,S,1> symmMatrixTimesVector(const Eigen::Matrix<T,S,S>& m, const Eigen::Matrix<T,S,1>& v)
{
    Eigen::Matrix<T,S,1> ret(Eigen::Matrix<T,S,1>::Zero());
    int counter(0);

    for (int c=0; c<S; ++c)
    {
        ret(c) += m(c,c)*v(c);
        for (int r=c+1; r<S; ++r)
        {
            ret(c) += m(r,c)*v(r);
            ret(r) += m(r,c)*v(c);
        }
    }
    return ret;
}

Performance gain for the product Matrix2d*Vector2d when compared to the product using the full matrix (2x2 = 4 elements) is 10% in my case.