Zastosowanie funkcji do macierzy odległości w R

Question

To pytanie pojawiło się dzisiaj na liście dyskusyjnej manipulatr.

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

Przepisuję.

Mając matrycę odległości (obliczoną za pomocą dist) zastosuj funkcję do rzędów macierzy odległości.

Kod:

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) 

Problem polega na tym, aby zastosować funkcję rzędzie trzeba przechowywać całą matrycę (zamiast tylko w dolnej części trójkątnej więc używa zbyt dużo pamięci dla dużych matrycach It.. nie w moim komputerze dla macierzy o wymiarach ~ 10000.

Jakieś pomysły?

Answer 1

Moje rozwiązanie jest dostać indeksy wektora odległości, ponieważ rząd i wielkość matrycy. mam to od 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; 
} 

Po przetłumaczeniu na R, przyjmując, że indeksy zaczynają się od 1, i zakładam, że otrzymałem niższą tri zamiast macierzy górnej tri.
EDIT: Używanie vectorized wersja wniesionego przez 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)]) 
    }) 
} 

Tak, z przykładu:

sumd2 <- getsum(d) 

To był znacznie wolniej niż as.matrix dla małych matrycach przed Wektoryzacja. Ale około 3 razy wolniej po wektoryzacji. W Intel Core2Duo 2ghz zastosowanie sumy po rzędzie matrycy wielkości 10000 zajęło nieco ponad 100 sekund. Metoda as.matrix nie działa. Dzięki rcs!

Answer 2

To vectorized wersja funkcji noeq (albo argumentu i lub j):

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 

Timings wykonywane są na 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

Przede wszystkim, dla każdego, kto nie widział jeszcze tego, gorąco polecam reading this article on the r-wiki o optymalizacji kodu.

Oto kolejna wersja bez użycia ifelse (to stosunkowo powolny funkcji):

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 
} 

A czasy na moim laptopie:

> 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 

I lapply jest szybszy niż 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