added function multiply_lower_tri_self_transpose

[Bob Carpenter]

The name's a mouthful, but the behavior is simple:

multiply_lower_tri_self_transpose(L)

= lower_tri(L) * lower_tri(L)'

where lower_tri(L) is the lower triangular portion of
L including the diagonal (thus the strictly upper
triangular portion of input L is ignored even if it
is nonzero).

I added some constructors to the dot_self_vari
and dot_product_vari in order to be able to share
variables in the dot products.

I tested values and gradients vs. the direct
implementation in Eigen, and everything looks good.

I didn't replace anything in the rest of the code
where this might get used, but I did add it to
function_signatures so you can use it in models and
also added it to the doc.

- Bob

[Ben Goodrich]

One place where this could be used is in some of the LKJ stuff, but I am guessing (without having looked) this implementation would calculate the diagonal and the associated gradients even though they are fixed in the LKJ case. Maybe add an option to fix the diagonal to 1.0?

Ben

[Ben Goodrich]

We could also use an optimized implementation of M * M' or M' * M where M is not necessarily triangular. See

http://comments.gmane.org/gmane.comp.lib.eigen/2832

Ben

On Thursday, September 27, 2012 5:19:04 PM UTC-4, Bob Carpenter wrote:

[Bob Carpenter]

On 9/27/12 8:00 PM, Ben Goodrich wrote:
> We could also use an optimized implementation of M * M' or M' * M where M is not necessarily triangular. See
>
> http://comments.gmane.org/gmane.comp.lib.eigen/2832

So that'd be a strictly lower triangular
and full matrix version. It's pretty easy
to add new operations like these -- the
code for L * L' is only 10 or 15 lines in
agrad/matrix.

- Bob

[Ben Goodrich]

On Thursday, September 27, 2012 5:19:04 PM UTC-4, Bob Carpenter wrote:

I tested values and gradients vs. the direct
implementation in Eigen, and everything looks good.

Something is fishy with the gradient tests here: test_mult_LLT_grad is only called by test_mult_LLT_grad, and test_mult_LLT_grad is not called. I noticed this because I was trying to test my implementation of M * M'.

Ben

[Bob Carpenter]

Thanks -- I'll fix it.

- Bob

> --
>
>

[Bob Carpenter]

On 9/27/12 11:57 PM, Ben Goodrich wrote:

I fixed the test. I calculated a 2 x 2 and
3 x 3 case by hand and tested those. I never
did figure out how to get the general test working.

- Bob