Индексирование элементов матрицы в R

Question

Проблема довольно глупа, но мне интересно, не хватает ли я чего-то. Допустим, что существует вектор k, который содержит некоторые цифры, скажем

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

Я хочу, чтобы преобразовать это в матрицу

> 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 

Моя первая идея была использовать что-то с upper.tri(), например, как m[upper.tri(m, diag = TRUE)] <- k, но это не даст матрицу выше.

Есть ли более интеллектуальное решение? Ниже есть мое решение, но давайте просто скажем, что я не слишком горжусь этим.

---

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

Вариантом ответа @docendodiscimus': Вместо транспонирования вы можете изменить строки и Col индексы, которые вы получаете, обернув lower.tri в 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 

Чтобы понять, как это работает, посмотрите на левая рука шагов:

• 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]

Я думаю, перенося может быть дорогостоящим, если матрица велика, поэтому я бы рассмотреть этот вариант. Примечание. Ответ Джозефа Вуда говорит о том, что я ошибаюсь, поскольку в его контрольном плане скорость переноса намного быстрее.

---

(Благодаря @JosephWood :) Вместо перечисления и суммирования с sum(seq(n)), вы можете использовать (n^2 - n)/2 + n.

Answer 2

Вот вариант с использованием lower.tri и t транспонировать результат:

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

Поскольку существует некоторая путаница с эталоном Иосифа, вот еще один. Я проверил три решения для матриц размером 10 * 10; 100 * 100; 1000 * 1000; 10000 * 10000.

Результаты:

Очевидно, что производительность сильно зависит от размера матрицы. Для больших матриц ответ Джозефа работает быстрее, а для меньших матриц мой - самый быстрый подход. Обратите внимание, что это не учитывает эффективность памяти.

Возпроизводимо тест:

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)") 

Проверить равенство результатов:

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 

Примечание: я не включал подход Джорджа в сравнении, так как, судя по результатам Иосифа , это, кажется, намного медленнее. Таким образом, все подходы, сравниваемые в моем тесте, написаны только в базе R.

Answer 3

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

Answer 4

Вот очень быстро база R решение:

Update

Я немного изменил код так, что я называю только vapply один раз вместо sapply/vapply комбо я имел раньше (я избавился от USE.NAMES=FALSE, поскольку он не имеет никакого значения). Хотя это немного чище, это не сильно меняло сроки на моей машине (я перепробовал тесты docendo с графиками, и это похоже почти одинаково).

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

Вот некоторые тайминги:

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 

Единственное, что немного озадачивает меня в том, что объект производства Triangle1 составляет половину размера всех других решений.

object.size(t1) 
400000200 bytes 

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

Когда я делаю некоторые проверки, это становится еще более запутанным.

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 

---

Как @Frank отметил в комментариях, t1 целая матрица, тогда как другие являются числовыми. Я должен был знать это, поскольку один из most important R functions сказал бы мне эту информацию с самого начала.

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 ...