使用data.table进行矩阵运算和分量加法

Question

如果事先不知道要求和的矩阵数,那么进行分量矩阵加法的最佳方法是什么？更一般地说,有没有一种很好的方法在data.table的上下文中执行矩阵(或多维数组)操作？我data.table通过几个固定变量或类别对数据进行排序和分组的效率,每个变量或类别包含不同数量的观察值.

例如:

1. 找到数据的每个观察(行)中给出的向量分量的外积,返回每行的矩阵.
2. 在每组数据类别的所有行上按组件顺序对结果矩阵求和.

这里用2x2矩阵说明,只有一个类别:

library(data.table)

# example data, number of rows differs by category t
N <- 5
dt <- data.table(t = rep(c("a", "b"), each = 3, len = N), 
                 x1 = rep(1:2, len = N), x2 = rep(3:5, len = N),
                 y1 = rep(1:3, len = N), y2 = rep(2:5, len = N))
setkey(dt, t)
> dt
   t x1 x2 y1 y2
1: a  1  3  1  2
2: a  2  4  2  3
3: a  1  5  3  4
4: b  2  3  1  5
5: b  1  4  2  2

我尝试了一个函数来计算外积的矩阵和, %o%

mat_sum <- function(x1, x2, y1, y2){
  x <- c(x1, x2) # x vector
  y <- c(y1, y2) # y vector
  xy <- x %o% y # outer product (i.e. 2x2 matrix)
  sum(xy)  # <<< THIS RETURNS A SINGLE VALUE, NOT WHAT I WANT.
  }

当然,这不起作用,因为sum在数组中添加了所有元素.

我看到这个答案使用,Reduce('+', .list)但似乎需要已经list添加了所有矩阵.我还没弄清楚如何在内部做到这一点data.table,所以相反,我有一个繁琐的解决方法:

# extract each outer product component first...
mat_comps <- function(x1, x2, y1, y2){
  x <- c(x1, x2) # x vector
  y <- c(y1, y2) # y vector
  xy <- x %o% y # outer product (i.e. 2x2 matrix)
  xy11 <- xy[1,1]
  xy21 <- xy[2,1]
  xy12 <- xy[1,2]
  xy22 <- xy[2,2]
  return(c(xy11, xy21, xy12, xy22))
}

# ...then running this function on dt, 
# taking extra step (making column 'n') to apply it row-by-row...
dt[, n := 1:nrow(dt)]
dt[, c("xy11", "xy21", "xy12", "xy22") := as.list(mat_comps(x1, x2, y1, y2)), 
   by = n]

# ...then sum them individually, now grouping by t
s <- dt[, list(s11 = sum(xy11),
               s21 = sum(xy21),
               s12 = sum(xy12),
               s22 = sum(xy22)),
        by = key(dt)]
> s
   t s11 s21 s12 s22
1: a   8  26  12  38
2: b   4  11  12  23

并给出了总和的组件,最终可以转换回矩阵.

Answer 1

通常,data.table设计用于列.您将问题转化为整体操作的次数越多,您就越能摆脱困境data.table.

这是尝试完成此操作.可能有更好的方法.这更像是一个模板,提供了解决问题的想法(尽管我知道在所有情况下可能都不可能).

xcols <- grep("^x", names(dt))
ycols <- grep("^y", names(dt))
combs <- CJ(ycols, xcols)
len <- seq_len(nrow(combs))
cols = paste("V", len, sep="")
for (i in len) {
    c1 = combs$V2[i]
    c2 = combs$V1[i]
    set(dt, i=NULL, j=cols[i], value = dt[[c1]] * dt[[c2]])
}

#    t x1 x2 y1 y2 V1 V2 V3 V4
# 1: a  1  3  1  2  1  3  2  6
# 2: a  2  4  2  3  4  8  6 12
# 3: a  1  5  3  4  3 15  4 20
# 4: b  2  3  1  5  2  3 10 15
# 5: b  1  4  2  2  2  8  2  8

这基本上适用于外部产品.现在只需聚合它.

dt[, lapply(.SD, sum), by=t, .SDcols=cols]

#    t V1 V2 V3 V4
# 1: a  8 26 12 38
# 2: b  4 11 12 23

HTH

编辑:修改cols, c1, c2了一个位以获得具有正确顺序的输出V2和V3.