2013-06-13 20 views
24

Ho una matrice m e un vettore v. Vorrei moltiplicare la prima colonna della matrice m dal primo elemento del vettore v e moltiplicare la seconda colonna della matrice m dal secondo elemento del vettore v e così via. Posso farlo con il seguente codice, ma sto cercando un modo che non richiede le due chiamate di trasposizione. Come posso farlo più velocemente in R?Il modo più veloce per moltiplicare le colonne di matrice con elementi vettoriali in R

m <- matrix(rnorm(120000), ncol=6) 
v <- c(1.5, 3.5, 4.5, 5.5, 6.5, 7.5) 

system.time(t(t(m) * v)) 

# user system elapsed 
# 0.02 0.00 0.02 
+0

correlati: http://stackoverflow.com/q/3643555/946850 – krlmlr

risposta

33

uso di alcune algebra lineare ed eseguire moltiplicazione di matrici, che è abbastanza veloce in R.

esempio

m %*% diag(v)

po 'di benchmarking

m = matrix(rnorm(1200000), ncol=6) 

v=c(1.5, 3.5, 4.5, 5.5, 6.5, 7.5) 
library(microbenchmark) 
microbenchmark(m %*% diag(v), t(t(m) * v)) 
## Unit: milliseconds 
##   expr  min  lq median  uq  max neval 
## m %*% diag(v) 16.57174 16.78104 16.86427 23.13121 109.9006 100 
##  t(t(m) * v) 26.21470 26.59049 32.40829 35.38097 122.9351 100 
+0

Proprio così, dovrebbe essere il microbenchmark (m% *% diag (v), t (t (m) * v)) – rose

+0

Effettivamente, Modifica apportata @rose – mnel

+1

Ho scoperto che il risultato dipende in gran parte dalla lunghezza 'V. Per 'v' più corti, l'opzione' diag() 'è più veloce, ma alla fine vince la doppia trasposizione. – krlmlr

3

Come sottolinea @Arun, non so che batterete la vostra soluzione in termini di efficienza temporale. In termini di codice comprensibilità, ci sono altre opzioni però:

Una possibilità:

> mapply("*",as.data.frame(m),v) 
     V1 V2 V3 
[1,] 0.0 0.0 0.0 
[2,] 1.5 0.0 0.0 
[3,] 1.5 3.5 0.0 
[4,] 1.5 3.5 4.5 

E un altro:

sapply(1:ncol(m),function(x) m[,x] * v[x]) 
+0

Dubito che questo sarò essere più veloce di lavorare su matrici (la tua prima soluzione, in particolare). – Arun

+0

quando controllo il system.time per un grande campione non c'è differenza tra loro e non è più veloce. – rose

+0

