R中用于上三角矩阵的外部函数

Question

我目前有这段代码。

s1=seq(0,10,length.out=3)
s2=seq(0,10,length.out=3)
d=outer(s1,s2,`-`)
I=outer(s1,s2,`==`)

然而我只需要 d 的上三角部分，所以目前我正在执行以下操作

d=d[upper.tri(d,diag=T)]
I=I[upper.tri(I,diag=T)]

有没有一种方法可以通过将上三角矩阵纳入外部函数来使此代码更快？ 请注意，这是一个可重现的示例，但我的代码要复杂得多，并且我需要尽可能地减少运行时间。

- raK1

---

a 和 b 是什么？ - akrun

@akrun，s1和s2是相同的，而且s=b=s1=s2。 - raK1

这些向量长度相同吗（在您的原始数据中）？也许使用combn(seq_along(s1), 2, FUN= function(i) s1[i[1]]- s2[i[2]])。 - akrun

@akrun 谢谢，您的方法无法考虑对角线元素，有没有一种方法可以将其纳入combn函数中？ - raK1

如果顺序不重要，则 c(rbind(rep(0,length(s1)),combn(seq_along(s1), 2, FUN= function(i) s1[i[1]]- s2[i[2]])))。 - akrun

Answer 1

outer函数在内部使用rep函数来拉伸其参数:

    Y <- rep(Y, rep.int(length(X), length(Y)))
    if (length(X)) 
        X <- rep(X, times = ceiling(length(Y)/length(X)))

您可以使用子集操作完成同样的任务，但是只需生成三角形索引。在这种情况下，sequence函数非常有用：

> i <- sequence(1:3)
> j <- rep(1:3, 1:3)
> d = s1[i] - s2[j]
> I = s1[i] == s2[j]
> d
[1]   0  -5   0 -10  -5   0
> I
[1]  TRUE FALSE  TRUE FALSE FALSE  TRUE

- Neal Fultz

---

你好@Neal，如果i和j的长度很大，我认为这比outer慢。如果是这种情况，你有什么建议吗？ - Shanks

我通常建议使用Rcpp以获得最大的性能。如果使用outer函数，它需要进行两倍的计算量，相比直接计算上三角矩阵，我会感到惊讶它是否更快，但我认为这是可能的。 - Neal Fultz