9

Let S be a symmetric n x n matrix and A be a m x n matrix.

Given: B = A * S * A_transpose (where "*" represents a matrix product operation)

B will also be a symmetric matrix.

Using the tuxfamily Eigen library, version 3, what is a clean and efficient way to implement this computation? (By efficient, I mostly mean that duplicate computations of the elements of B are not performed where symmetry makes them unnecessary.)

I'm guessing it will make use of SelfAdjointView, but I have searched high and low and not found a clean example of this.

The application is a Kalman filter, which depends heavily on operations involving (symmetric) covariance matrices, so I want to be sure get the implementation/design correct.

Thank you!

1 Answer 1

6

This should be quite simple. As you say yourself, you can make Eigen aware of the fact that your matrix is a symmetric matrix through the SelfAdjointView. There is also another view, which is the TriangularView, which you can use to store your results. According to the reference if you assign to a TriangularView, only the relevant parts of the rhs are evaluated. So

B.triangularView<Upper>() = A * S.selfadjointView<Upper>() * A.transpose();

will store the result in the upper triangle of B. You can then use B.selfadjointView<Upper> in any further calculations. I am not sure if this is optimal in terms of required operations and you might do some benchmarking to verify.

2
  • Thank you! That makes sense. I had a feeling it was something simple like that, but was having trouble pulling the pieces together. Looking at the reference it appears to me what B.selfAdjointView would accomplish this as well. It seems that for lhs, selfAdjointView and triangularView would behave exactly the same... but I'm guessing, as I'm still figuring this out.
    – arr_sea
    Nov 5, 2012 at 0:50
  • 6
    By the way, my benchmarking efforts (on an older intel i5 chip) suggested that it is significantly more efficient to perform full matrix multiplication, rather than using the triangular and selfadjoint views. I tested with matrix S ranging from 4x4 to 7x7 and matrix A having from 1 to 4 rows. I'm guessing vectorization is what's making the difference in performance. If A is quite sparse, doing the multiplication "by hand" is faster (and using a sparse matrix representation is slower by a factor of 10).
    – arr_sea
    Nov 27, 2013 at 20:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.