1

我试图在以下函数中分解以下代码：

paired <- function(x) crossprod(x[x]-1:length(x))==0

对于上下文，这出现在 CV 的这个答案中。

我们有一个由十个 8 维向量组成的矩阵，对应于向量的随机排列1:8

n <- 8
x <- replicate(10, sample(1:n, n))  

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    1    3    5    2    5    1    8    8     3
[2,]    4    3    5    4    3    8    5    2    1     8
[3,]    5    5    4    3    8    2    6    3    6     7
[4,]    6    8    1    7    4    3    8    6    2     2
[5,]    7    4    8    2    7    4    4    1    3     4
[6,]    2    6    2    6    5    1    3    4    7     6
[7,]    8    2    7    8    1    6    7    5    5     1
[8,]    3    7    6    1    6    7    2    7    4     5

我们按列应用该函数paired以apply(x, 2, paired)获取布尔运算的结果，我猜这会评估是否有任何结果function(x) crossprod(x[x]-1:length(x))等于零：

intermediate <- function(x) crossprod(x[x]-1:length(x))
apply(x, 2, intermediate)
 [1]  74  80  88  88 100 108  26  90 124 106

apply(x, 2, paired)
 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

调用中的值intermediate不是按列中的列向量x[x] - 1与它们自身的点积，对于第一列，它是：

> w <- c(1,4,5,6,7,2,8,3)
> v <- w[w] - 1
> t(v)%*%v
     [,1]
[1,]  140

而不是74.

那么问题来了length(x)，修改的作用是什么crossprod()？

编辑：按照评论（谢谢），如果我尝试以下代码：

w <- c(1,4,5,6,7,2,8,3)
v <- w[w] - 1:length(x)
t(v)%*%v

我crossprod(v) [1,] 146660也变得不同74。

4

1 回答 1

2

不要太挂crossprod在paired函数中。它所做的只是检查x[x] - 1:length(x)一个零向量（即“完美配对”的条件）。它可以以不同的方式更快地编码（参见paired2或paired3）：

> n <- 8
> set.seed(17)
> x <- replicate(1e6, sample(1:n, n))
>
> paired <- function(x) crossprod(x[x] - 1:length(x))==0
> paired2 <- function(x) sum(x[x]==1:length(x))==length(x)
> paired3 <- function(x) sum(abs(x[x]-1:length(x)))==0
>
> system.time(i.paired  <- apply(x, 2, paired))
   user  system elapsed
  9.812   0.000   9.821
> system.time(i.paired2 <- apply(x, 2, paired2))
   user  system elapsed
  4.548   0.000   4.550
> system.time(i.paired3 <- apply(x, 2, paired3))
   user  system elapsed
  4.617   0.000   4.617
>
> all.equal(i.paired,i.paired2)
[1] TRUE
> all.equal(i.paired,i.paired3)
[1] TRUE
于 2015-12-17T23:19:36.730 回答