Indicizzazione degli elementi di una matrice in R

Question

Il problema è piuttosto stupido, ma mi chiedo se mi manca qualcosa. Diciamo che c'è un vettore k che contiene alcuni numeri, dire

> k 
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 

voglio trasformare questo ad una matrice

> m 
    [,1] [,2] [,3] [,4] [,5] 
[1,] 1 2 3 4 5 
[2,] 0 6 7 8 9 
[3,] 0 0 10 11 12 
[4,] 0 0 0 13 14 
[5,] 0 0 0 0 15 

La mia prima idea era quella di utilizzare qualcosa con upper.tri(), ad esempio come m[upper.tri(m, diag = TRUE)] <- k, ma questo non darà la matrice sopra.

C'è una soluzione più intelligente a questo? Di seguito c'è la mia soluzione, ma diciamo che non ne sono troppo orgoglioso.

---

rows <- rep(1:5, 5:1) 

cols1 <- rle(rows)$lengths 


cols <- do.call(c, lapply(1:length(cols1), function(x) x:5)) 

for(i in 1:length(k)) { 
    m[rows[i], cols[i]] <- k[i] 
} 

Answer 1

Una variante di risposta @docendodiscimus': Invece di recepimento è possibile modificare gli indici di riga e di Col, che si ottiene avvolgendo lower.tri in which:

n = 5 
m = matrix(0, n, n) 

m[ which(lower.tri(m, diag=TRUE), arr.ind=TRUE)[, 2:1] ] = seq(sum(seq(n))) 


    [,1] [,2] [,3] [,4] [,5] 
[1,] 1 2 3 4 5 
[2,] 0 6 7 8 9 
[3,] 0 0 10 11 12 
[4,] 0 0 0 13 14 
[5,] 0 0 0 0 15 

Per capire come funziona, guardate il lato sinistro in fasi:

• lower.tri(m, diag=TRUE)
• which(lower.tri(m, diag=TRUE), arr.ind=TRUE)
• which(lower.tri(m, diag=TRUE), arr.ind=TRUE)[, 2:1]

Credo che recepiscono potrebbe essere costoso se la matrice è di grandi dimensioni, che è il motivo per cui mi piacerebbe prendere in considerazione questa opzione. Nota: la risposta di Joseph Wood suggerisce che ho torto, dal momento che il modo di trasporre è più veloce nel suo benchmark.

---

(Grazie a @JosephWood :) Invece di enumerare e sommando con sum(seq(n)), è possibile utilizzare (n^2 - n)/2 + n.

Answer 2

Ecco un'opzione con lower.tri e t di trasporre il risultato:

k <- 1:15 
m <- matrix(0, 5,5) 
m[lower.tri(m, diag = TRUE)] <- k 
m <- t(m) 
m 
#  [,1] [,2] [,3] [,4] [,5] 
#[1,] 1 2 3 4 5 
#[2,] 0 6 7 8 9 
#[3,] 0 0 10 11 12 
#[4,] 0 0 0 13 14 
#[5,] 0 0 0 0 15 

---

microbenchmark

Dato che c'era un po 'di confusione con riferimento di Giuseppe, ecco un altro uno. Ho testato le tre soluzioni per matrici di dimensioni 10 * 10; 100 * 100; 1000 * 1000; 10000 * 10000.

Risultati:

A quanto pare, la performance dipende fortemente dalla dimensione della matrice. Per le matrici di grandi dimensioni, la risposta di Joseph è più veloce, mentre per le matrici più piccole, la mia è stata l'approccio più veloce. Si noti che ciò non tiene conto dell'efficienza della memoria.

riferimento riproducibile:

Joseph <- function(k, n) { 
    y <- 1L 
    t <- rep(0L,n) 
    j <- c(y, sapply(1:(n-1L), function(x) y <<- y+(n+1L)-x)) 
    t(vapply(1:n, function(x) c(rep(0L,x-1L),k[j[x]:(j[x]+n-x)]), t, USE.NAMES = FALSE)) 
} 

Frank <- function(k, n) { 
    m = matrix(0L, n, n) 
    m[ which(lower.tri(m, diag=TRUE), arr.ind=TRUE)[, 2:1] ] = k 
    m 
} 

docendo <- function(k,n) { 
    m <- matrix(0L, n, n) 
    m[lower.tri(m, diag = TRUE)] <- k 
    t(m) 
} 

