Наружное / тензорное произведение в R - PullRequest
Новая
       36

Наружное / тензорное произведение в R

7 голосов
/ 07 января 2012

Учитывая p векторов x1,x2,...,xp каждого измерения d, каков наилучший способ вычислить их тензорное / внешнее / произведение Крускала (массив p с записями X[i1,i2,..ip] = x1[i1]x2[i2]...xp[ip])? Циклирование тривиально, но Глупо. Использование повторных вызовов outer работает нормально, но не кажется оптимальным решением (и будет увеличиваться с ростом p, очевидно). Есть ли лучший способ?

Edit:

Мой лучший результат 

array(apply(expand.grid(x1, x2, x3), 1, prod), dim=rep(d, 3))

который, по крайней мере, чувствует себя лучше ...

Редактировать 2: В ответ на @Dwin вот завершено пример 

d=3
x1 = 1:d
x2 = 1:d+3
x3 = 1:d+6
array(apply(expand.grid(x1, x2, x3), 1, prod), dim=rep(d, 3))

, , 1

     [,1] [,2] [,3]
[1,]   28   35   42
[2,]   56   70   84
[3,]   84  105  126

, , 2

     [,1] [,2] [,3]
[1,]   32   40   48
[2,]   64   80   96
[3,]   96  120  144

, , 3

     [,1] [,2] [,3]
[1,]   36   45   54
[2,]   72   90  108
[3,]  108  135  162

Ответы [ 3 ]

7 голосов
/ 08 января 2012

Трудно превзойти производительность outer.В результате выполняется умножение матриц, которое выполняется библиотекой BLAS.Повторный вызов outer также не имеет значения, так как последний вызов будет доминировать как по скорости, так и по памяти.Например, для векторов длины 100 последний вызов по крайней мере в 100 раз медленнее, чем предыдущий ...

Лучшим вариантом для достижения максимальной производительности является получение лучшей библиотеки BLAS для R.по умолчанию один не очень хорошо.В Linux вы можете довольно легко настроить R для использования ATLAS BLAS.На Windows это сложнее, но возможно.См. R для Windows FAQ . 

# multiple outer
mouter <- function(x1, ...) { 
    r <- x1
    for(vi in list(...)) r <- outer(r, vi)
    r
}

# Your example
d=3
x1 = 1:d
x2 = 1:d+3
x3 = 1:d+6 
mouter(x1,x2,x3)

# Performance test
x <- runif(1e2)
system.time(mouter(x,x,x))   # 0 secs (less than 10 ms)
system.time(mouter(x,x,x,x)) # 0.5 secs / 0.35 secs (better BLAS)

Я заменил свою Windows Rblas.dll на версию GOTO BLAS DYNAMIC_ARCH в это место , что улучшило время с 0,5до 0,35 с, как показано выше.

1 голос
/ 07 января 2012

Я задаюсь вопросом, подходит ли вам продукт kronecker.Я не могу сказать из вашего описания проблемы, что именно нужно, но элементы этого в небольшом наборе аргументов одинаковы (хотя и в другом расположении, чем те, что были получены в результате решения Чаласани, которое вы критиковали как медленное: 

kronecker( outer(LETTERS[1:2], c(3, 4, 5),FUN=paste), letters[6:8] ,FUN=paste)
     [,1]    [,2]    [,3]   
[1,] "A 3 f" "A 4 f" "A 5 f"
[2,] "A 3 g" "A 4 g" "A 5 g"
[3,] "A 3 h" "A 4 h" "A 5 h"
[4,] "B 3 f" "B 4 f" "B 5 f"
[5,] "B 3 g" "B 4 g" "B 5 g"
[6,] "B 3 h" "B 4 h" "B 5 h"

Если вы хотите продукты, то замените либо prod, либо "*". В любом случае, предлагая выборочный набор векторов и желаемый результат, лучше всего ставить вопросы.

1 голос
/ 07 января 2012

Вы можете использовать тензор пакет.

А также %o% функция 

A <- matrix(1:6, 2, 3)
D <- A %o% A
...