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).
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.
i,jelements. It is pretty simple in fact to unfold a double for loop over two indicesi=1...N, j=i...N, into a single loop. – Discomfortable