[Rd] (PR#3427)

[maechler at stat.math.ethz.ch]

>>>>> "PD" == Peter Dalgaard BSA <p.dalgaard at biostat.ku.dk>
>>>>>     on 07 Jul 2003 23:59:59 +0200 writes:

    PD> Dursun.Bulutoglu at afit.edu writes:
    >> Hi;
    >> 
    >> I am having problems inverting matrices using the function
    >> solve()
    >> 
    >> For example R can not invert the following matrix
    >> 
    >> =20
    >> 
    >> [,1]           [,2]              [,3]              [,4]
    >> [,5]
    >> 
    >> [1,]      25            500              11250           275000
    >> 7.106250e+06
    >> 
    >> [2,]      500          11250           275000         7106250
    >> 1.906250e+08
    >> 
    >> [3,]      11250       275000         7106250       190625000
    >> 5.247656e+09
    >> 
    >> [4,]      275000     7106250       190625000     5247656250
    >> 1.471719e+11
    >> 
    >> [5,]      7106250   190625000    5247656250   147171875000
    >> 4.184754e+12
    >> 
    >> solve(t(xxmodel)%*%(xxmodel))
    >> 
    >> Yields the following massage:
    >> 
    >> Error in solve.default(t(xxmodel) %*% (xxmodel)) : singular matrix `a'
    >> in solve
    >> 
    >> The above 5X5 matrix is invertible. It has non-zero eigenvalues. Could
    >> someone explain whether there is a problem in R's solve() function.

    PD> Please use the r-help list unless you are sure that you have found a
    PD> bug (and check up on the definition of a bug in the FAQ, and also
    PD> details required when reporting).

yes!!

    PD> You seem to be running into a generic problem with matrix inverters,
    PD> namely that they are unhappy when the columns are on widely different
    PD> scales. Your diagonal elements vary by a factor of about 2e11. Since
    PD> detection of linear dependency has to operate with a small tolerance
    PD> to account for roundoff, this can become indistinguishable from a
    PD> nonsingular matrix with a large range of eigenvalues.

    PD> The typical workaround would to rescale the columns of xxmodel.

2 more things:

1) my version of R  *does* invert your example matrix.
   Probably because newer versions of R have been using LAPACK
   as opposed to LINPACK.

2) If you do use a pre-LAPACK version of solve(), i.e., solve.default();
   there's an argument `tol = 1e-7' which allows to set
   a tolerance about how singularity is determined.

Martin