preferences for symmetrizing matrices (e.g., cov_matrix inverse transform)?

[Bob Carpenter]

Three questions.

Andrew filed the issue

https://github.com/stan-dev/math/issues/342

because he found cov_matrix parameters that were asymmetric
and failing when plugged into multi_normal.

QUESTION I: Does anyone have a preference as to how I do
the symmetrization? The two obvious options are

(1) copying below-diagonal elements to above-diagonal

(2) adding matrix to its own transpose and dividing by 2

Both versions pass existing unit tests and provide answers matching
finite diffs when I run

parameters {
real<lower=0> a;
real b;
cov_matrix[10] Sigma;
}
model {
target += b * sum(a * Sigma);
}

I'm inclined to go with (1) because it's more efficient, even though
(2) may be a bit more numerically stable.

Also, the end of the function looks like so:

Matrix<T, Dynamic, Dynamic> result
= multiply_lower_tri_self_transpose(L);
symmetrize(result);
return result;

QUESTION 2: Should I push the symmetrization down to multiply_lower_tri_self_transpose?

QUESTION 3: Assuming we go with option (1), does anyone have a better
suggestion than:

for (index_t i = 1; i < a.rows(); ++i)
for (index_t j = 0; j < i; ++j)
a(i, j) = a(j, i);

I'm thinking I could transpose up front then do some fancy head operations
to do the assignment, but I'm not sure it's worth the complexity here.

- Bob

P.S. Don't worry, I'll turn the model test into a unit test with:

test/unit/math/mix/mat/functor/finite_diff_test.cpp

The testing for the transforms is minimal at the moment. The rev
test only tests the value, not the gradients, but I'll add the
gradient tests.

[Ben Goodrich]

On Wednesday, August 10, 2016 at 5:56:24 PM UTC-4, Bob Carpenter wrote:

QUESTION I:  Does anyone have a preference as to how I do
the symmetrization?  The two obvious options are

  (1)  copying below-diagonal elements to above-diagonal

  (2)  adding matrix to its own transpose and dividing by 2

multiply_lower_tri_self_transpose() already does (1); the bug is that cov_matrix_constrain.hpp calls multiply_lower_tri_self_transpose() only in the comments in the version that increments lp__.

Ben

[Bob Carpenter]

Thanks. I missed the double assignment in multiply_lower_tri_self_transpose:

LLt(n, m) = LLt(m, n) = Lt.col(m).head(k).dot(Lt.col(n).head(k));

That makes it trivial to fix the transform asymmetry, which I think you
were trying to tell me from the get go.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.