Jean Salvati <[email protected]> wants to know more about how syminv() works:
> Specifically, I don't understand how looking for zeroes on the diagonal of
> syminv(X'X) allows you to determine which variables cause singularity
> problems.
The inverse of matrix A is defined as the matrix X solving the equation
A*X = I (1)
where I is the identity matrix. Some matrices are singular, meaning there is
no solution to (1). In that case, it is often sufficient to obtain a
generalized inverse, defined as matrix X meeting the restrictions
A*X*A = A
(2)
X*A*X = X
When A is singular, (2) is not sufficient to identify uniquely an X -- there
are lots of different X's that suffice. For most problems, one X will work
as well as another.
syminv() produces generalized inverses for symmetric, positive definite
matrices, and chooses among the X's the X that minimizes roundoff error given
its approach, its approach being to ignore rows and columns that look
collinear to it. Ignoring rows (columns) amounts to dropping variables. When
syminv() decides to ignore a row (column), it substitutes 0 for the row and
column in the result. Thus, the easy way to determine whether syminv()
ignored a row (column) is to check the diagonal element for zero.
That makes it easy to calculate the rank. If
+- -+
| 1 0 0 |
A * syminv(A) = | 0 1 0 |
| 0 0 1 |
+- -+
then A is full rank. If
+- -+
| 1 0 0 |
A * syminv(A) = | -1 0 1 |
| 0 0 1 |
+- -+
then A has rank 2. More correctly, we should say that syminv() acted as if
A had rank 2. If we tightened down one some of syminv()'s secret tolerances,
we could probably convince syminv() to treat A as if it had full rank.
Substituting zeros in the result also means that if you calculate A*syminv(A),
you can read the linear dependencies from the rows of the result. For
instance, if you got
+- -+
| 1 0 0 |
A * syminv(A) = | -1 0 1 |
| 0 0 1 |
+- -+
then that means that
x2 = -x1 + x3
> Do you use LU or QR decomposition?
No. LU is a popular and accurate way to to obtain inverses of square
matrices, but it does not exploit symmetry (i.e., LU can calculate the inverse
of nonsymmetric matrices). The more numerically intensive QR is even more
general, in that it can obtain inverses of nonsquare matrices (see equation
(2) above).
syminv() is implemented in terms of a sweep operator. Some researchers think
the sweep operator is not adequately accurate, but the sweep can be made
accurate if one puts the same effort in its development and implementation as,
say, LU decomposition.
What is sometimes taught as LU decomposition is not the LU decomposition
actually used. LU decomposition is often taught as finding L and U such
that A = LU, where L is lower triangular and U and upper triangular.
Numerically speaking, however, finding L and U such that A = LU cannot be done
very accurately, and so what is actually used is A = PLU, where P is a
permutation matrix. The reordering performed by P makes all the difference,
accuracywise, but it complicates how one uses the results from the
decomposition.
Similarly, the sweep operator, performed in its simplest form, is inaccurate.
There are two issues, scaling and ordering, the latter being similar to the
reordering performed by P in LU decomposition. To handle scaling, syminv()
scales the original matrix by the median of the diagonal. As for reordering,
that cannot be handled by pivot methods (permuation matrices) without losing
the symmetric structure of the problem. The counterpart is to determine the
order in which the rows are to be swept rather than sweeping them in order in
which they appear, and syminv() makes calculations concerning the roundoff
error that are very similar to those made in constructing P.
-- Bill
[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/
