[GenABEL-dev] faster polygenic

[Yurii Aulchenko]

One extra thought on how we can seep up 'polygenic' a bit (thanks to
Gulya to provoke these thoughts)

see line 23 of 'polylik'

ginvsig <- ervec %*% diag(es,ncol=length(qt)) %*% t(ervec)

this is very stupid way to do this, as time for computations goes
quadratic because of vector-matrix-vector product -- see example

[code]
nseq <- c(2,seq(1000,7500,500))
tt <- c()
i <- 1
for (N in nseq) {
v1 <- rnorm(N)
diagel <- rnorm(N)
sgm <- diag(diagel)
tt[i] <- system.time(t(v1) %*% sgm %*% v1)[3]
i<-i+1
}
plot(nseq,tt)
[/code]

but because the matrix has only diagonal elements non-zero, we can get
to same results with vec-vec-vec product, which goes much faster:

[code]
nseq <- c(2,seq(1000,7500,500))
tt <- c()
ttN <- c()
i <- 1
for (N in nseq) {
v1 <- rnorm(N)
diagel <- rnorm(N)
sgm <- diag(diagel)
tt[i] <- system.time(res <- t(v1) %*% sgm %*% v1)[3]
ttN[i] <- system.time(resB <- sum(v1*diagel*v1))[3]
print(c(N,res-resB))
i<-i+1
}
plot(nseq,tt,col="red")
points(nseq,ttN,col="green")
[/code]

will try to fix this line + run the tests -- probably we get x1.5
speed-up or so :)



another place where we can gain a bit is probably line 15:

sigma <- h2*tvar*relmat + (1-h2)*tvar*diag(x=1,ncol=nids,nrow=nids)

which is better re-written as

sigma <- h2*tvar*relmat + diag(x=(1-h2)*tvar,ncol=nids,nrow=nids)


Yurii