Как найти и назвать непрерывные ненулевые элементы в разреженной матрице в R?

Моя проблема концептуально проста. Я ищу эффективное с вычислительной точки зрения решение (мое собственное я прикрепляю в конце).

Предположим, у нас есть потенциально очень большая разреженная матрица, подобная той, что слева внизу, и мы хотим «назвать» каждую область смежных ненулевых элементов отдельным кодом (см. матрицу справа)

1 1 1 . . . . .          1 1 1 . . . . .
1 1 1 . 1 1 . .          1 1 1 . 4 4 . .
1 1 1 . 1 1 . .          1 1 1 . 4 4 . .
. . . . 1 1 . .   --->   . . . . 4 4 . .
. . 1 1 . . 1 1          . . 3 3 . . 7 7
1 . 1 1 . . 1 1          2 . 3 3 . . 7 7
1 . . . 1 . . .          2 . . . 5 . . .
1 . . . . 1 1 1          2 . . . . 6 6 6

В моем приложении смежные элементы будут образовывать прямоугольники, линии или отдельные точки, и они могут касаться друг друга только вершинами (т.е. в матрице не будет неправильных/непрямоугольных областей).

Решение, которое я придумал, состоит в том, чтобы сопоставить индексы строк и столбцов разреженного матричного представления с вектором с соответствующими значениями (кодами «имени»). Мое решение использует несколько for loops и отлично работает для малых и средних матриц, но быстро застревает в циклах, когда размеры матрицы становятся большими (> 1000). Вероятно, это зависит от того, что я не настолько продвинут в программировании на R - я не смог найти какой-либо вычислительный трюк/функцию, чтобы решить ее лучше.

Может ли кто-нибудь предложить более эффективный в вычислительном отношении способ сделать это в R?

Мое решение:

mySolution <- function(X){

  if (class(X) != "ngCMatrix") {stop("Input must be a Sparse Matrix")}
  ind <- which(X == TRUE, arr.ind = TRUE)
  r <- ind[,1]
  c <- ind[,2]

  lr <- nrow(ind)
  for (i in 1:lr) {
    if(i == 1) {bk <- 1}
    else {
      if (r[i]-r[i-1] == 1){bk <- c(bk, bk[i-1])}
      else {bk <- c(bk, bk[i-1]+1)}
    }
  }

  for (LOOP in 1:(lr-1)) {
    tr <- r[LOOP]
    tc <- c[LOOP]
    for (j in (LOOP+1):lr){
      if (r[j] == tr) {
        if(c[j] == tc + 1) {bk[j] <- bk[LOOP]} 
      }
    }
  }

  val <- unique(bk)
  for (k in 1:lr){
    bk[k] <- which(val==bk[k])
  }

  return(sparseMatrix(i = r, j = c, x = bk))
}

Заранее спасибо за любую помощь или указатель.


person GiuGe    schedule 03.02.2017    source источник
comment
Может быть хорошим вариантом использования для Rcpp. Ваш код кажется достаточно простым для перевода   -  person Aurèle    schedule 03.02.2017
comment
@apom может ты и прав, но я сейчас просматриваю Rcpp виньетки, и если мне нужно знать C++ язык, как кажется, это не для меня, к сожалению...   -  person GiuGe    schedule 03.02.2017
comment
Предполагая, что m — это ваш ngCMatrix, вы можете попытаться найти способ заставить что-то вроде sm = summary(m); sparseMatrix(i = sm$i, j = sm$j, x = cutree(hclust(dist(sm, "maximum")), h = 2)) работать должным образом.   -  person alexis_laz    schedule 03.02.2017
comment
@alexis_laz, я впечатлен - так элегантно. Шапо, я бы никогда не нашел такого решения. Я запускаю несколько симуляций, но, похоже, это работает. Спасибо за просвещение. :)   -  person GiuGe    schedule 04.02.2017
comment
... ну, на самом деле я говорил слишком рано. Я все равно попытаюсь заставить его работать должным образом, но кластеризация областей, соединенных вершиной, или даже больше, чем расстояние = 2 от границы до границы, просто плохо себя ведет...   -  person GiuGe    schedule 05.02.2017
comment
@GiuGe: я знаю; не мог найти способ заставить его работать. В том же примечании, принимая во внимание условие только прямоугольников и вычисление пользовательского расстояния sm = as.matrix(summary(m)); d = as.dist(sapply(1:nrow(sm), function(i) rowSums(abs(sm[i, col(sm)] - sm)))), использование sparseMatrix(i = sm[, "i"], j = sm[, "j"], x = cutree(hclust(d, "single"), h = 1)), кажется, отлично работает для текущего примера, хотя я склонен полагать, что будет много случаев, которые его нарушат.   -  person alexis_laz    schedule 05.02.2017
comment
@alexis_laz Бинго! Это оно! Вы можете изменить вычисление d на d <- dist(sm, method = "manhattan"), потому что это в основном то, что вы там делаете. Агрегация с одной связью и разрез из трех на высоте 1 эффективно различают области, которые «соприкасаются» друг с другом в вершине (и, следовательно, находятся на расстоянии 2 шагов друг от друга): как я сказал в своем вопросе, это ситуация, которую я имею в мои данные, так что это нормально. Пара симуляций с моей функцией и вашим методом говорят, что ваш метод быстрее в 4 раза с матрицами 200x200. Почему бы вам не опубликовать свой ответ ниже? Спасибо   -  person GiuGe    schedule 06.02.2017


