Pushed check_pos_definite specialization

[Marcus Brubaker]

Just for kicks this afternoon I implemented a check_pos_definite for matrices of vars which does the positive definiteness test on a double matrix instead to avoid the wasteful creation of vars.  This should provide a moderate speedup to anyone who has cov_matrix or corr_matrix parameters.  As it turns out, it won't impact the multi_normal case because the positive definiteness test is done differently there.

Cheers,
Marcus

[Daniel Lee]

Thanks!

If there isn't a good reason for multi_normal to be different, would you mind changing the implementation so it matches?

--
 
 

[Ben Goodrich]

On Saturday, October 13, 2012 6:56:41 PM UTC-4, Daniel Lee wrote:

Thanks!

If there isn't a good reason for multi_normal to be different, would you mind changing the implementation so it matches?

Our check for positive definiteness in the multi_normal case is whether the Cholesky factorization works, in which case it calls multi_normal_cholesky(). So, if Sigma is a parameter, we want that Cholesky factorization to return a matrix of vars and not check positive definiteness twice. All of this will hopefully be refactored in Stan 2.x once I implement the CovarianceMatrix class(es).

Ben

[Ben Goodrich]

On Saturday, October 13, 2012 6:47:42 PM UTC-4, Marcus Brubaker wrote:

Just for kicks this afternoon I implemented a check_pos_definite for matrices of vars which does the positive definiteness test on a double matrix instead to avoid the wasteful creation of vars.  This should provide a moderate speedup to anyone who has cov_matrix or corr_matrix parameters.  As it turns out, it won't impact the multi_normal case because the positive definiteness test is done differently there.

 
So, I guess this specialization won't get called that often, but just to make sure I understand (and to bring together the two threads), how is what you just pushed different from the compiler's perspective than just changing src/stan/math/matrix_error_handling.hpp to

    template <typename T_y, typename T_result, class Policy>
    inline bool check_pos_definite(const char* function,
                  const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
                  const char* name,
                  T_result* result,
                  const Policy&) {
      typedef
        typename Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>::size_type
        size_type;
      if (y.rows() == 1 && y(0,0) <= CONSTRAINT_TOLERANCE) {
        std::ostringstream message;
        message << name << " is not positive definite. "
                << name << "(0,0) is %1%.";
        T_result tmp = policies::raise_domain_error<T_y>(function,
                                                         message.str().c_str(),
                                                         y(0,0), Policy());
        if (result != 0)
          *result = tmp;
        return false;
      }
      if ((y.cast<double>().ldlt().vectorD().array() <= CONSTRAINT_TOLERANCE).any()) {     
        std::ostringstream message;
        message << name << " is not positive definite. "
                << name << "(0,0) is %1%.";
        T_result tmp = policies::raise_domain_error<T_y>(function,
                                                         message.str().c_str(),
                                                         y(0,0), Policy());
        if (result != 0)
          *result = tmp;
        return false;
      }
      return true;
    }

Thanks,
Ben

[Marcus Brubaker]

That would probably work, but is it possible to cast from var to double?  I thought you couldn't convert to double without calling .val() which would break the cast call.  The other thing (although less relevant) is that what you posted would also convert floats, ints, long doubles, etc to double.

Cheers,
Marcus

[Bob Carpenter]

On 10/13/12 7:59 PM, Marcus Brubaker wrote:

> On Sat, Oct 13, 2012 at 7:18 PM, Ben Goodrich <goodri...@gmail.com <mailto:goodri...@gmail.com>> wrote:
>
> On Saturday, October 13, 2012 6:47:42 PM UTC-4, Marcus Brubaker wrote:
>
> Just for kicks this afternoon I implemented a check_pos_definite for matrices of vars which does the positive
> definiteness test on a double matrix instead to avoid the wasteful creation of vars. This should provide a
> moderate speedup to anyone who has cov_matrix or corr_matrix parameters. As it turns out, it won't impact the
> multi_normal case because the positive definiteness test is done differently there.