@rose - Nonostante fornisca alternative, sono d'accordo con Arun qui. Non sono sicuro di cosa non funzioni con la soluzione 't (t (..' – thelatemail

15

Se si dispone di un maggior numero di colonne vostra t (t (m) * v) supera la soluzione di matrice moltiplicazione per una vasta margine. Tuttavia, esiste una soluzione più veloce, ma ha un costo elevato nell'utilizzo della memoria. Crei una matrice grande quanto m usando rep() e moltiplica elementwise. Ecco il confronto, modificando l'esempio di mnel:

m = matrix(rnorm(1200000), ncol=600) 
v = rep(c(1.5, 3.5, 4.5, 5.5, 6.5, 7.5), length = ncol(m)) 
library(microbenchmark) 

microbenchmark(t(t(m) * v), 
    m %*% diag(v), 
    m * rep(v, rep.int(nrow(m),length(v))), 
    m * rep(v, rep(nrow(m),length(v))), 
    m * rep(v, each = nrow(m))) 

# Unit: milliseconds 
#         expr  min   lq  mean  median   uq  max neval 
#       t(t(m) * v) 17.682257 18.807218 20.574513 19.239350 19.818331 62.63947 100 
#       m %*% diag(v) 415.573110 417.835574 421.226179 419.061019 420.601778 465.43276 100 
# m * rep(v, rep.int(nrow(m), ncol(m))) 2.597411 2.794915 5.947318 3.276216 3.873842 48.95579 100 
#  m * rep(v, rep(nrow(m), ncol(m))) 2.601701 2.785839 3.707153 2.918994 3.855361 47.48697 100 
#    m * rep(v, each = nrow(m)) 21.766636 21.901935 23.791504 22.351227 23.049006 66.68491 100 

Come si può vedere, l'utilizzo di "ogni" in rep() sacrifici velocità per chiarezza. La differenza tra rep.int e rep sembra essere trascurabile, entrambe le implementazioni scambiano i posti su ripetute serie di microbenchmark(). Tieni presente che ncol (m) == lunghezza (v).

autoplot

+0

Si noti che la doppia trasposizione duplica anche la matrice almeno una volta, non è sicuro che l'utilizzo della memoria sia molto meglio della semplice espansione della matrice, L'espansione stessa può essere eseguita leggermente più leggibile usando 'matrix (v, nrow = nrow (m) , ncol = ncol (m), byrow = TRUE) '. – krlmlr

+0

Informazioni sulla soluzione 'rep' che scrivi" ... ha un costo elevato nell'utilizzo della memoria ". Non 't (m)' comporta lo stesso costo, perché questo crea una nuova matrice con lo stesso numero di elementi di 'm'? – jochen

1

Come fatto da bluegrue, un semplice rep sarebbe sufficiente anche per effettuare la moltiplicazione elemento saggio.

Il numero di moltiplicazioni e sommatorie viene ridotto di un margine ampio come se si eseguisse una semplice moltiplicazione di matrice con diag(), dove per questo caso si possono evitare molte moltiplicazioni di zero.

m = matrix(rnorm(1200000), ncol=6) 
v=c(1.5, 3.5, 4.5, 5.5, 6.5, 7.5) 
v2 <- rep(v,each=dim(m)[1]) 
library(microbenchmark) 
microbenchmark(m %*% diag(v), t(t(m) * v), m*v2) 

Unit: milliseconds 
      expr  min  lq  mean median  uq  max neval cld 
m %*% diag(v) 11.269890 13.073995 16.424366 16.470435 17.700803 95.78635 100 b 
    t(t(m) * v) 9.794000 11.226271 14.018568 12.995839 15.010730 88.90111 100 b 
     m * v2 2.322188 2.559024 3.777874 3.011185 3.410848 67.26368 100 a 
1

Per ragioni di completezza, ho aggiunto sweep al benchmark. Nonostante i suoi nomi degli attributi un po 'fuorviante, penso che potrebbe essere più leggibile rispetto ad altre alternative, e anche abbastanza veloce:

n = 1000 
M = matrix(rnorm(2 * n * n), nrow = n) 
v = rnorm(2 * n) 

microbenchmark::microbenchmark(
    M * rep(v, rep.int(nrow(M), length(v))), 
    sweep(M, MARGIN = 2, STATS = v, FUN = `*`), 
    t(t(M) * v), 
    M * rep(v, each = nrow(M)), 
    M %*% diag(v) 
) 

Unit: milliseconds 
             expr   min   lq  mean 
    M * rep(v, rep.int(nrow(M), length(v))) 5.259957 5.535376 9.994405 
sweep(M, MARGIN = 2, STATS = v, FUN = `*`) 16.083039 17.260790 22.724433 
           t(t(M) * v) 19.547392 20.748929 29.868819 
       M * rep(v, each = nrow(M)) 34.803229 37.088510 41.518962 
           M %*% diag(v) 1827.301864 1876.806506 2004.140725 
     median   uq  max neval 
    6.158703 7.606777 66.21271 100 
    20.479928 23.830074 85.24550 100 
    24.722213 29.222172 92.25538 100 
    39.920664 42.659752 106.70252 100 
1986.152972 2096.172601 2432.88704 100 
Problemi correlati