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: [email protected] (Vince Wiggins, StataCorp)
To: [email protected]
Subject: Re: st: Matrix inversion "bug"
Date sent: Tue, 10 Jun 2003 15:11:27 -0500
Send reply to: [email protected]
> Mark Schaffer <[email protected]> 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 <[email protected]> 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
> [email protected]
>
> *
> * 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/
