chol() dimnames

[Douglas Bates]

Peter Dalgaard BSA <p.dalgaard@biostat.ku.dk> writes:

> However: Somewhere in the back of my mind lies a note to the effect
> that chol() should really handle the rank-deficient cases better, via
> some sort of pivoting. Will this muck things up (possibly not if it
> only removes singularities and plugs in rows/cols of zeroes)? I seem to
> recall that the current chol() got lme() in trouble. 

I believe the trouble there was not due to R's chol not handling the
semi-definite case.  (Sorry about the double negative.)  I think it is
because internally chol is called differently than in S.  We have code
like

#ifdef R_S_H
	  F77_CALL(chol)(auxRes, &ncol, &ncol, auxRes, &l);
#else
	  F77_CALL(chol)(auxRes, &ncol, res, &zero, &zero, &l);
#endif /* R_S_H */

in the C code that implements lme.

I believe there is a pivoting option available for the Linpack dchdc
routine but not for the dpodc routine.  (My copy of the manual is at
the office so I can't check this.)  That pivoting is similar to the
pivoting in dqrdc - you specify which columns (and the corresponding
rows in the case of dchdc) are available to be re-ordered and it
re-orders everything in the optimal ordering.  This doesn't work well
for model matrices or quantities derived from model matrices where you
want to keep the columns associated with the same terms adjacent if
possible.

Notice that in your example the diagonal elements of the decomposition
are strictly decreasing after pivoting.  That is what I mean by the
optimal ordering.

> > chol(x,piv=T)
>          [,1]       [,2]       [,3]          [,4]          [,5] 
> [1,] 2.065148 -0.7414994 -0.8465894  1.583980e+00 -8.629924e-01
> [2,] 0.000000  1.4051462  0.8172503 -5.339291e-01  3.893657e-01
> [3,] 0.000000  0.0000000  0.5126111  3.296278e-01  4.426312e-01
> [4,] 0.000000  0.0000000  0.0000000  2.107342e-08  5.268356e-09
> [5,] 0.000000  0.0000000  0.0000000  0.000000e+00  7.850462e-17
> attr(, "pivot"):
> [1] 2 4 3 1 5
> attr(, "rank"):
> [1] 5
>     ^
>     !
> oops....
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._