[This is a continuation of the discussion on GitHub in a PR. It made sense to start discussing it on a forum instead of a closed PR. Start of the discussion: stan-dev/math#2096]
@bbbales2 and I were talking about symmetric matrices in the Stan Math library.
There is an absolute tolerance in the
check_symmetric() function. We don’t use
== equality to check for symmetry (which wouldn’t make sense for IEEE floating point numbers anyway). @bbbales2 suggested we use relative tolerance and I wholeheartedly agree. Note: we don’t have a symmetric matrix type. This is just an Eigen matrix that has the property of being close enough to symmetric.
I mentioned that because numerically it’s not truly symmetric, we may end up with operations that result in non-symmetric output even if the linear algebra dictates that it should be symmetric. And this may fail if the output of one function is being used as an input to another function. I had some ideas on how to deal with this and wanted to just write some of them down. (These are just ideas… no prototypes or anything associated with it yet.)
- Create a symmetric matrix type in Stan. We could either use a different C++ type or use Eigen’s triangular views. This would be a bit of work.
- Rewrite functions to use Eigen’s triangular views after checking for symmetry. I’m not sure how this will interact with the autodiff stack… (do we need to copy
stan::math::varto the off-diagonal right before returns?)