[Rd] Cholesky update/downdate

[Yves Deville]

Dear R-devel members,

I am looking for a fast Cholesky update/downdate. The matrix A being 
symmetric positive definite (n, n) and factorized as
A = L %*% t(L), the goal is to factor the new matrix  A +- C %*% t(C) 
where C is (n, r). For instance, C is 1-column when adding/removing an 
observation in a linear regression. Of special interest is the case 
where A is sparse.


Looking at the 'Matrix' package (help and source code), it seems that 
the CHOLMOD library shipped with 'Matrix' allows this,
but is not (yet?) interfaced in 'Matrix', where the 'update' method for 
Cholesky decomposition objects seems limited to a new matrix A + m*I 
with a scalar (diagonal) modification.


If this is true: are there plans to implement such up/downdates?

Best,

Yves

Yves Deville, statistical consultant, France.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Re: [Rd] Cholesky update/downdate

[Yves Deville]

Hi Douglas,

thanks for your answer.

My question indeed arises from a sparse matrix context: 'A' is sparse 
symmetric, and 'C' must also be sparse since it would otherwise fill.


It comes from a Bayes regression with a very large number of parameters; 
a loop over blocks will do the job in my specific case. Yet  I wondered 
about this since similar need for "covariance updating" may arise from 
linear filtering or kriging.


Douglas Bates wrote

The CHOLMOD library provides sparse matrix methods, especially the
Cholesky decomposition and modifications to that decomposition, which
is where the name comes from.  Do you expect to work with sparse
matrices?

I haven't seem too much code for update/downdate operations on the
Cholesky decomposition for dense matrices.  There were rank-1
update/downdate methods in Linpack but they didn't make it through to
Lapack.


__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

[Rd] ordering in 'gnls' with 'corCompSymm' corStruct

[Yves Deville]

Dear R-devel members,

While writing a new correlation structure similar to 'corCompSymm' and 
intended to be used with 'gnls', I got puzzled with the 'Initialize' method.


Using 'Initialize' before 'gnls' may be regarded as a mean to set an 
initial value for the corStruct parameter. However 'gnls' does not work 
properly with a 'corCompSymm' correlation structure when the data frame 
is not suitably ordered by group and when 'Initialize' is used for the 
corStruct before calling 'gnls'. Maybe the documentation could warn 
about this problem?


As an example, the results of 'fit1' and 'fit2' below are different, and 
only 'fit1' gives the good results. The problem remains if we drop the 
0.2  (for the 'value' formal) when calling Initialize.


My explanation is that when 'Initialize' is used with a 'corCompSymm' 
object and a data frame which is not sorted by the group variable, the 
'start' indices computed by the 'Dim' method are misleading. These 
'start' values will remain attached to the object in 'gnls' even though 
the data frame is sorted therein, since the call to 'Initialize' within 
'gnls' will not then change them. If this explanation is true, 
'Initialize.corCompSymm' should not accept an unordered 'data' formal, 
or should sort the data before computing 'start' with 'match'.


Best

Yves Deville, statistical consultant. France

##---
## generate grouped data
M <- 20
set.seed(123)
Ni <- 1 + rpois(M, lambda = 4); N <- sum(Ni)
grp <- factor(rep(1:M, times = Ni))
x <- lapply(Ni, function(n) 1:n)
y <- lapply(x, function(x) 1 + 2*x + rnorm(1) + rnorm(length(x)))
dataOrd <- data.frame(grp = grp, x = unlist(x), y = unlist(y))
## change order
ind <- sample(1:N); dataUnOrd <- dataOrd[ind, ]

library(nlme)
cS1 <- corCompSymm(0.2, form = ~ 1 | grp)
fit1 <- gnls(model = y ~ a + b*x, data = dataUnOrd, start = c(a = 0.1, b 
= 0.3),
correlation = cS1) 
## the same with 'Initialize'

cS2 <- corCompSymm(0.2, form = ~ 1 | grp)
cS2 <- Initialize(cS2, data = dataUnOrd)
fit2 <- gnls(model = y ~ a + b*x, data = dataUnOrd, start = c(a = 0.1, b 
= 0.3),
correlation = cS2)  


## Dim does not give the good 'start' on 'cS2'
cS3 <- corCompSymm(0.2, form = ~ 1 | grp)
cS3 <- Initialize(cS3, data = dataOrd)
Dim(cS2)$start
Dim(cS3)$start
##--

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel