From: Douglas Bates <dmbates_at_gmail.com>

Date: Fri 03 Jun 2005 - 19:45:40 GMT

R-devel@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Sat Jun 04 05:52:04 2005

Date: Fri 03 Jun 2005 - 19:45:40 GMT

On 6/3/05, Susanne Heim <susanne.heim@stat.uni-muenchen.de> wrote:

> Dear R developers,

*>
**> The trace of the hat matrix H~(n,n) is computed as follows:
**>
**> tr(H) = tr(BS^-1B') = tr(S^-1B'B) := tr(X) = sum(diag(X))
**>
**> with B~(n,p), S~(p,p).
**> Since p is of the order 10^3 but S is sparse I would like to employ
**> Taucs linear solver ( http://www.tau.ac.il/~stoledo/taucs/ ) on
**>
**> SX = B'B.
**>
**> (Further improvement by implying a looping over i=1,...,p, calling
**> taucs_linsolve(S, X[,i], (B'B)[,i]) and saving X[i,i] only is pending.)
**>
**> For this purpose I compiled the C code "hattrace.c" to a shared object
**> using:
**> gcc -g -Wall -I/usr/local/taucs/src -I/usr/local/taucs/build/linux -c
**> hattrace.c -o hattrace.o
**> gcc -g -L/usr/local/taucs/external/lib/linux
**> -L/usr/local/taucs/lib/linux -L/usr/local/lib -L/opt/gnome/lib
**> -L/usr/lib/R/lib -shared -fpic -o hattrace.so hattrace.o -ltaucs
**> -llapack -lf77blas -lcblas -latlas -lmetis -lm -lg2c -lR
**>
**> I tried the following test commands:
**> library(splines)
**> library(SparseM)
**> B <- splineDesign(knots = 1:10, x = 4:7)
**> D <- diff(diag(dim(B)[2]), differences = 1)
**> BB <- t(B) %*% B
**> S <- as.matrix.ssc(BB + t(D) %*% D)
**> if (!is.loaded(symbol.C("hattrace"))) { dyn.load(paste("hattrace",
**> .Platform$dynlib.ext, sep = "")) }
**> out <- 0
**> spur <- (.C("hattrace", as.double(as.vector(slot(S, "ra"))),
**> as.integer(as.vector(slot(S, "ja") - 1)),
**> as.integer(as.vector(slot(S, "ia") - 1)),
**> as.integer(dim(S)[1]), as.double(as.vector(BB)),
**> as.double(out), PACKAGE = "hattrace"))[[6]]
**>
**> Unfortunately, I get an R process segmentation fault although the C Code
**> outputs the correct trace value to /tmp/hattrace.log which I checked by
**> a equivalent R routine. Since this segmentation fault does not occur
**> every time, I assume a pointer problem. Any help on how to solve it is
**> greatly appreciated.
*

Is S positive definite? If so, it may be more effective to take the Cholesky decomposition of S and solve the system S^(1/2)X = B then take the sum of the squares of the elements of X.

If you wish to provide me off-list with examples of the matrices S and B, I can check how best to do this with the Matrix package.

R-devel@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Sat Jun 04 05:52:04 2005

*
This archive was generated by hypermail 2.1.8
: Mon 24 Oct 2005 - 22:26:57 GMT
*