# [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

```