Выберите подмножество столбцов, которые минимизируют критерий в R

У меня есть разреженный двоичный файл data.frame, который выглядит так

set.seed(123)
dat <- as.data.frame(matrix(rep(round(runif(40,0,0.9),0),5),ncol = 20))

#  > dat
#    V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V20
# 1   0  0  1  1  0  0  1  1  0   0   1   1   0   0   1   1   0   0   1   1
# 2   0  0  0  1  0  0  0  1  0   0   0   1   0   0   0   1   0   0   0   1
# 3   0  1  0  1  0  1  0  1  0   1   0   1   0   1   0   1   0   1   0   1
# 4   0  0  0  0  0  0  0  0  0   0   0   0   0   0   0   0   0   0   0   0
# 5   0  1  1  0  0  1  1  0  0   1   1   0   0   1   1   0   0   1   1   0
# 6   0  0  0  0  0  0  0  0  0   0   0   0   0   0   0   0   0   0   0   0
# 7   0  0  1  0  0  0  1  0  0   0   1   0   0   0   1   0   0   0   1   0
# 8   0  1  1  1  0  1  1  1  0   1   1   1   0   1   1   1   0   1   1   1
# 9   0  1  1  0  0  1  1  0  0   1   1   0   0   1   1   0   0   1   1   0
# 10  1  0  0  0  1  0  0  0  1   0   0   0   1   0   0   0   1   0   0   0

Мне нужно найти 3 столбца, которые минимизируют количество нулей, полученных при вызове rowSums для этих столбцов.

Пример:

 # > rowSums(dat[,1:3])
 # [1] 2 2 2 3 2 2 0 2 0 1
 # 
 # > rowSums(dat[,2:4])
 # [1] 3 2 3 3 1 2 1 1 0 1

Здесь, когда я вызываю rowSums для первых 3 столбцов, я получаю 2 нуля, а когда я вызываю rowSums для столбцов 2:4, я получаю только один 0, поэтому второе решение было бы предпочтительнее.

Конечно, мне не нужно, чтобы столбцы располагались рядом друг с другом, когда я применяю rowSums, поэтому мне нужно изучить все возможные комбинации (например, я хочу, чтобы rowSums учитывал также случай ov V1+V5+V17, ...), и если есть несколько «оптимальных» решений, я могу просто оставить одно из них.

Обратите внимание, что мои реальные data.frame составляют 220 000 строк x 200 столбцов, поэтому мне нужен эффективный подход с точки зрения потребляемого времени/памяти.


person hellter    schedule 08.07.2016    source источник


Ответы (1)


Это наиболее очевидное решение, хотя, вероятно, оно не очень хорошо масштабируется:

which.min(combn(dat,3L,function(x) sum(rowSums(x)==0)));
## [1] 2

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

cis <- combn(seq_along(dat),3L)[,2L];
cis;
## [1] 1 2 4

И тогда получить имена столбцов легко:

names(dat)[cis];
## [1] "V1" "V2" "V4"

Получить количество нулей в решении можно следующим образом:

sum(rowSums(dat[,cis])==0);
## [1] 1

Я написал гораздо более быстрое решение в Rcpp.

Чтобы сделать функцию более универсальной, я написал ее так, чтобы она брала логическую матрицу, а не data.frame, с целью поиска комбинации столбцов с наименьшим количеством полностью истинных строк. Таким образом, для вашего случая вы можете вычислить аргумент как dat==0. Я также параметризовал количество столбцов в комбинации как второй параметр r, который будет равен 3 для вашего случая.

library(Rcpp);
Sys.setenv('PKG_CXXFLAGS'='-std=c++11');

cppFunction('
    IntegerVector findColumnComboWithMinimumAllTrue(LogicalMatrix M,int r) {
        std::vector<int> rzFull(M.nrow()); std::iota(rzFull.begin(),rzFull.end(),0);
        std::vector<int> rzErase;
        std::vector<std::vector<int>> rzs(M.ncol(),std::vector<int>(M.nrow()));
        std::vector<std::vector<int>*> rzps(M.ncol());
        std::vector<int>* rzp = &rzFull;
        std::vector<int> com(r);
        int bestAllTrueCount = M.nrow()+1;
        std::vector<int> bestCom(r);
        int pmax0 = M.ncol()-r;
        int p = 0;
        while (true) {
            rzErase.clear();
            for (int rzi = 0; rzi < rzp->size(); ++rzi)
                if (!M((*rzp)[rzi],com[p])) rzErase.push_back(rzi);
            if (p+1==r) {
                if (rzp->size()-rzErase.size() < bestAllTrueCount) {
                    bestAllTrueCount = rzp->size()-rzErase.size();
                    bestCom = com;
                }
                if (com[p]==pmax0+p) {
                    do {
                        --p;
                    } while (p >= 0 && com[p]==pmax0+p);
                    if (p==-1) break;
                    ++com[p];
                    rzp = p==0 ? &rzFull : rzps[p-1];
                } else {
                    ++com[p];
                }
            } else {
                if (rzErase.empty()) {
                    rzps[p] = rzp;
                } else {
                    rzs[p].clear();
                    int rzi = -1;
                    for (int ei = 0; ei < rzErase.size(); ++ei)
                        for (++rzi; rzi < rzErase[ei]; ++rzi)
                            rzs[p].push_back((*rzp)[rzi]);
                    for (++rzi; rzi < rzp->size(); ++rzi)
                        rzs[p].push_back((*rzp)[rzi]);
                    rzp = rzps[p] = &rzs[p];
                }
                ++p;
                com[p] = com[p-1]+1;
            }
        }
        IntegerVector res(bestCom.size());
        for (int i = 0; i < res.size(); ++i)
            res[i] = bestCom[i]+1;
        return res;
    }
');

Вот демонстрация вашего примера ввода:

set.seed(123L);
dat <- as.data.frame(matrix(rep(round(runif(40,0,0.9),0),5),ncol=20L));
findColumnComboWithMinimumAllTrue(dat==0,3L);
## [1] 1 2 4

А вот полноразмерный тест, который на моей системе занимает почти 10 минут:

set.seed(1L); NR <- 220e3L; NC <- 200L;
dat <- as.data.frame(matrix(sample(0:1,NR*NC,T),NR,NC));
system.time({ findColumnComboWithMinimumAllTrue(dat==0,3L); });
##    user  system elapsed
## 555.641   0.328 556.401
res;
## [1] 28 64 89
person bgoldst    schedule 08.07.2016
comment
Большое спасибо за ваш ответ, вы были очень полезны, и я бы никогда не пришел к решению, подобному вашему. Ваша функция работает с r=3, но, к сожалению, она слишком медленная с r=5, который мне нужен. Я не написал это в вопросе, потому что не думал, что это будет критично, но на самом деле это так, так как с r=3 у нас есть ~1,3 миллиона возможных комбинаций, а с r=5 это число увеличивается до ~2,5 млрд (почти в 2000 раз больше). Извините за ошибку. Если вы видите способ улучшить функцию, это было бы здорово. В остальном все равно спасибо! - person hellter; 08.07.2016
comment
@hellter Всегда пожалуйста. Из любопытства, смогли ли вы найти какое-либо решение, которое могло бы справиться с делом r=5 за достаточно короткий период времени? - person bgoldst; 08.07.2016
comment
Я обдумываю это, но я еще не пришел к решению, и я не вижу простого способа сделать это. - person hellter; 08.07.2016