Is it possible to disable symmetry checks?


#1

Is it possible to turn off the symmetry check?
Commonly, BLAS calls only look at the lower or upper triangle of the input matrix, treating the other triangle as the mirror image.

I ask because of errors like:

Chain 11: Rejecting initial value:
Chain 11:   Error evaluating the log probability at the initial value.
Chain 11: Exception: Exception: cholesky_decompose: m is not symmetric. m[1,2] = -4.37401, but m[2,1] = -4.37401  (in 'modela8982e43d111_model6_marginalize_interp' at line 178)

For numerical reasons, the matrix will often not be precisely symmetric, although within rounding, as the exception itself makes clear: obviously -4.37401 == -4.37401.

It’s also worth pointing out that I am calculating the symmetric matrix via:

non_missing_block - quad_form(inverse_spd(missing_block), association_block)

instead of

non_missing_block - quad_form_sym(inverse_spd(missing_block), association_block)

because quad_form_sym would complain of similar numerical differences, eg

Exception: Exception: quad_form_sym: A is not symmetric. A[72,78] = 1.87279e+07, but A[78,72] = 1.87279e+07

Floating point errors mean that matrices that ought to be precisely symmetric will often fail the symmetry checks.
For example:

SAMPLING FOR MODEL 'model6_marginalize_interp' NOW (CHAIN 13).
Chain 13: Rejecting initial value:
Chain 13:   Error evaluating the log probability at the initial value.
Chain 13: Exception: Exception: quad_form_sym: A is not symmetric. A[37,68] = 1.03544e+07, but A[68,37] = 1.03544e+07  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 178)
  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 335)

Chain 13: Rejecting initial value:
Chain 13:   Error evaluating the log probability at the initial value.
Chain 13: Exception: Exception: quad_form_sym: A is not symmetric. A[1,67] = -13273.6, but A[67,1] = -13273.6  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 178)
  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 335)

Chain 13: Rejecting initial value:
Chain 13:   Error evaluating the log probability at the initial value.
Chain 13: Exception: Exception: quad_form_sym: A is not symmetric. A[8,68] = 275203, but A[68,8] = 275203  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 178)
  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 335)

Chain 13: Rejecting initial value:
Chain 13:   Error evaluating the log probability at the initial value.
Chain 13: Exception: Exception: quad_form_sym: A is not symmetric. A[1,68] = -97994.3, but A[68,1] = -97994.3  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 178)
  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 335)

Chain 13: Rejecting initial value:
Chain 13:   Error evaluating the log probability at the initial value.
Chain 13: Exception: Exception: quad_form_sym: A is not symmetric. A[1,73] = -54295.9, but A[73,1] = -54295.9  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 178)
  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 335)

Chain 13: Rejecting initial value:
Chain 13:   Error evaluating the log probability at the initial value.
Chain 13: Exception: Exception: quad_form_sym: A is not symmetric. A[56,74] = -4.43735e+07, but A[74,56] = -4.43735e+07  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 178)
  (in 'modelb0b43b4c742c_model6_marginalize_interp' at line 335)

Chain 13 failed on all 100 initial values, causing the chain to abort. Every single one of the error messages looked something like A[56,74] = -4.43735e+07, but A[74,56] = -4.43735e+07 – clearly, the error is smaller than the number of significant figures the print method keeps.

I can of course replace quad_form_sym with quad_form. This seemingly solves the problem.
After another operation, I do a cholesky_decompose, which does sometimes throw an exception, but rarely enough that I’ve not seen it cause a chain to crash.
That is the approach I have generally been forced to take. An alternative is trying to force symmetry in the input.

Is there a way to disable the symmetry check? Or, at least decrease its sensitivity.

EDIT:
For what it is worth, this was the offending bit of code:

quad_form_sym(inverse_spd(missing_block), association_block)

I replaced that with this, which I assume is more stable than just using quad_form, and ought to guarantee symmetry of the output:

crossprod(mdivide_left_tri_low(cholesky_decompose(missing_block), association_block))

which does not seem to throw any errors. I am guessing the inverse_spd produced large magnitudes, which may have tripped an absolute tolerance? Or large magnitudes reduce stability?
Either way – code works, but it would be nice to have more functions with an api like mdivide_left_tri_low, that cast triangular (or symmetric) views by focusing on only a lower/upper half of a matrix.

EDIT:
Runtime was/is a major issue with that model.
So I was shocked to see that switching to the crossprod(mdivide_left_tri_low(cholesky(...))) version decreased the 1000 transition-10-leapfrog time roughly 3-fold (~ 800 seconds to ~250 seconds)!
I am pleased.
Seems like trying multiple variations can be worthwhile. And that I should probably favor “div” methods, even over specialized inverses (like inverse_spd) in general?
EDIT:
Performance difference is closer to 5x based on sampling rate.
Maybe I should be glad I couldn’t disable a symmetry check, otherwise I probably wouldn’t have tried this work-around.


#2

If you want to know why the difference occur it may be that the more numerically stable version is more accurate and results in a higher acceptance rate for a trajectory which in turn allows adaptation to push for a higher stepsize in the integrator. Pretty easy to look for this in diagnostics and it’s generally relevant for speeding up sampling