Filtering a logical and symmetric matrix

Viewed 147

I have a logical and symmetric matrix and I want to create blocks depending on the value I have on the diagonal. Here is my matrix:

MAT <- data.frame(matrix(data = c(rep(TRUE, 9), rep(FALSE, 3), rep(TRUE, 3), FALSE, rep(TRUE, 3), rep(FALSE, 2), TRUE, 
                FALSE, rep(TRUE, 2), FALSE, TRUE, FALSE, rep(TRUE, 3), FALSE, rep(TRUE, 4)), 6))

 [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
[1,] TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[2,] TRUE  TRUE  TRUE FALSE FALSE FALSE
[3,] TRUE  TRUE  TRUE FALSE  TRUE  TRUE
[4,] TRUE FALSE FALSE  TRUE FALSE  TRUE
[5,] TRUE FALSE  TRUE FALSE  TRUE  TRUE
[6,] TRUE FALSE  TRUE  TRUE  TRUE  TRUE

Here is the desired matrix: enter image description here We can see that a block is based on TRUE values. All columns and rows must be related to each other by TRUE in order to be a block. Once we have even one FALSE there is no block. We can have a block constituted of only one column as in the fourth column.

My goal is to put the same number in each column belonging to one block as seen in the image.

2 Answers

The question states that MAT is a matrix but, in fact, it is defined there as a data.frame. This is important because all(MAT) where MAT is as defined in the question gives an error in R 4.0 but not in R 4.1 so be sure you are using R 4.1 with the code below. Alternately convert MAT to a matrix using MAT <- as.matrix(MAT) in which case the code below works with both 4.1 and 4.0 and older.

1) Loop is.complete(i, j) is TRUE if all cells in MAT[i:j, i:j] are TRUE. If p to i is complete and p to i+1 is not complete then p to i is a block so record i in d and then at the end find the length of each block and use rep to create the desired vector placing it in the column names since one cannot mix logical and integer values together.

is.complete <- function(i, j) all(MAT[i:j, i:j])

i <- p <- 1
d <- c()
for(i in 1:nrow(MAT)) {
  ok <- is.complete(p, i) && (i == nrow(MAT) || !is.complete(p, i+1) )
  if (ok) { p <- i+1; d <- c(d, i) }
}
colnames(MAT) <- rep(seq_along(d), diff(c(0, d)))
MAT

giving:

     1     1     1     2     3     3
1 TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
2 TRUE  TRUE  TRUE FALSE FALSE FALSE
3 TRUE  TRUE  TRUE FALSE  TRUE  TRUE
4 TRUE FALSE FALSE  TRUE FALSE  TRUE
5 TRUE FALSE  TRUE FALSE  TRUE  TRUE
6 TRUE FALSE  TRUE  TRUE  TRUE  TRUE

2) Reduce Another possibility is Reduce. We accumulate the current start of block, p. At each step we consider the next index i and if the block from p to i is entirely TRUE then we continue p for consideration onwards; otherwise we start a new block at i. In the end Reduce returns a vector of the start of block repeated for each row (or column) of the block. We can then convert to factor and take the integer levels to get 1, 2, ... This approach is compact and uses no explicit loops.

f <- function(p, i) if (all(MAT[p:i, p:i])) p else i
colnames(MAT) <- as.integer(factor(Reduce(f, 1:nrow(MAT), acc = TRUE)))
MAT

giving:

     1     1     1     2     3     3
1 TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
2 TRUE  TRUE  TRUE FALSE FALSE FALSE
3 TRUE  TRUE  TRUE FALSE  TRUE  TRUE
4 TRUE FALSE FALSE  TRUE FALSE  TRUE
5 TRUE FALSE  TRUE FALSE  TRUE  TRUE
6 TRUE FALSE  TRUE  TRUE  TRUE  TRUE

3) maximize cells in blocks The solutions so far are greedy in nature. They take the largest block and when no larger one can be found start a new block; however, if it were desired to maximize the number of cells in blocks then that might not give the maximum value.

If we regard MAT as an adjacency matrix we can convert it to a graph in which case blocks correspond to complete subgraphs, also known as cliques. We generate all cliques and form an nrow(MAT) row matrix K with one clique per column such that each column is a 0/1 vector indicating which rows in MAT are in that clique. If there are p columns in that matrix then we form an objective p-vector v such that v[i] equals sum(K[, i])^2 which is the number of cells in the block that corresponds to clique i. From this we form the integer linear program to maximize the number of cells in blocks by choosing the 0/1 vector which satisfies:

max v'x such that Kx = 1 
x is 0/1 vector

where 1 is a vector of nrow(K) ones. This is known as a set partitioning problem.

