# [R] Problem with new(er) R version's matrix package

Martin Maechler maechler at stat.math.ethz.ch
Mon Apr 28 13:06:24 CEST 2014

```>>>>> Werner W <pensterfuzzer at yahoo.de>
>>>>>     on Sun, 27 Apr 2014 17:48:23 +0100 writes:

> Dear Martin, Arne, David,

> Thanks for considering my problem. However, I have a bit of a problem in making a small reproducible example. So far, I tried to drill down into the code and to quick fix this by setting the particular tolerance to zero which switches off the check in ".solve.sparse.dgC". The obstacle I got stuck with doing that is, seemingly, the "solvetol" parameter from systemfit.control is not handed through all the way.
> ".calcGLS" from systemfit still hands it over but it does not reach ".solve.sparse.dgC"

> This is the stack at that point before the error is thrown:
> Browse[2]> where
> where 1: .solve.sparse.dgC(as(a, "dgCMatrix"), b = b, tol = tol)
> where 2: solve.dsC.dC(a, b)
> where 3: solve(a, as(b, "CsparseMatrix"))
> where 4: solve(a, as(b, "CsparseMatrix"))
> where 5: eval(expr, envir, enclos)
> where 6: eval(call, sys.frame(sys.parent()))
> where 7: callGeneric(a, as(b, "CsparseMatrix"))
> where 8: solve(forceCspSymmetric(a, isTri = FALSE), b = Diagonal(nrow(a)))
> where 9: solve(forceCspSymmetric(a, isTri = FALSE), b = Diagonal(nrow(a)))
> where 10: .local(a, b, ...)
> where 11: solve(W, tol = solvetol)
> where 12: solve(W, tol = solvetol)
> where 13: as.matrix(solve(W, tol = solvetol)[1:ncol(xMat), 1:ncol(xMat)])
> where 14: .calcGLS(xMat = xMatAll, R.restr = R.restr, q.restr = q.restr,
> sigma = rcov, validObsEq = validObsEq, useMatrix = control\$useMatrix,
> solvetol = control\$solvetol)
> where 15 at ILLS.R#465: systemfit(eqsys, method = "SUR", data = dat, restrict.matrix = restrmat,
> control = sysfit.contr)

> I will now further try to create a small, reproducible example.

Thank you, Werner,  the above is already very useful!

Currently, I'd concluded that there must be a buglet in either
systemfit or Matrix  as you assert that the 'tol' argument is
not passed all the way down the Matrix package sparse solver.
[after 10 more minutes ..]
My current guess is that the buglet is in 'Matrix'...
I see the phenomenon also, using   example(systemfit)
so we have a reproducible example for
"tol is not passed all the way to the Matrix solve(..)",
and I will investigate that later.
There are a few more urgent things currently, and other
duties...

> I also have to better understand what the procedure is doing to be able to decide whether the previous results were wrong or this newly occurring error is wrong.

Indeed, the above is the real question you should address and
you have to do yourself alone.

> Many thanks,
> Werner

> Martin Maechler <maechler at stat.math.ethz.ch> schrieb am 22:00 Samstag, 26.April 2014:

>>>>> Arne Henningsen <arne.henningsen at gmail.com>
>>
>>>>>>>      on Sat, 26 Apr 2014 08:15:37 +0200 writes:
>>
>>     > On 25 April 2014 20:15, David Winsemius
>>     > <dwinsemius at comcast.net> wrote:
>>     >>
>>     >> On Apr 25, 2014, at 9:17 AM, Werner W. wrote:
>>     >>
>>     >>> Dear Rs,
>>     >>>
>>     >>> I am re-executing some older code. It does work in the
>>     >>> ancient R 2.12.0 which I still have on my PC but with
>>     >>> the new version R 3.1.0 it does not work any more (but
>>     >>> some other new stuff, which won't work with 2.12).
>>     >>>
>>     >>> The problem arises in context with the systemfit package
>>     >>> using the matrix package. In R 3.1.0 the following error
>>     >>> is thrown: Error in as.matrix(solve(W, tol =
>>     >>> solvetol)[1:ncol(xMat), 1:ncol(xMat)]) : error in
>>     >>> evaluating the argument 'x' in selecting a method for
>>     >>> function 'as.matrix': Error in .solve.sparse.dgC(as(a,
>>     >>> "dgCMatrix"), b = b, tol = tol) : LU computationally
>>     >>> singular: ratio of extreme entries in |diag(U)| =
>>     >>> 7.012e-39
>>     >>>
>>     >>> there some change in the defaults of the matrix package?
>>     >>> I couldn't find anything apparent in the changelog. As
>>     >>> the same code works in R 2.12.0, I suppose that the
>>     >>> problem is not my data.
>>     >>
>>     >> You have not told us what version of the Matrix package
>>     >> you were using.  As such I would suggest that you review
>>     >> the Changelog which is a link for the CRAN page for
>>     >> pkg:Matrix and go back 4 years or so since R major
>>     >> versions change about once a year.
>>     >>
>>     >> http://cran.r-project.org/web/packages/Matrix/ChangeLog
>>
>>     > reproducible example.
>>
>> Yes, please do.   As maintainer of the Matrix package, I'm
>> willing to look into the situation of course.
>>
>> As was mentioned, many things have changed in 4 years.
>> The error message above looks like you'd want to invert a
>> (very close to) singular matrix, and there could be quite few
>> reasons why parts of the older code gave slightly different
>>
>> Without a reproducible example, we can't get started though.
>>
>> Best regards,
>> Martin Maechler, ETH Zurich
>>
>>
>>
>>

```