Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: Matrix inversion "bug"


From   "Mark Schaffer" <[email protected]>
To   [email protected]
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:           	[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/



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