R矩阵运算

Question

我有一个矩阵 (15000 x 3000)。目标是根据原始矩阵和初始值生成一个新的矩阵。例如，我想执行的标准是：

目前我的代码是这样设置的。

DF1[1,]=1

for( i in 2:2000 ) {
    for( j in 1:15000 ) {

              if(DF[j,i] == 1 && DF1[j-1,i] == 0)
                DF1[j,i] = 1
              else if(DF[j,i] == 0 && DF1[j-1,i] == 1)
                DF1[j,i] = 0
              else DF[j,i,1] = DF1[j-1,i]

    }
}

DF为原矩阵

DF1是新形成的矩阵

我的问题：还有其他方法吗？更快的方法？

因为嵌套循环不好用，我尝试用apply，但是不知道函数怎么写，因为涉及到两个矩阵

一个例子

x <- structure(c(1L, 0L, 0L, NA, NA, 0L, NA, 0L, 1L, 0L, 1L, 0L, 0L, 
NA, 0L, NA, 1L, NA, 1L, 0L, 1L, 0L, 1L, 0L), .Dim = c(4L, 6L), .Dimnames = list(
    NULL, NULL))
x
#     [,1] [,2] [,3] [,4] [,5] [,6]
#[1,]    1   NA    1    0    1    1
#[2,]    0    0    0   NA   NA    0
#[3,]    0   NA    1    0    1    1
#[4,]   NA    0    0   NA    0    0

循环（不起作用）

for( i in 1:4 ) { 
     for( j in 2:4 ) { 
         if(x[j,i] == 1 && y[j-1,i] == 0) { 
            y[j,i] = 1 
         }else{
            if(x[j,i] == 0 && y[j-1,i] == 1) {
               y[j,i] = 0 
         }else{ 
            y[j,i] = y[j-1,i]
         }
     }
   }

Answer 1

函数f1使用嵌套循环。（为了摆脱与 NA 比较导致非逻辑值 NA 的问题，我将 NA 替换为 Inf。）仔细阅读由循环表示的算法会得出替代方案 f2:

f1 <- function( x, initialValues = 1 )
{
  x[which(is.na(x))] <- Inf
  y <- matrix(NA,nrow(x),ncol(x))
  y[1,] <- initialValues

  for( i in 1:ncol(x) ) { 
    for( j in 2:nrow(x) ) { 
      if(x[j,i] == 1 && y[j-1,i] == 0) { 
        y[j,i] = 1 
      }else{
        if(x[j,i] == 0 && y[j-1,i] == 1) {
          y[j,i] = 0 
        }else{ 
          y[j,i] = y[j-1,i]
        }
      }
    }
  }
  return(y)
}

f2 <- function( x, initialValues = 1 )
{  
  g <- function(v)
  {
    m <- cumsum(!is.na(v))
    v[which(!is.na(v))[m]]
  }

  x[which(!(x %in% 0:1))] <- NA
  x[1,] <- initialValues
  return( apply(x,2,g) )
}

函数 g 填充向量 v 中的 NA-gaps：g(v)[i] 等于 v[j] 其中 j 是最大索引 j<=i 和 v[j]!=NA。（归纳证明：v[which(!is.na(v))] 包含 v 中的非 NA 值。如果 v[i]==NA 则 m[i]==m[i-1] 和 g(v)[i]==v[which(!is.na(v))[m[i]]]==v[which(!is.na(v))[m[i-1]]==g(v)[i-1]。否则 m[i]==m[i-1]+1，因此 g(v)[i-1]==v[which(!is.na(v))[m[i-1]]]==v[which(!is.na(v))][m[i-1]] 和 g(v)[i]==v[which(!is.na(v))[m[i]]]==v[which(!is.na(v))][m[i]]==v[which(!is.na(v))][m[i-1]+1]，下一个非 NA 值。）

f2 比 f1 快，特别是对于大矩阵。问题的小矩阵：

> library(microbenchmark)

> x <- structure(c(1L, 0L, 0L, NA, NA, 0L, NA, 0L, 1L, 0L, 1L, 0L, 0L, 
+                  NA, 0L, NA, 1L, NA, 1L, 0L, 1L, 0L, 1L, 0L), .Dim = c(4L, 6 .... [TRUNCATED] 

> microbenchmark( f1(x), f2(x) )
Unit: microseconds
  expr     min       lq     mean   median      uq     max neval
 f1(x) 433.864 461.2645 482.9120 471.6805 480.059 920.716   100
 f2(x) 379.518 387.6700 402.9235 391.7465 414.617 620.453   100

> all(f1(x)==f2(x))
[1] TRUE

更大的矩阵：

> set.seed(1)

> n <- 200

> m <- 300

> big_x <- matrix(sample(0:10,n*m,replace=TRUE),n,m)

> big_x[sample(1:(n*m),floor((n*m)/3))] <- NA

> microbenchmark( f1(big_x), f2(big_x) )
Unit: milliseconds
      expr       min        lq      mean    median        uq      max neval
 f1(big_x) 360.42174 495.63713 662.54576 772.42981 778.18100 890.0092   100
 f2(big_x)  33.54202  38.65849  62.25661  67.82429  72.42288 188.2729   100

> all(f1(big_x)==f2(big_x))
[1] TRUE
>

更大：

> set.seed(1)

> n <- 800

> m <- 1000

> huge_x <- matrix(sample(0:10,n*m,replace=TRUE),n,m)

> huge_x[sample(1:(n*m),floor((n*m)/3))] <- NA

> microbenchmark( f1(huge_x), f2(huge_x) )
Unit: milliseconds
       expr       min        lq     mean    median       uq       max neval
 f1(huge_x) 4002.4121 7759.2438 8149.821 8466.4698 8950.172 10087.251   100
 f2(huge_x)  311.4259  520.5374  751.874  774.2699 1010.188  1228.504   100

> all(f1(huge_x)==f2(huge_x))
[1] TRUE
>

一个大小为15000乘以3000的矩阵，问题中提到：

> set.seed(1)

> n <- 15000

> m <- 3000

> x_15k.3k <- matrix(sample(0:1,n*m,replace=TRUE),n,m)

> x_15k.3k[sample(1:(n*m),floor((n*m)/3))] <- NA

> microbenchmark( f1(x_15k.3k), f2(x_15k.3k), times=1 )
Unit: seconds
         expr       min        lq      mean    median        uq       max
 f1(x_15k.3k) 389.47262 389.47262 389.47262 389.47262 389.47262 389.47262
 f2(x_15k.3k)  19.97606  19.97606  19.97606  19.97606  19.97606  19.97606
 neval
     1
     1

> all(f1(x_15k.3k)==f2(x_15k.3k))
[1] TRUE
>

R矩阵运算

R matrix operation

r

matrix-multiplication