[Rd] (PR#3427)
maechler at stat.math.ethz.ch
maechler at stat.math.ethz.ch
Tue Jul 8 09:47:23 MEST 2003
>>>>> "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
More information about the R-devel
mailing list