Aplicar una función a una matriz de distancia en R

Question

Esta pregunta ha llegado hoy en la lista de distribución de manipulatr.

http://groups.google.com/group/manipulatr/browse_thread/thread/fbab76945f7cba3f 

Estoy reformulando.

Dada una matriz de distancia (calculada con dist) aplique una función a las filas de la matriz de distancias.

Código:

library(plyr) 
N <- 100 
a <- data.frame(b=1:N,c=runif(N)) 
d <- dist(a,diag=T,upper=T) 
sumd <- adply(as.matrix(d),1,sum) 

El problema es que para aplicar la función por la fila que tiene que almacenar toda la matriz (en lugar de sólo la parte triangular inferior por lo que utiliza demasiada memoria para grandes matrices Es.. fracasa en mi equipo para matrices de dimensiones ~ 10000.

¿Alguna idea?

Answer 1

mi solución es conseguir que los índices del vector de distancia, dieron una fila y el tamaño de la matriz. tengo esto desde codeguru

int Trag_noeq(int row, int col, int N) 
{ 
    //assert(row != col); //You can add this in if you like 
    if (row<col) 
     return row*(N-1) - (row-1)*((row-1) + 1)/2 + col - row - 1; 
    else if (col<row) 
     return col*(N-1) - (col-1)*((col-1) + 1)/2 + row - col - 1; 
    else 
     return -1; 
} 

Después de traducir a R, suponiendo que los índices comiencen en 1, y asumiendo un tri inferior en lugar de tri matrix superior que obtuve.
EDIT: El uso de la versión vectorizada aportado por rcs

noeq.1 <- function(i, j, N) { 
    i <- i-1 
    j <- j-1 
    ix <- ifelse(i < j, 
       i*(N-1) - (i-1)*((i-1) + 1)/2 + j - i, 
       j*(N-1) - (j-1)*((j-1) + 1)/2 + i - j) * ifelse(i == j, 0, 1) 
    ix 
} 

## To get the indexes of the row, the following one liner works: 

getrow <- function(z, N) noeq.1(z, 1:N, N) 

## to get the row sums 

getsum <- function(d, f=sum) { 
    N <- attr(d, "Size") 
    sapply(1:N, function(i) { 
     if (i%%100==0) print(i) 
     f(d[getrow(i,N)]) 
    }) 
} 

lo tanto, con el ejemplo:

sumd2 <- getsum(d) 

Este fue mucho más lento que as.matrix para las pequeñas matrices antes de la vectorización. Pero casi 3 veces más lenta después de vectorizar. En un Intel Core2Duo 2ghz aplicando la suma por fila del tamaño 10000, la matriz tardó algo más de 100 segundos. El método as.matrix falla. Gracias rcs!

Answer 2

Ésta es una versión vectorizada de la función noeq (alguno de los argumentos o ij):

noeq.1 <- function(i, j, N) { 
    i <- i-1 
    j <- j-1 
    ifelse(i < j, 
      i*(N-1) - ((i-1)*i)/2 + j - i, 
      j*(N-1) - ((j-1)*j)/2 + i - j) * ifelse(i == j, 0, 1) 
} 

> N <- 4 
> sapply(1:N, function(i) sapply(1:N, function(j) noeq(i, j, N))) 
    [,1] [,2] [,3] [,4] 
[1,] 0 1 2 3 
[2,] 1 0 4 5 
[3,] 2 4 0 6 
[4,] 3 5 6 0 
> sapply(1:N, function(i) noeq.1(i, 1:N, N)) 
    [,1] [,2] [,3] [,4] 
[1,] 0 1 2 3 
[2,] 1 0 4 5 
[3,] 2 4 0 6 
[4,] 3 5 6 0 

sincronizaciones se realizan en un 2,4 GHz Intel Core 2 Duo (Mac OS 10.6.1):

> N <- 1000 
> system.time(sapply(1:N, function(j) noeq.1(1:N, j, N))) 
    user system elapsed 
    0.676 0.061 0.738 
> system.time(sapply(1:N, function(i) sapply(1:N, function(j) noeq(i, j, N)))) 
    user system elapsed 
14.359 0.032 14.410 

Answer 3

En primer lugar, para cualquier persona que no haya visto esto todavía, recomiendo reading this article on the r-wiki sobre la optimización de código.

Aquí es otra versión sin usar ifelse (que es una función relativamente lenta):

noeq.2 <- function(i, j, N) { 
    i <- i-1 
    j <- j-1 
    x <- i*(N-1) - (i-1)*((i-1) + 1)/2 + j - i 
    x2 <- j*(N-1) - (j-1)*((j-1) + 1)/2 + i - j 
    idx <- i < j 
    x[!idx] <- x2[!idx] 
    x[i==j] <- 0 
    x 
} 

y los tiempos en mi portátil:

> N <- 1000 
> system.time(sapply(1:N, function(i) sapply(1:N, function(j) noeq(i, j, N)))) 
    user system elapsed 
    51.31 0.10 52.06 
> system.time(sapply(1:N, function(j) noeq.1(1:N, j, N))) 
    user system elapsed 
    2.47 0.02 2.67 
> system.time(sapply(1:N, function(j) noeq.2(1:N, j, N))) 
    user system elapsed 
    0.88 0.01 1.12 

Y lapply es más rápido que sapply:

> system.time(do.call("rbind",lapply(1:N, function(j) noeq.2(1:N, j, N)))) 
    user system elapsed 
    0.67 0.00 0.67