loosen constraints on simplexes, covariance, etc.

[bob-carpenter]

Requested by Dan Stowell (most recently, and ourselves in the past):

We're gettig too many errors from floating-point rounding, so we should probably relax our constraints a bit. I'm not sure exactly how to pull together a list of functions that check constraints other than reading the manual closely and pulling them out.

• math functions
• distribution checks
• error function prints at full precision (to avoid embarssing error messages like "should be 1, but is 1")

[bgoodri]

In transform.hpp, CONSTRAINT_TOLERANCE is defined as 1E-8 . How would well-written code with double-precision variables end up at single-precision?

[bob-carpenter]

On 11/5/13, 8:23 PM, bgoodri wrote:

In transform.hpp, CONSTRAINT_TOLERANCE is defined as 1E-8 .
How would well-written code with double-precision variables
end up at single-precision?

I can never tell what's a rhetorical question in e-mail, so
I'll answer this one.

For a start, I think our CONSTRAINT_TOLERANCE is probably (improperly)
being used to measure absolute error rather than relative error, and
of course floating-point precision is relative, not absolute.

I believe the main place we're having problems is with correlation
and covariance matrices violating positive-definiteness
constraints. See the related (if not identical) issue

#240

The function doing the positive-definiteness testing is here:

src/stan/math/error_handling/matrix/check_pos_definite.hpp

So you (Ben) probably have a better handle on the issue than the rest of us!

Also, what should we print out in case of error? I don't see the point of
printing out just entry (0,0) at default precision.

And if somebody fixes this function, can they change this mess:

   if(cholesky.info() != Eigen::Success ||
      cholesky.isNegative() ||
      (cholesky.vectorD().array() <= CONSTRAINT_TOLERANCE).any())

to something with the continuing operators at the start of a line to
follow both standard math and programming practice (R requires the ops
at the end because of its line-based parser):

      if (cholesky.info() != Eigen::Success
          || cholesky.isNegative()
          || (cholesky.vectorD().array() <= CONSTRAINT_TOLERANCE).any())

vectorD is the diagonals of LDLT (I have no idea what that is or why it should be
less than the tolerance), and any() just checks if any of the elements of
a matrix are true. In this case, it would almost certainly be faster not
to use this kind of vectorization since we don't need the intermediate vector
that gets constructed. But it is convenient to have it in a single line. I
don't know if there's any utility to indicating which of these conditions fails.

Can't we just replace the last two conditions with (!cholesky.isPositive())
since isPositive() is the test for positive (semi) definiteness?

• Bob

[betanalpha]

For a start, I think our CONSTRAINT_TOLERANCE is probably (improperly)
being used to measure absolute error rather than relative error, and
of course floating-point precision is relative, not absolute.

Between this and the possible subtraction of two large numbers followed by a multiply,
it's pretty easy to lose double precision.

[bob-carpenter]

May need to write more of the LKJ transforms.

[bob-carpenter]

Copied from #240:

Jo(h)n Smith reported the following issue on the mailing list.

// File: corr-init.stan

parameters {
  corr_matrix[100] Omega;
}
transformed parameters {
  corr_matrix[100] OmegaCopy;
  OmegaCopy <- Omega;
}
model {
}

This is still an issue in Stan 2, where running this:

./corr-init sample

produces a stream of errors of the form:

Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:
Invalid value of OmegaCopy: Error in function validate transformed params N4stan5agrad3varE: yy is not positive definite. y(0,0) is 1:0.
If this warning occurs sporadically then the sampler is fine,
but if this warning occurs often then your model is either severely ill-conditioned or misspecified.
Rejecting proposed initial value with zero density.

and then finally dies with

Initialization between (-2, 2) failed after 100 attempts. 
 Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

Maybe we can (a) initialize randomly a little better by shrinking to identity matrix, or (b) turn down the stringency of our positive-definiteness tests.

[bob-carpenter]

Does anyone have any idea how to fix any of this? I don't, but I'd really like someone to look at it. Ben? I'm pushing it off to 2.4.0+ so it doesn't hold us up for 2.4.

[bob-carpenter]

See #817 for an issue with the Wishart that we should be able to fix.

[syclik]

@bgoodri or @mbrubake, do you have any suggestions?

[syclik]

@bgoodri, @mbrubake: I am closing this issue. If you have thoughts, please reopen.