We should go back and check where the tests are getting done
with the transforms. It should help there if we're testing.
That would wind up testing redundantly in the case of a cov_matrix
that was unmodified and used in the normal.

Also, anything that's just constant double shouldn't need to
be continually checked.

These are the things that Ben's new design for positive def
matrices with caching should address. (They're orthogonal
to the issue of speeding up from pulling the value out of
var and calling the tests with double -- that should be a
pretty big win in terms of auto-dif vars generated, so that
the copy will be of negligible cost in comparison.)

No -- it won't work for just that reason.

We don't want to add an implicit cast from var to double, either.

You probably want two functions to make the templates easier
to go down.

I'd just pull the test code into a function validate_pos_def()
and use:

void validate_pos_def(const Matrix<var,Dynamic,Dynamic>& m)
// convert to double matrix
// call shared test code

void validate_pos_def(const Matrix<double,Dynamic,Dynamic>& m) {
// call shared test code

Trying to share more than that in a single function is going to
lead to trouble because both branches of any if statement must
compile.

> I thought you couldn't convert to double
> without calling .val() which would break the cast call.

Correct.

> The other thing (although less relevant) is that what you
> posted would also convert floats, ints, long doubles, etc to double.

It'd be best to avoid this, but I'm not too worried about
making our libs play nicely with float or long double types.
Many things that should probably be more general are already
somewhat or fully hard-coded for double.

- Bob

[Ben Goodrich]

On Saturday, October 13, 2012 8:29:05 PM UTC-4, Bob Carpenter wrote:

On 10/13/12 7:59 PM, Marcus Brubaker wrote:
> On Sat, Oct 13, 2012 at 7:18 PM, Ben Goodrich <goodri...@gmail.com <mailto:goodri...@gmail.com>> wrote:
>
>     On Saturday, October 13, 2012 6:47:42 PM UTC-4, Marcus Brubaker wrote:
>
>         Just for kicks this afternoon I implemented a check_pos_definite for matrices of vars which does the positive
>         definiteness test on a double matrix instead to avoid the wasteful creation of vars.  This should provide a
>         moderate speedup to anyone who has cov_matrix or corr_matrix parameters.  As it turns out, it won't impact the
>         multi_normal case because the positive definiteness test is done differently there.

We should go back and check where the tests are getting done
with the transforms.  It should help there if we're testing.
That would wind up testing redundantly in the case of a cov_matrix
that was unmodified and used in the normal.

Also, anything that's just constant double shouldn't need to
be continually checked.

Actually, check_positive_definite is currently not being called at all:

goodrich@CYBERPOWERPC:/opt/stan$ git grep -n "check_pos_definite" src/stan/ | cat

src/stan/agrad/matrix.hpp:1980:    inline bool check_pos_definite(const char* function,

src/stan/agrad/matrix.hpp:1992:        return stan::math::check_pos_definite(function,y_d,name,result,Policy());

src/stan/math/matrix_error_handling.hpp:136:    inline bool check_pos_definite(const char* function,

src/stan/math/matrix_error_handling.hpp:173:    inline bool check_pos_definite(const char* function,

src/stan/math/matrix_error_handling.hpp:177:      return check_pos_definite(function,y,name,result,default_policy());

src/stan/math/matrix_error_handling.hpp:182:    inline bool check_pos_definite(const char* function,

src/stan/math/matrix_error_handling.hpp:186:      return check_pos_definite(function,y,name,result,default_policy());

src/stan/math/matrix_error_handling.hpp:217:      if (!check_pos_definite(function, y, name, result, Policy()))

src/stan/math/matrix_error_handling.hpp:285:      if (!check_pos_definite(function, y, "y", result, Policy()))

We do check for positive definiteness in multi_normal and multi_student_t, but as part of the density evaluation rather than a separate call to check_pos_definite(). It is a shame that fixed matrices get checked every time the density is evaluated, but I am not sure there is much we can do about it with the Stan 1.0.x design.

 

> That would probably work, but is it possible to cast from var to double?

No -- it won't work for just that reason.

We don't want to add an implicit cast from var to double, either.

Ah, I had forgotten that casting doesn't go that direction.

Ben

[Marcus Brubaker]

It gets called from check_cov_matrix and check_corr_matrix which get called when there are cov_matrix and corr_matrix parameters and one or two other places.
 

Cheers,
Marcus

[Ben Goodrich]

On Sunday, October 14, 2012 12:00:01 PM UTC-4, Marcus Brubaker wrote:

It gets called from check_cov_matrix and check_corr_matrix which get called when there are cov_matrix and corr_matrix parameters and one or two other places.

Aha!.

goodrich@CYBERPOWERPC:/opt/stan$ git grep -n "check_corr_matrix" src/stan | cat
src/stan/io/reader.hpp:905:       * <p>See <code>stan::math::check_corr_matrix(Matrix)</code>.
src/stan/io/reader.hpp:913:        stan::math::check_corr_matrix("stan::math::corr_matrix(%1%)", x, "Constrained matrix");
src/stan/io/writer.hpp:334:        stan::math::check_corr_matrix("stan::io::corr_matrix_unconstrain(%1%)", y, "Matrix");
src/stan/math/matrix_error_handling.hpp:258:    inline bool check_corr_matrix(const char* function,
src/stan/math/matrix_error_handling.hpp:292:    inline bool check_corr_matrix(const char* function,
src/stan/math/matrix_error_handling.hpp:296:      return check_corr_matrix(function,y,name,result,default_policy());
src/stan/math/matrix_error_handling.hpp:300:    inline bool check_corr_matrix(const char* function,
src/stan/math/matrix_error_handling.hpp:304:      return check_corr_matrix(function,y,name,result,default_policy());
src/stan/prob/distributions/multivariate/continuous/lkj_corr.hpp:131:      using stan::math::check_corr_matrix;
src/stan/prob/distributions/multivariate/continuous/lkj_corr.hpp:141:      if (!check_corr_matrix(function, y, "Correlation matrix", &lp, Policy())) {

We do need to check in the lkj_ case, but if the shape is 1.0 and you know it is a valid correlation matrix, then it is a waste of time to call lkj_corr(1.0).

goodrich@CYBERPOWERPC:/opt/stan$ git grep -n "check_cov_matrix" src/stan | cat
src/stan/io/reader.hpp:950:       * <p>See <code>stan::math::check_cov_matrix(Matrix)</code>.
src/stan/io/reader.hpp:959:        stan::math::check_cov_matrix("stan::io::cov_matrix(%1%)", y, "Constrained matrix");
src/stan/math/matrix_error_handling.hpp:206:    inline bool check_cov_matrix(const char* function,
src/stan/math/matrix_error_handling.hpp:223:    inline bool check_cov_matrix(const char* function,
src/stan/math/matrix_error_handling.hpp:227:      return check_cov_matrix(function,y,name,result,default_policy());
src/stan/math/matrix_error_handling.hpp:232:    inline bool check_cov_matrix(const char* function,
src/stan/math/matrix_error_handling.hpp:236:      return check_cov_matrix(function,y,name,result,default_policy());
src/stan/math/matrix_error_handling.hpp:320:    inline bool check_cov_matrix(const char* function,
src/stan/math/matrix_error_handling.hpp:334:    inline bool check_cov_matrix(const char* function,
src/stan/math/matrix_error_handling.hpp:337:      return check_cov_matrix(function,Sigma,result,default_policy());
src/stan/math/matrix_error_handling.hpp:341:    inline bool check_cov_matrix(const char* function,
src/stan/math/matrix_error_handling.hpp:344:      return check_cov_matrix(function,Sigma,result,default_policy());
src/stan/prob/distributions/multivariate/continuous/multi_normal.hpp:47:      using stan::math::check_cov_matrix;
src/stan/prob/distributions/multivariate/continuous/multi_normal.hpp:135:      using stan::math::check_cov_matrix;
 
Looks like we can git rid of those last two using statements.

Ben