search

stan-dev / stan

Stan development repository. The master branch contains the current release. The develop branch contains the latest stable development. See the Developer Process Wiki for details.
https://mc-stan.org
BSD 3-Clause "New" or "Revised" License
2.59k stars 369 forks source link

Cholesky factor type for correlation matrices #413

Closed bob-carpenter closed 10 years ago

bob-carpenter commented 10 years ago

I've already started a branch:

feature/corr-cholesky-type

and what remains is:

Unit testing's going to be a big hurdle here because nothing that came before has proper unit tests.

I have a question about whether we want to do the following, and if so if there's a requirement that M < N or something like that (for this and covariance matrix factors, which we'll rename):

bob-carpenter commented 10 years ago

Moving this discussion from stan-dev; questions are mine, answers are Ben's:

Q: I pushed a branch from develop that has everything stubbed out that we'll need to add a Cholesky factor of correlation matrices type:

feature/corr-cholesky-type

Sorry there's not an issue, but I started doing this on the plane w/o an internet connection. Is there a way to change the name of a branch?

A: 

git branch -m

Q: What we need to add is:

A: read_corr_L already does the transformation from CPCs to L I just pushed factor_L which does the transformation from L to CPCs We already have the infrastructure in place to transform CPCs to or from the unconstrained space

Q: Do I need to do anything more in transform.hpp?

A: We also need a density for L such that LL' is LKJ(eta).

Q: Does this also need to allow rectangular sizes?

A: I suppose, but we can do the square ones first.

Q: And are there constraints on the rectangular sizes for the ordinary Cholesky factors for symmetric, positive definite matrices?

A: At least as many rows as there are columns

bob-carpenter commented 10 years ago
bob-carpenter commented 10 years ago

Why is there a cov_matrix_free_lkj() function in prob/transforms.hpp? It's not getting used anywhere I can see in src/stan.

bob-carpenter commented 10 years ago

Should we have cholesky_corr_free (perhaps via factor_L) throws an exception if L is not a Cholesky factor for a correlation matrix?

This is another example of a function that's being called only from our own code when we know it will work because it's going to operate on a variable of type cholesky_corr.

@bgoodri Why go through an array in factor_L the body of which is L.diagonal().array().inverse().diagonal() ?

bgoodri commented 10 years ago

For a large matrix, verifying that the rows of L are normalized would be somewhat expensive for an internal function. I think L.diagonal() is an expression template for a vector, which doesn't have an .inverse() method. Maybe that is wrong, but in any event, .array() and .matrix() are no-ops.

bgoodri commented 10 years ago

I think cov_matrix_free_lkj() is left over from the days when the mapping from a covariance matrix to the unbounded space first decomposed the covariance matrix into a diagonal matrix of standard deviations and a correlation matrix, which were separately mapped to the unbounded space. Now that the primitive of a covariance matrix is just a Cholesky factor, I don't think cov_matrix_free_lkj() is needed

bob-carpenter commented 10 years ago

It's not technically a no-op --- it's an expression template and hence requires a constructor call, which requires a bit of stack-based storage (not heap, so no malloc penalty). The indirection in calls through the expression template should be inlineable away, though.

I'm not sure about whether the expression template for .diagonal() supports .inverse() --- it's the Wild West with template specializations in C++!

On 12/7/13, 6:14 PM, bgoodri wrote:

For a large matrix, verifying that the rows of L are normalized would be somewhat expensive for an internal function. I think L.diagonal() is an expression template for a vector, which doesn't have an .inverse() method. Maybe that is wrong, but in any event, .array() and .matrix() are no-ops.

bgoodri commented 10 years ago

It inherits from MatrixBase so .inverse() for something not square yields an error novel

#include<Eigen/Dense>
#include<iostream>

int main() {
  Eigen::Matrix3d m = Eigen::Matrix3d::Random();
  std::cout << m.diagonal().inverse() << std::endl;
  return 0;
}

Ben

On Sat, Dec 7, 2013 at 8:53 PM, Bob Carpenter notifications@github.comwrote:

It's not technically a no-op --- it's an expression template and hence requires a constructor call, which requires a bit of stack-based storage (not heap, so no malloc penalty). The indirection in calls through the expression template should be inlineable away, though.

I'm not sure about whether the expression template for .diagonal() supports .inverse() --- it's the Wild West with template specializations in C++!

  • Bob

On 12/7/13, 6:14 PM, bgoodri wrote:

For a large matrix, verifying that the rows of L are normalized would be somewhat expensive for an internal function. I think L.diagonal() is an expression template for a vector, which doesn't have an .inverse() method. Maybe that is wrong, but in any event, .array() and .matrix() are no-ops.

— Reply to this email directly or view it on GitHubhttps://github.com/stan-dev/stan/issues/413#issuecomment-30072859 .

bob-carpenter commented 10 years ago

@bgoodri Could you add a test for factor_L in transform.hpp? I can't get it instantiated:

src/stan/prob/transform.hpp:118:54: error: no member named 'diagonal' in
      'Eigen::ArrayBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const
      Eigen::ArrayWrapper<const Eigen::Diagonal<const Eigen::Matrix<double, -1, -1, 0, -1, -1>, 0> > >
      >'
        = L * L.diagonal().array().inverse().array().diagonal();
              ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ^
src/stan/prob/transform.hpp:1499:7: note: in instantiation of function template specialization
      'stan::prob::factor_L<double>' requested here
      factor_L(cpcs,L);
      ^
src/stan/io/writer.hpp:407:13: note: in instantiation of function template specialization
      'stan::prob::cholesky_corr_free<double>' requested here
          = stan::prob::cholesky_corr_free(y);

Also, what's the right drill w.r.t. Git here? I can probably stub things out to where they don't work, but at least they'll compile. Should I do that and push it?

bob-carpenter commented 10 years ago

We can worry about efficiency later -- I can't even get it to compile.

I have everything else compiling now, but am blocked on how to get the factor_L thing working. (In my branch, I've just stubbed out the bit that won't compile.)

Here's the minimal example of the expression you're trying to evaluate in factor_L() in prob/transform.hpp:

#include <Eigen/Dense>

int main() {
   Eigen::MatrixXd L(2,2);
   L << 1, 0, 0, 1;
   Eigen::MatrixXd foo = L * ( L.diagonal().array().inverse().diagonal() );
}

It produces this error:

$ clang++ -I ~/stan/lib/eigen_3.2.0/ test.cpp

test.cpp:6:62: error: no member named 'diagonal' in 'Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const
       Eigen::ArrayWrapper<Eigen::Diagonal<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 0> > >'
   Eigen::MatrixXd foo = L * ( L.diagonal().array().inverse().diagonal() );
                               ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ^
1 error generated.
bgoodri commented 10 years ago

I fixed it. It needed to be asDiagonal() rather than diagonal() and with a vector that inherits from MatrixBase rather than ArrayBase. The test I added might conflict with your unpushed stuff but it should be easy to harmonize it.

bob-carpenter commented 10 years ago

Thanks. I'll merge and give it a go.

On 12/8/13, 9:43 AM, bgoodri wrote:

I fixed it. It needed to be asDiagonal() rather than diagonal() and with a vector that inherits from MatrixBase rather than ArrayBase. The test I added might conflict with your unpushed stuff but it should be easy to harmonize it.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/413#issuecomment-30082895.