symmetry in Wishart RNG

[syclik]

From @bob-carpenter on December 16, 2015 19:28

Andrew reported that the Wishart RNG was returning non-symmetric returns, so we need to either add it to its transpose and divide by 2. You can see in the code that there's no guarantee that everything's going to be symmetric:

     MatrixXd B = MatrixXd::Zero(k, k);
     for (int j = 0; j < k; ++j) {
        for (int i = 0; i < j; ++i)
          B(i, j) = normal_rng(0, 1, rng);
        B(j, j) = std::sqrt(chi_square_rng(nu - j, rng));
      }

      return stan::math::crossprod(B * S.llt().matrixU());

• [x] fix wishart_rng
• [x] fix inv_wishart_rng
• [x] check lkj_corr_rng and see if that also needs to be symmetrized

Copied from original issue: stan-dev/stan#1746

[syclik]

From @andrewgelman on December 16, 2015 23:43

Hi, I was also getting the occasional simulated matrix that was not positive definite. I assume this is roundoff error as well. I don’t know if this is important enough to worry about, I just wanted to let you know it’s happening on rare occasions. A

On Dec 16, 2015, at 2:38 PM, Bob Carpenter notifications@github.com wrote:

Andrew reported that the Wishart RNG was returning non-symmetric returns, so we need to either add it to its transpose and divide by 2. You can see in the code that there's no guarantee that everything's going to be symmetric:

 MatrixXd B = MatrixXd::Zero(k, k);
 for (int j = 0; j < k; ++j) {
    for (int i = 0; i < j; ++i)
      B(i, j) = normal_rng(0, 1, rng);
    B(j, j) = std::sqrt(chi_square_rng(nu - j, rng));
  }

  return stan::math::crossprod(B * S.llt().matrixU());

fix wishart_rng fix inv_wishart_rng check lkj_corr_rng and see if that also needs to be symmetrized — Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/1746.

[syclik]

From @bgoodri on February 9, 2016 1:55

I believe lkj_corr_rng uses multiply_lower_tri_self_transpose or whatever it is called, so it should be symmetric. Perhaps wishart_rng should too, although it would require its code to be transposed relative to what it is now.

[bob-carpenter]

The wishart_rng seems to be OK---it uses the Bartlett decomposition to generate a Choleksy factor for a unit covariance Wishart matrix then Cholesky decomposes the parameter, multiples, then takes a cross-product (see Wikipedia page on . But the inverse Wishart then just generates a Wishart then tries to invert that, which is where the non-positive-definiteness problem arises. I think we should just do the inversion of the Cholesky factors, but I don't know how to do that in Eigen efficiently.