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

Проблема довольно глупа, но мне интересно, не хватает ли я чего-то. Скажем, что существует вектор 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]
}

Ответы

Ответ 1

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

Ответ 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.

Результаты:

pic

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

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

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.

Ответ 3

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

Ответ 4

Вот действительно быстрое базовое решение R:

Обновление

Я немного изменил код, поэтому я вызывал только 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 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 going on here
identical(t1,t2)
[1] FALSE

identical(t2,t3)
[1] TRUE

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

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