library(igraph)
cli <- cliques(graph_from_adjacency_matrix(as.matrix(MAT), mode = "undirected"))
K <- +sapply(cli, function(x) colnames(MAT) %in% names(x))

library(lpSolve)
obj <- colSums(K)^2
res <- lp("max", obj, K, "=", 1, all.bin = TRUE)
Ksoln <- K[, res$solution == 1]
o <- order(apply(Ksoln, 2, which.max))
colnames(MAT) <- Ksoln[, o] %*% 1:ncol(Ksoln)
MAT

giving:

     1     2     1     3     1     1
1 TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
2 TRUE  TRUE  TRUE FALSE FALSE FALSE
3 TRUE  TRUE  TRUE FALSE  TRUE  TRUE
4 TRUE FALSE FALSE  TRUE FALSE  TRUE
5 TRUE FALSE  TRUE FALSE  TRUE  TRUE
6 TRUE FALSE  TRUE  TRUE  TRUE  TRUE

or sorting by block:

o <- order(as.numeric(colnames(MAT)))
MAT2 <- as.matrix(MAT)[o, o]
colnames(MAT2) <- sub("\\.*", "", colnames(MAT2))
MAT2

giving this reordered matrix:

        1     1     1     1     2     3
[1,] TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[2,] TRUE  TRUE  TRUE  TRUE  TRUE FALSE
[3,] TRUE  TRUE  TRUE  TRUE FALSE FALSE
[4,] TRUE  TRUE  TRUE  TRUE FALSE  TRUE
[5,] TRUE  TRUE FALSE FALSE  TRUE FALSE
[6,] TRUE FALSE FALSE  TRUE FALSE  TRUE

We see that with this solution there are 18 cells in blocks

res$objval
## [1] 18

whereas for the greedy solutions there are only 3^2 + 1^2 + 2^2 = 14 cells in blocks.

4) In the comments the poster indicated they are interested in the variation of (3) in which the groups are constrained to be consecutive. To do that simply restrict the K matrix to those columns that represent consecutive groups, i.e. add the two lines of code involving is.consec marked with ## below. We also use the modified input that the poster has defined in the comments below and have also marked those lines with ##. The code is otherwise identical to (3).

MAT <- data.frame(matrix(data = c(rep(TRUE, 9), rep(FALSE, 3), 
  rep(TRUE, 3), FALSE, rep(TRUE, 3), rep(FALSE, 2), TRUE, 
  FALSE, rep(TRUE, 2), FALSE, TRUE, FALSE, rep(TRUE, 3), FALSE, 
  rep(TRUE, 4)), 6))
MAT<-cbind(MAT,c(TRUE,FALSE,rep(TRUE,2),FALSE,TRUE)) ##
MAT<-rbind(MAT,c(TRUE,FALSE,rep(TRUE,2),FALSE,rep(TRUE,2))) ##
MAT<-cbind(MAT,c(rep(FALSE,5),rep(TRUE,2))) ##
MAT<-rbind(MAT,c(rep(FALSE,5),rep(TRUE,3))) ##

library(igraph)
cli <- cliques(graph_from_adjacency_matrix(as.matrix(MAT), mode = "undirected"))
K <- +sapply(cli, function(x) colnames(MAT) %in% names(x))
 
is.consec <- function(x) sum(x == 1) == 1 ||  all(diff(which(x == 1)) == 1) ##
K <- K[, apply(K, 2, is.consec)] ##
 
library(lpSolve)
obj <- colSums(K)^2
res <- lp("max", obj, K, "=", 1, all.bin = TRUE)
Ksoln <- K[, res$solution == 1]
o <- order(apply(Ksoln, 2, which.max))
colnames(MAT) <- Ksoln[, o] %*% 1:ncol(Ksoln)
MAT

giving:

      1     1     1     2     3     4     4     4
1  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
2  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE
3  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE
4  TRUE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE
5  TRUE FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE
6  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
7  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE
8 FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE

You can try the code below

mat <- MAT <- as.matrix(MAT)
r <- 0
while (length(mat) != 0) {
    k <- max(which(sapply(seq_along(diag(mat)), function(k) all(mat[1:k, 1:k]))))
    r <- c(r, rep(tail(r,1) + 1, k))
    mat <- as.matrix(mat[-(1:k), ][, -(1:k)])
}
colnames(MAT) <- tail(r, -1)

which gives

> MAT
        1     1     1     2     3     3
[1,] TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[2,] TRUE  TRUE  TRUE FALSE FALSE FALSE
[3,] TRUE  TRUE  TRUE FALSE  TRUE  TRUE
[4,] TRUE FALSE FALSE  TRUE FALSE  TRUE
[5,] TRUE FALSE  TRUE FALSE  TRUE  TRUE
[6,] TRUE FALSE  TRUE  TRUE  TRUE  TRUE
Related