在特定条件下创建二元矩阵答案

【问题标题】：Create a binary matrix under certain conditions在特定条件下创建二元矩阵
【发布时间】：2017-03-27 14:20:54
【问题描述】：

我正在尝试创建一个给定 m 和 p 的函数，它返回一个包含 m 行和 mxp 列的矩阵。矩阵应该有0，除了p 位置，从p（行数）开始。

例如，给定m=4 和p=2，矩阵应如下所示：

1    1    0    0    0    0    0    0
0    0    1    1    0    0    0    0
0    0    0    0    1    1    0    0
0    0    0    0    0    0    1    1

我想处理大矩阵。我知道如何使用其他编程语言（如 python）中的循环来执行此操作，但我确信它应该是在 R 中执行此操作的一种更简单、更优雅的方式。我在diag() 上玩了一段时间，但没有找到解决方案。

【问题讨论】：

标签： r matrix

【解决方案1】：

apply()ing rep() 函数到对角矩阵的每一行（或列，都是一样的）：

t(apply(diag(m), 2, rep, each = p))

#      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
# [1,]    1    1    0    0    0    0    0    0
# [2,]    0    0    1    1    0    0    0    0
# [3,]    0    0    0    0    1    1    0    0
# [4,]    0    0    0    0    0    0    1    1

【讨论】：

【解决方案2】：

p=2的这个解决方案使用了行数的变化：

m <- 4
d <- diag(m)
matrix(rbind(d,d), m)
#      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
# [1,]    1    1    0    0    0    0    0    0
# [2,]    0    0    1    1    0    0    0    0
# [3,]    0    0    0    0    1    1    0    0
# [4,]    0    0    0    0    0    0    1    1

对于p的其他值（来自A5C1D2H2I1M1N2O1R2T1的评论）：

p <- 3; m <- 4
matrix(rep(diag(m), each = p), nrow = m, byrow = TRUE)

【讨论】：

【解决方案3】：

这个怎么样：

f <- function(m, p){
     a <- diag(m)
     a[,rep(seq_len(m), each=p)]
}

> f(m = 4, p = 2)

#     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#[1,]    1    1    0    0    0    0    0    0
#[2,]    0    0    1    1    0    0    0    0
#[3,]    0    0    0    0    1    1    0    0
#[4,]    0    0    0    0    0    0    1    1

> f(m = 3, p = 4)

#     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
#[1,]    1    1    1    1    0    0    0    0    0     0     0     0
#[2,]    0    0    0    0    1    1    1    1    0     0     0     0
#[3,]    0    0    0    0    0    0    0    0    1     1     1     1

我们的想法是首先创建一个大小为m（我们将其命名为a）的对角矩阵，然后将该矩阵的每一列重复p 次（即m*p 矩阵）。

【讨论】：

【解决方案4】：

此方法使用矩阵子集填充1。

myMatFunc <- function(m, p) {
  # initialize matrix of correct size, filled with 0s
  myMat <- matrix(0L, m, m * p)
  #fill in 1s using matrix subsetting
  myMat[cbind(rep(seq_len(m), each=p), seq_len(m * p))] <- 1L

  myMat
}

那么，

myMatFunc(4, 2)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    1    0    0    0    0    0    0
[2,]    0    0    1    1    0    0    0    0
[3,]    0    0    0    0    1    1    0    0
[4,]    0    0    0    0    0    0    1    1

感谢来自@joseph-wood、@jogo 和 @A5C1D2H2I1M1N2O1R2T1 的 cmets，我提高了删除对 nrow 的调用和对 ncol 的调用的效率，通过转换将矩阵的大小减半转换为整数，并修复了初始测试拼写错误。

【讨论】：

不错的答案...但是我认为应该是myMat <- matrix(0, m, m*p)
感谢您发现这个错字。
更短的代码：myMat[cbind(rep(1:m, each=p), 1:(m*p))] <- 1（但最终不会更快）
@jogo，更长的代码，但最快的选项在这里：myfun <- function(m, p) {M <- matrix(0L, ncol = m*p, nrow = m); C <- seq.int(m*p); R <- rep(seq.int(m), each = p);M[cbind(R, C)] <- 1L;M}。虽然如果效率真的是这样一个功能的关注点，我会感到惊讶......

【解决方案5】：

这是一个非常快的基本 R 解决方案：

Joseph <- function(m, p) {
  mat <- matrix(0L, nrow = m, ncol = m*p)
  for (i in 1:m) {mat[i, p*(i-1L) + 1:p] <- 1L}
  mat
}

以下是一些相等比较：

fun989 <- function(m, p){
  a <- diag(m)
  a[,rep(seq_len(m), each=p)]
}

IMO <- function(m, p) {
  myMat <- matrix(0L, m, m*p)
  myMat[cbind(rep(seq_len(nrow(myMat)), each=p), seq_len(ncol(myMat)))] <- 1
  myMat
}

JOGO <- function(m, p) {matrix(rep(diag(m), each = p), nrow = m, byrow = TRUE)}
APOM <- function(m, p) {t(apply(diag(m), 2, rep, each = p))}

library(compiler)
enableJIT(3)  ## compiling each function
all.equal(Joseph(100, 50), fun989(100, 50))
[1] TRUE
all.equal(Joseph(100, 50), APOM(100, 50))
[1] TRUE
all.equal(Joseph(100, 50), JOGO(100, 50))
[1] TRUE
all.equal(Joseph(100, 50), IMO(100, 50))
[1] TRUE
enableJIT(0)  ## return to standard setting

以下是基准：

library(microbenchmark)

microbenchmark(Joseph(100, 50), JOGO(100, 50), fun989(100, 50), APOM(100, 50), IMO(100, 50), unit = "relative")
Unit: relative
           expr       min        lq     mean    median        uq      max neval cld
Joseph(100, 50)  1.000000  1.000000 1.000000  1.000000  1.000000 1.000000   100  a 
  JOGO(100, 50) 33.388929 20.892988 6.593804 22.365625 19.161056 1.167957   100   b
fun989(100, 50)  7.192071  4.577225 2.044973  4.432824  4.129563 1.029050   100  a 
  APOM(100, 50) 40.244128 28.176729 8.805715 27.785985 23.966477 1.209582   100   b
   IMO(100, 50)  6.119685  3.898451 2.712222  6.192030  6.033916 1.044422   100  a

【讨论】：

【解决方案6】：

这是另一种方法，但我会选择@989 答案而不是我的答案；

 cadv.func = function(m,p)
{

  cmat <- matrix(data=NA,nrow=m,ncol=m*p)
  cmat[is.na(cmat)] <- 0

  for (i in 1:m){
    for (j in 1:p){

    cmat[i,j+p*(i-1)] = 1

  } 
  }

  return(cmat)
}

cadv.func(4,2)


 #       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
 # [1,]    1    1    0    0    0    0    0    0
 # [2,]    0    0    1    1    0    0    0    0
 # [3,]    0    0    0    0    1    1    0    0
 # [4,]    0    0    0    0    0    0    1    1

【讨论】：

这不会返回正确的结果。我认为问题在于cmat[i,j + 2*(i-j)] = 1
应该是cmat[i,j + p*(i-j)] = 1
@JosephWood 感谢您的评论。我犯了一个错误。 2 仅适用于p=2。但是(i-1) 应该保持不变。