[ADMB Users] ADMB-RE fitting failure in 11 but not in 10.1

[Mark Wolters]

Hi All,

My two cents here--given with limited knowledge of the problem at
hand, and of how ADMB handles matrix computations internally:

I infer from the previous emails that singularity of a matrix is being
checked by looking at the determinant, which is calculated as the
product of eigenvalues.  Naturally when the determinant is zero to
machine precision, this should cause an error.  But even if you change
the way you compute the determinant to make it "not quite" machine
zero, you probably still have an ill-conditioned matrix, and the
system should probably throw a warning. Just echoing what Ben said
here.

My understanding is that the usual "condition number" used to check
for problems with invertibility is the ratio of the largest to the
smallest eigenvalue.  If this ratio is close to 1, then things are
good.  As it gets larger, accuracy decreases.  Not being a numerical
analyst, I don't know what would be a typical cutoff for the condition
number to consider a matrix seriously ill-conditioned and worthy of a
warning. Maybe one could inspect how this is done in, e.g., R or
MATLAB?

Final comment, if the problem at hands yields matrices so
ill-conditioned that the manner of calculating the determinant becomes
significant, maybe it's better to re-parametrize the model, do a
transformation, start from different starting values, etc.?

This problem is a bit tricky since, if the singularity problem arises
early in a search, you probably just want to survive it and keep going
until the search moves to a better region of the parameter space.

Mark

On Thu, Dec 6, 2012 at 4:38 PM, Ben Bolker <bbolker at gmail.com> wrote:
> On 12-12-06 11:47 AM, Tim Miller wrote:
>> Dear Users,
>>
>> During the fitting of a random effects model in version 11 (Windows 7,
>> MinGW GCC-4.5.2), I get a failure with error message:
>>
>> Error in matrix inverse -- matrix singular in solve(dmatrix)
>>
>> This occurs before a .par file is written. In 10.1 convergence is fine.
>> The error doesn't occur in 11 when I use the sparse hessian option
>> -shess and in that case I get the same result as 10.1 (with or without
>> -shess).
>>
>> I think I have traced it to the solve function in ludcmp_solve.cpp just
>> by searching for the error text in all the source code files and I think
>> the error occurs when checking the determinant.  In version 11 the
>> function appears to check whether exp(sum(log(fabs(determinant
>> components)))) == 0.0 whereas in 10.1 it checks product(determinant
>> components) == 0.0.
>>
>> Assuming I have traced the source of the error correctly, why is the
>> newer determinant checking method better than that in 10.1? If it is
>> better, maybe there is a way to continue on rather than exiting on the
>> error?
>>
>> I could provide the tpl/dat files, but fitting the model without -shess
>> is annoying because it takes a while to get to the error.
>>
>> Tim
>>
>
>   Hmm.  I'm trying to track these changes in the SVN repository.
>
> I did
>
>    find . -name "*.cpp" -exec grep -H singular {} \;
>
> to find every instance of "singular" in the code ...
>
> found ludcmp_solve.cpp in contrib/ludcmp/linad99 (took me a little while
> since originally I was looking just in src/)
>
>   As Mollie pointed out,
>
> ================
> r303 | mbrooks | 2012-03-15 21:31:26 -0400 (Thu, 15 Mar 2012) | 1 line
>
> sum log det instead of product of det to check invertibility
> ===============
>
>   but a less stringent (possibly better??) check is (as implemented
> elsewhere in the code) to check individual elements of the products to
> make sure none of them are zero.  (The product can underflow even if
> none of the individual elements is numerically equal to zero.)
>
>   This is fundamentally a difficult problem (I think; I would be happy
> to be corrected if someone has an easy answer!).  The rules
>
> a. check whether product is zero, or close to zero
> b. check whether exp(sum(log(...))) is zero, or close to zero
> c. check whether individual elements to be multiplied are zero, or close
> to zero
>
>   are _mathematically_ equivalent, but they aren't computationally
> equivalent. It may be possible to prove that one of them is always more
> sensitive (i.e. more likely to detect singularity), but I wouldn't be
> surprised if they differ depending on the scenario.  Furthermore (and
> this is the important part), **it's not clear what the 'right' answer
> is** in terms of sensitivity.  There is a very sensible argument that
> says that the software should detect when the user is trying to invert a
> matrix that is very ill-conditioned/close to singularity, and at least
> warn and possibly give an error -- the argument being that the answers
> may be extremely unreliable in these situations, and that it isn't
> ethical to allow the user to proceed (thinking they got an answer when
> in fact it's rubbish) in this case.
>
> PS: checking (a) vs (b) in R shows that the answers differ for different
> sets of random numbers -- sometimes (a) is bigger, sometimes (b) is bigger.
>
>> tmpf <- function() { r <- 10^runif(100,min=-3,max=-1);
> prod(r)-exp(sum(log(r))) }
>> replicate(20,tmpf())
>  [1]  7.556839e-210 -2.392981e-220 -5.453231e-218  2.266101e-208
> -1.133139e-216
>  [6] -1.728037e-217 -3.965986e-218 -1.867355e-221 -3.146840e-221
> 9.885676e-210
> [11] -9.282674e-214 -1.682003e-218  1.822128e-205 -1.469539e-218
> 5.105471e-214
> [16]  1.827526e-214 -7.008184e-207  4.061170e-215 -1.594751e-202
> 4.641337e-213
>
>   Which one is bigger is pretty much random:
>
>> table(sign(replicate(1000,tmpf())))
>
>  -1   0   1
> 507   6 487
>
>
>   My only argument would be that it might be a good idea if the 23
> places that singularity is checked in the ADMB code (find . -name
> "*.cpp" -exec grep -H "Error in matrix inverse.*singular" {} \; | wc )
> behaved *consistently*, to the extent that was possible.  It might be
> nice, *if* it seemed to make sense, to allow some degree of user control
> (either provide parallel sloppy versions of the functions, or a global
> argument -danger_danger_ignore_singularity_checks ...) to allow the
> computation to go forward under these circumstances.
>
>   However, since I'm not volunteering to do the work myself, this is
> just a suggestion ...
>
>   cheers
>     Ben Bolker
>
> _______________________________________________
> Users mailing list
> Users at admb-project.org
> http://lists.admb-project.org/mailman/listinfo/users