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

From |
"Mark Schaffer" <M.E.Schaffer@hw.ac.uk> |

To |
statalist@hsphsun2.harvard.edu |

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

Date |
Wed, 11 Jun 2003 10:08:25 +0100 |

Vince, David, Nick, Al, Thanks for the helpful replies. With respect to Nick's suggestion, > There is an -issym()- function in Stata 8. This makes the hack a bit cleaner - add "if not(issym(M))" before symmetrizing. It doesn't address my main complaint, though, which is the need for the hack. David's suggestion is my favourite: > Instead, I'd suggest that Stata develop an internally created and > maintained property for matrices that indicates whether or not it is > symmetric. Perhaps this is already being done for diagonal matrices. > This property would "know" that the result of particular matrix > expressions should create a symmetric matrix and enforce this by > whatever internal operations are necessary. For example, Stata should > know that if S is symmetric, then X*S*X' is also symmetric. The only insurmountable problem I can see with this is that if the matrix multiplication takes place in two lines, say M=X*S followed by M=M*X', then it's not reasonable to expect Stata to "know" that the result is a symmetric matrix. I admit, it's been a few years since I wrote a parser (closer to a couple of decades is more like it) but in principle it ought to be possible for Stata to look at a matrix multiplication expression and detect whether or not it's a "matrix palindrome". (Take the multiplication, reverse the order of the matrices, transpose each, and end up with the same matrix multiplication you started with.) I am not sure I am entirely convinced by Vince: From: vwiggins@stata.com (Vince Wiggins, StataCorp) To: statalist@hsphsun2.harvard.edu Subject: Re: st: Matrix inversion "bug" Date sent: Tue, 10 Jun 2003 15:11:27 -0500 Send reply to: statalist@hsphsun2.harvard.edu > 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.) I can see why this would cause Vince to shiver. I've often programmed some matrix expression, used syminv to invert it, watched syminv crash, and discovered a bug in my code - the matrix was supposed to be symmetric, and wasn't. But.... > > 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. The problem with Vince's suggestion is that I cannot tell when such a problem is likely to arise. In the case of ivreg2, the bug cropped up, someone sent me a way of replicating it, I found the bug and squashed it in that location in the code, and then a few months later the bug cropped up again but at a different point in the code. My reaction this time was to say "sod it" (local lingo), and I decided to pre-symmetrize every matrix in the code before inverting so the bug would never come up again. Is this a mistake? With respect to Al's point, > Vince - > > But on what basis can you say that A-1 is more "accurate" than > [(A+A')/2 ]-1 (where A is the finite representation of the supposedly > symmetric matrix stored in Stata) when it is already known that A is > supposed to be symmetric? The very fact that you are using -symimv- in > these 55 ado files suggests that you are assuming A IS symmetric. I would like to chime in and ask the following: if we know that the result of some matrix multiplication A is supposed to be symmetric, and it (very slightly) isn't, then does pre-symmetrizing it with (A+A')/2 make it a more accurate in terms of resembling the correct answer to the matrix multiplcation? --Mark > > > -- 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/ Prof. Mark E. Schaffer Director Centre for Economic Reform and Transformation Department of Economics School of Management & Languages Heriot-Watt University, Edinburgh EH14 4AS UK 44-131-451-3494 direct 44-131-451-3008 fax 44-131-451-3485 CERT administrator http://www.som.hw.ac.uk/cert * * 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/

**References**:**Re: st: Matrix inversion "bug"***From:*vwiggins@stata.com (Vince Wiggins, StataCorp)

- Prev by Date:
**Re: st: Fixed effects and time dummy variables** - Next by Date:
**st: Default precision - was Why does Stata do it wrong?** - Previous by thread:
**Re: st: Matrix inversion "bug"** - Next by thread:
**RE: st: Matrix inversion "bug"** - Index(es):

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