Операция R-матрицы

Функция 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 заполняет NA пробелов в векторе v: g(v)[i] равно v[j], где j — наибольший индекс, такой что j<=i и v[j]!=NA. (Доказательство по индукции: v[which(!is.na(v))] содержит не-NA значения в v. Если 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
>