[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

From |
vwiggins@stata.com (Vince Wiggins, StataCorp) |

To |
statalist@hsphsun2.harvard.edu |

Subject |
Re: st: Matrix inversion "bug" |

Date |
Tue, 10 Jun 2003 15:11:27 -0500 |

Mark Schaffer <M.E.Schaffer@hw.ac.uk> notes that occasionally Stata will produce a non-symmetric matrix after computations that are known mathematically to result in a symmetric matrix, and that the resulting matrix cannot be inverted by -syminv(), > Every now and then I encounter a matrix inversion "bug" in Stata. > It cropped up yet again last week when an ivreg2 user wrote to me, > and I finally decided to see what my fellow Statalisters think of > it. > > Say we have a matrix M which is calculated as X*S*X'. S is > symmetric and so is M. We want to invert M and so the best thing to > do is to use -syminv-. The Stata programming manual tells us to use > -syminv- instead of -inv- wherever possible. > > The problem is this. Although M is guaranteed to be symmetric in > principle, it is not guaranteed to be symmetric in practice. Tiny > rounding errors can arise when Stata multiplies matrices. Once in a > while, M will be very slightly non-symmetric, but when it is, > -syminv- will exit with an error. Mark also provided a work-around for the problem, > When I first encountered this, I wrote to Stata Technical Support, > and I received a very helpful and simple fix: prior to the call to > syminv, simply do the following: > > mat M = (M+M')/2 > > This makes the matrix symmetric. Works, no problem. But... He then suggests that this problem be covered up with an official alternate inverter that pre-symmetrize the matrix, > Arguably, this is just a fix for what is really a Stata bug. Ought > it be possible for Stata's code to recognise when the result of a > matrix multiplication is supposed to be a symmetric matrix? > [...] > (3) Neither of the above is very pretty. An alternative is to ask > our Stata Corp friends to, say, add a variant of the -syminv- > function, say -syminv2-. This function would automatically > symmetrize a square matrix before inverting it. And Al Feiveson <alan.h.feiveson@nasa.gov> voiced agreement. > Mark - I too have been bitten by this one - and I agree - it is a > nuisance to have to "symmetrize" the matrix by hand before calling > -syminv- or to use -capture-. > > I like your idea for a syminv2 function that always returns [(A + > A')/2 ]-1 no matter what A is as long as A is square and A+A' is > nonsingular. While I almost always nod in agreement when reading posts by Mark or Al, this suggestion actually gave me a cold shiver. (OK, admittedly I should get out more and experience some truly scary things.) When programming statistics on a computer, we can ignore the precision of the computations about 99.99% of the time. We can just pretend we are implementing mathematical relations precisely. Occasionally, however, the computations go a little, or maybe a lot, awry and we are forced to remember that we are really programmers, working with finite-precision representations of our mathematical relations. We probably should never really forget this, but it is a convenient fiction and I spend as much time as I can in that fiction. When the computer stops us and yells "your assumptions don't hold in my finite precision world!", we need to stop and ask ourselves if our precision problem is small enough to be ignored or if it is likely to lead to such imprecision that our inferences may be affected. It is possible that (1) we have grouped our computations in a way that compounds precision problems, or (2) that we have a computation that is inherently imprecise, or (3) that our computation is imprecise for many datasets. If it is any of these three, then we have real work to do in determining how to repair our computation or when to allow our computation. Automatically symmetrizing a matrix in such cases would hide a real problem. Adding an extra statement in to symmetrize the matrix, is a message that we have thought about the computational issue and deemed that there is no real problem with our computation. As a side note, Stata has a built-in tolerance for matrix symmetry. For a matrix to be considered symmetric each off-diagonal symmetric element pair -- d = X[i,j] - X[j,i] must meet the condition d^2 < 1e-16*abs(X[i,i])*abs(X[j,j]) with a different, somewhat tighter criterion when x[i,ik] or x[j,j] is exactly zero. We occasionally use the trick that Mark mentions to symmetrize a matrix before inversion, but we never do it without first considering other options and evaluating whether we believe it is computationally justified in the specific context. The function -syminv()- is used 127 times in 55 ado-files in official Stata. In all those uses, the matrix is pre-symmetrized in only 9 cases, and three of those are repeats for older version of code retained for version control. So official Stata has only 6 computations where pre-symmetrizing the matrix was necessary and deemed safe. If Mark's computation is what I think it is, I would be inclined to symmetrize the matrix. I don't think there is any inherent problem with the computation. I just don't want to make it too easy to make that decision. -- Vince vwiggins@stata.com * * For searches and help try: * http://www.stata.com/support/faqs/res/findit.html * http://www.stata.com/support/statalist/faq * http://www.ats.ucla.edu/stat/stata/

**Follow-Ups**:**Re: st: Matrix inversion "bug"***From:*"Mark Schaffer" <M.E.Schaffer@hw.ac.uk>

- Prev by Date:
**RE: st: New package -ingap- on SSC** - Next by Date:
**st: is ordering with -bysort- unique?** - Previous by thread:
**st: RE: RE: Matrix inversion "bug"** - Next by thread:
**Re: st: Matrix inversion "bug"** - Index(es):

© Copyright 1996–2015 StataCorp LP | Terms of use | Privacy | Contact us | What's new | Site index |