Ответы (1)


В значительной степени опираясь на тот факт, что все соседние элементы, которые необходимо сгруппировать, образуют только прямоугольники/линии/точки, мы видим, что элементы матрицы могут быть агрегированы на основе их [row, col] индексов в матрице по отношению (abs(row1 - row2) + abs(col1 - col2)) < 2.

Итак, начиная с индексов [row, col]:

sm = as.matrix(summary(m))

Мы вычисляем их расстояние, которое, как отмечает GiuGe, на самом деле является «манхэттенским» методом:

d = dist(sm, "manhattan")

Здесь полезно свойство одиночной связи группировать элементы по их ближайшим соседям. Кроме того, мы можем получить группировку элементов с помощью cutreeing на «h = 1» (где расстояние индексов равно «‹ 2»):

gr = cutree(hclust(d, "single"), h = 1)

Наконец, мы можем обернуть вышеуказанное в новую разреженную матрицу:

sparseMatrix(i = sm[, "i"], j = sm[, "j"], x = gr)
#8 x 8 sparse Matrix of class "dgCMatrix"
#                    
#[1,] 1 1 1 . . . . .
#[2,] 1 1 1 . 4 4 . .
#[3,] 1 1 1 . 4 4 . .
#[4,] . . . . 4 4 . .
#[5,] . . 3 3 . . 7 7
#[6,] 2 . 3 3 . . 7 7
#[7,] 2 . . . 5 . . .
#[8,] 2 . . . . 6 6 6

«м» используется:

library(Matrix)
m = new("ngCMatrix"
    , i = c(0L, 1L, 2L, 5L, 6L, 7L, 0L, 1L, 2L, 0L, 1L, 2L, 4L, 5L, 4L, 
5L, 1L, 2L, 3L, 6L, 1L, 2L, 3L, 7L, 4L, 5L, 7L, 4L, 5L, 7L)
    , p = c(0L, 6L, 9L, 14L, 16L, 20L, 24L, 27L, 30L)
    , Dim = c(8L, 8L)
    , Dimnames = list(NULL, NULL)
    , factors = list()
)

ИЗМЕНИТЬ 10 фев. 2017

Другая идея (опять же, принимая во внимание тот факт, что соседние элементы образуют только прямоугольники/линии/точки) состоит в том, чтобы перебирать - в восходящих столбцах - по [row, col] индексам и на каждом шаге находить расстояние каждого элемента от его ближайшего соседа в текущий столбец и строка. Если найдено расстояние «‹ 2», то элемент группируется со своим соседом, в противном случае начинается новая группа. Завернутый в функцию:

