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.)
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
[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/
