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.
I am running R 2.0.0 on SuSE Linux.
Susanne Heim
--
Susanne Heim
Institute of Statistics, Ludwig-Maximilians-University
Ludwigstr.33/II, D-80539 Munich, Germany.
Phone: +49-89-2180-2226
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 103 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.
I am running R 2.0.0 on SuSE Linux.
Susanne Heim
--
Susanne Heim
Institute of Statistics, Ludwig-Maximilians-University
Ludwigstr.33/II, D-80539 Munich, Germany.
Phone: +49-89-2180-2226
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.