Anwenden einer Funktion auf eine Distanzmatrix in R

Question

Diese Frage kam heute in der Mailingliste von manipulatr.

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

Ich umformuliere.

Bei einer Abstandsmatrix (berechnet mit dist) eine Funktion auf die Zeilen der Abstandsmatrix anwenden.

Code:

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) 

Das Problem ist, dass die Funktion von Zeile wenden Sie die gesamte Matrix (statt nur den unteren Dreiecksteil zu speichern, haben also es zu viel Speicher für große Matrizen verwendet es.. in meinem Computer nicht für Matrizen der Dimensionen ~ 10000.

Irgendwelche Ideen?

Answer 1

Meine Lösung ist, die Indizes des Abstandsvektors zu erhalten, da eine Reihe und die Größe der Matrix. ich habe diese von codeguru

Nach der Übersetzung in R, unter der Annahme, Indizes beginnen bei 1, und unter der Annahme, dass ein niedriger tri anstelle der oberen tri-Matrix ich bekam.
EDIT: Mit der vektorisierte Version von rcs beigetragen

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)]) 
    }) 
} 

Also, mit dem Beispiel:

sumd2 <- getsum(d) 

Diese viel langsamer als as.matrix für kleine Matrizen vor Vektorisierung war. Aber nur etwa 3x so langsam nach dem Vektorisieren. In einem Intel Core2Duo 2ghz die Anwendung der Summe von Zeile der Größe 10000 Matrix dauerte etwas über 100s. Die as.matrix-Methode schlägt fehl. Danke rcs!

Answer 2

Dies ist eine vektorisierte Version der Funktion noeq (entweder Argument i oder 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 auf einem 2,4 GHz Intel Core 2 Duo (Mac OS 10.6.1) durchgeführt werden:

> 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

Vor allem für alle, die das noch nicht gesehen haben, empfehle ich dringend reading this article on the r-wiki über Code-Optimierung.

Hier ist eine andere Version ohne ifelse zu verwenden (das eine relativ langsame Funktion ist):

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 
} 

und Timings auf meinem Laptop:

> 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 

Und lapply ist schneller als 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