ff = function(x) 
{
    sm = as.matrix(summary(x))

    gr = integer(nrow(sm)); ngr = 0L ; gr[1] = ngr 

    lastSeenRow = integer(nrow(x))
    lastSeenCol = integer(ncol(x))

    for(k in 1:nrow(sm)) {
        kr = sm[k, 1]; kc = sm[k, 2]
        i = lastSeenRow[kr]
        j = lastSeenCol[kc]

        if(i && (abs(kc - sm[i, 2]) == 1)) gr[k] = gr[i]
        else if(j && (abs(kr - sm[j, 1]) == 1)) gr[k] = gr[j]  
             else { ngr = ngr + 1L; gr[k] = ngr } 

        lastSeenRow[kr] = k
        lastSeenCol[kc] = k        
    }

    sparseMatrix(i = sm[, "i"], j = sm[, "j"], x = gr)                 
}                  

И нанес на "м":

ff(m)
#8 x 8 sparse Matrix of class "dgCMatrix"
#                    
#[1,] 1 1 1 . . . . .
#[2,] 1 1 1 . 4 4 . .
#[3,] 1 1 1 . 4 4 . .
#[4,] . . . . 4 4 . .
#[5,] . . 3 3 . . 7 7
#[6,] 2 . 3 3 . . 7 7
#[7,] 2 . . . 5 . . .
#[8,] 2 . . . . 6 6 6

Также удобно, что обе функции возвращают группы в одном и том же порядке, в чем мы можем убедиться:

identical(mySolution(m), ff(m))
#[1] TRUE

На, казалось бы, более сложном примере:

mm = new("ngCMatrix"
    , i = c(25L, 26L, 27L, 25L, 29L, 25L, 25L, 17L, 18L, 26L, 3L, 4L, 5L, 
14L, 17L, 18L, 25L, 27L, 3L, 4L, 5L, 17L, 18L, 23L, 26L, 3L, 
4L, 5L, 10L, 17L, 18L, 9L, 11L, 17L, 18L, 10L, 17L, 18L, 3L, 
17L, 18L, 21L, 17L, 18L, 17L, 18L, 1L, 2L, 3L, 4L, 16L, 8L, 17L, 
18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 7L, 9L, 10L, 11L, 26L, 
8L, 27L, 1L, 2L, 28L, 1L, 2L, 15L, 27L, 1L, 2L, 21L, 22L, 1L, 
2L, 7L, 21L, 22L, 1L, 2L, 6L, 24L, 1L, 2L, 5L, 11L, 16L, 25L, 
26L, 27L, 4L, 15L, 17L, 19L, 25L, 26L, 27L, 3L, 16L, 25L, 26L, 
27L, 2L, 28L, 1L)
    , p = c(0L, 0L, 3L, 3L, 5L, 6L, 7L, 7L, 10L, 18L, 25L, 31L, 35L, 38L, 
42L, 44L, 46L, 51L, 61L, 66L, 68L, 71L, 75L, 79L, 84L, 88L, 96L, 
103L, 108L, 110L, 111L)
    , Dim = c(30L, 30L)
    , Dimnames = list(NULL, NULL)
    , factors = list()
)
identical(mySolution(mm), ff(mm))
#[1] TRUE

И простой бенчмарк на большей матрице:

times = 30 # times `dim(mm)`
MM2 = do.call(cbind, rep_len(list(do.call(rbind, rep_len(list(mm), times))), times))
dim(MM2)
#[1] 900 900

system.time({ ans1 = mySolution(MM2) })
#   user  system elapsed 
# 449.50    0.53  463.26

system.time({ ans2 = ff(MM2) })
#   user  system elapsed 
#   0.51    0.00    0.52

identical(ans1, ans2)
#[1] TRUE
person alexis_laz    schedule 08.02.2017
comment
Спасибо @alexis_laz. Я пробовал этот код в последние дни. Результаты правильные, и это быстрее, чем мое решение. Тем не менее, с очень большими матрицами (> 1000x1000) он все еще может забивать оперативную память при построении матрицы расстояний. Так что новые подходы тоже приветствуются. - person GiuGe; 10.02.2017
comment
@GiuGe: Да, хотя dist сохраняет не всю матрицу расстояний, а ее lower.tri, она все еще требует значительной памяти. Я добавил еще один подход, который, протестировав его на некоторых случаях, кажется правильным; надеюсь, я ничего не пропустил в процессе. - person alexis_laz; 10.02.2017
comment
Фантастика! @alexis_laz ты гений! Большое спасибо. - person GiuGe; 10.02.2017