library(microbenchmark) 
library(data.table) 
library(ggplot2) 
n <- c(10L, 100L, 1000L, 10000L) 
k <- lapply(n, function(x) seq.int((x^2 + x)/2)) 

b <- lapply(seq_along(n), function(i) { 
    bm <- microbenchmark(Joseph(k[[i]], n[i]), Frank(k[[i]], n[i]), docendo(k[[i]], n[i]), times = 10L) 
    bm$n <- n[i] 
    bm 
}) 

b1 <- rbindlist(b) 

ggplot(b1, aes(expr, time)) + 
    geom_violin() + 
    facet_wrap(~ n, scales = "free_y") + 
    ggtitle("Benchmark for n = c(10L, 100L, 1000L, 10000L)") 

Controllo parità di risultati:

all.equal(Joseph(k[[1]], n[1]), Frank(k[[1]], n[1])) 
#[1] TRUE 
all.equal(Joseph(k[[1]], n[1]), docendo(k[[1]], n[1])) 
#[1] TRUE 

Nota: non ho incluso l'approccio di George nel confronto dal momento che, a giudicare dai risultati di Giuseppe , sembra essere molto più lento. Quindi tutti gli approcci comparati nel mio benchmark sono scritti solo nella base R.

Answer 3

library(miscTools) 
k <- 1:15 
triang(k, 5) 

Answer 4

Ecco una soluzione di base R veramente veloce:

Aggiornamento

Ho modificato il codice leggermente in modo che io chiamo solo vapply una volta al posto del sapply/vapply combo che avevo prima (ho anche sbarazzarsi di USE.NAMES=FALSE in quanto sembra non fare alcuna differenza). Anche se questo è un po 'più pulito, non ha cambiato significativamente i tempi sul mio computer (ho rirquisito i benchmark di docendo con i grafici e sembra essere quasi lo stesso).

Triangle1 <- function(k,n) { 
    y <- -n 
    r <- rep(0L,n) 
    t(vapply(1:n, function(x) {y <<- y+n+2L-x; c(rep(0L,x-1L),k[y:(y+n-x)])}, r)) 
} 

Ecco alcuni tempi:

Triangle2 <- function(k,n) { 
    m <- matrix(0, n,n) 
    m[lower.tri(m, diag = TRUE)] <- k 
    t(m) 
} 

Triangle3 <- function(k, n) { 
    m = matrix(0, n, n) 
    m[ which(lower.tri(m, diag=TRUE), arr.ind=TRUE)[, 2:1] ] = k ## seq(sum(seq(n))) for benchmarking 
    m 
} 

k2 <- 1:50005000 
n2 <- 10^4 

system.time(t1 <- Triangle1(k2,n2)) 
user system elapsed   ## previously user system elapsed 
2.29 0.08 2.41   ##    2.37 0.13 2.52 

system.time(t2 <- Triangle2(k2,n2)) 
user system elapsed 
5.40 0.91 6.30 

system.time(t3 <- Triangle3(k2,n2)) 
user system elapsed 
7.70 1.03 8.77 

system.time(t4 <- triang(k2,n2)) 
user system elapsed 
433.45 0.20 434.88 

Una cosa che è un po 'sconcertante per me è che l'oggetto prodotto da Triangle1 è grande la metà di tutte le altre soluzioni.

object.size(t1) 
400000200 bytes 

object.size(t2) ## it's the same for t3 and t4 
800000200 bytes 

Quando eseguo alcuni controlli, diventa più confuso.

all(sapply(1:ncol(t1), function(x) all(t1[,x]==t2[,x]))) 
[1] TRUE 

class(t1) 
[1] "matrix" 
class(t2) 
[1] "matrix" 

attributes(t1) 
$dim 
[1] 10000 10000 
attributes(t2) 
$dim 
[1] 10000 10000 

## not sure what's going on here 
identical(t1,t2) 
[1] FALSE 

identical(t2,t3) 
[1] TRUE 

---

Come @Frank sottolineato nei commenti, t1 è una matrice integer mentre gli altri sono numerici. Avrei dovuto sapere che uno dei most important R functions mi avrebbe detto questa informazione fin dall'inizio.

str(t1) 
int [1:10000, 1:10000] 1 0 0 0 0 0 0 0 0 0 ... 
str(t2) 
num [1:10000, 1:10000] 1 0 0 0 0 0 0 0 0 0 ...