Determinante máximo da submatriz

Question

Assumindo que temos uma matriz quadrada M, por exemplo,

set.seed(1)
M <- matrix(rnorm(5*5), 5, 5)

> M
           [,1]       [,2]       [,3]        [,4]        [,5]
[1,] -0.6264538 -0.8204684  1.5117812 -0.04493361  0.91897737
[2,]  0.1836433  0.4874291  0.3898432 -0.01619026  0.78213630
[3,] -0.8356286  0.7383247 -0.6212406  0.94383621  0.07456498
[4,]  1.5952808  0.5757814 -2.2146999  0.82122120 -1.98935170
[5,]  0.3295078 -0.3053884  1.1249309  0.59390132  0.61982575

Estou me perguntando se existe uma maneira eficiente de encontrar a submatriz a de modo que seu determinante seja o máximo entre todas as submatrizes. O tamanho da matriz deve ser maior, 1x1mas menor ou igual a 5x5. Alguns exemplos de submatriz são como abaixo

> M[c(1,5),c(2,3)]
           [,1]     [,2]
[1,] -0.8204684 1.511781
[2,] -0.3053884 1.124931

> M[c(1,2,4),c(1,4,5)]
           [,1]        [,2]       [,3]
[1,] -0.6264538 -0.04493361  0.9189774
[2,]  0.1836433 -0.01619026  0.7821363
[3,]  1.5952808  0.82122120 -1.9893517

> M[1:4,2:5]
           [,1]       [,2]        [,3]        [,4]
[1,] -0.8204684  1.5117812 -0.04493361  0.91897737
[2,]  0.4874291  0.3898432 -0.01619026  0.78213630
[3,]  0.7383247 -0.6212406  0.94383621  0.07456498
[4,]  0.5757814 -2.2146999  0.82122120 -1.98935170

Posso fazer isso de uma maneira de força bruta, ou seja, iterando por todas as submatrizes possíveis, mas acredito que deve haver alguma abordagem de otimização para facilitar.

Prefiro ver soluções com, CVXRmas não tenho certeza se esse problema de otimização pode ser formulado de maneira convexa. Alguém pode ajudar? Caso contrário, outros pacotes de otimização também são bem-vindos!

Answer 1

Como já se passaram quatro dias sem uma resposta, pensei em fazer a bola rolar com uma solução generalizável funcional. Infelizmente, ele se enquadra na categoria de força bruta, embora para uma matriz 5 x 5 seja bastante rápido, completando em cerca de 5 ms:

max_det <- function(M) {
  if(diff(dim(M)) != 0) stop("max_det requires a square matrix")
  
  s  <- lapply(seq(dim(M)[1])[-1], function(x) combn(seq(dim(M)[1]), x))
  
  all_dets <- lapply(s, function(m) {
    apply(m, 2, function(i) apply(m, 2, function(j) det(M[j, i])))
    })
  
  i <- which.max(sapply(all_dets, max))
  subs <- which(all_dets[[i]] == max(all_dets[[i]]), arr.ind = TRUE)

  sub_M <- M[s[[i]][,subs[1]], s[[i]][,subs[2]]]
  
  list(max_determinant = det(sub_M),
       indices = list(rows = s[[i]][,subs[1]], columns = s[[i]][,subs[2]]),
       submatrix = sub_M)
}

O formato da saída é:

max_det(M)
#> $max_determinant #> [1] 4.674127 #> #> $indices
#> $indices$rows
#> [1] 3 4 5
#> 
#> $indices$columns
#> [1] 1 3 4
#> 
#> 
#> $submatrix
#>            [,1]       [,2]      [,3]
#> [1,] -0.8356286 -0.6212406 0.9438362
#> [2,]  1.5952808 -2.2146999 0.8212212
#> [3,]  0.3295078  1.1249309 0.5939013

O problema, claro, é que isso não é bem dimensionado para matrizes maiores. Embora ainda funcione:

set.seed(1)
M <- matrix(rnorm(10 * 10), 10, 10)

#> max_det(M)
#> $max_determinant
#> [1] 284.5647
#> 
#> $indices #> $indices$rows #> [1] 1 3 4 5 6 8 9 10 #> #> $indices$columns #> [1] 2 3 4 6 7 8 9 10 #> #> #> $submatrix
#>             [,1]        [,2]        [,3]       [,4]        [,5]         [,6]
#> [1,]  1.51178117  0.91897737  1.35867955  0.3981059  2.40161776  0.475509529
#> [2,] -0.62124058  0.07456498  0.38767161  0.3411197  0.68973936  0.610726353
#> [3,] -2.21469989 -1.98935170 -0.05380504 -1.1293631  0.02800216 -0.934097632
#> [4,]  1.12493092  0.61982575 -1.37705956  1.4330237 -0.74327321 -1.253633400
#> [5,] -0.04493361 -0.05612874 -0.41499456  1.9803999  0.18879230  0.291446236
#> [6,]  0.94383621 -1.47075238 -0.05931340 -1.0441346  1.46555486  0.001105352
#> [7,]  0.82122120 -0.47815006  1.10002537  0.5697196  0.15325334  0.074341324
#> [8,]  0.59390132  0.41794156  0.76317575 -0.1350546  2.17261167 -0.589520946
#>            [,7]       [,8]
#> [1,] -0.5686687 -0.5425200
#> [2,]  1.1780870  1.1604026
#> [3,] -1.5235668  0.7002136
#> [4,]  0.5939462  1.5868335
#> [5,]  0.3329504  0.5584864
#> [6,] -0.3041839 -0.5732654
#> [7,]  0.3700188 -1.2246126
#> [8,]  0.2670988 -0.4734006

Estou demorando um segundo para encontrar essa solução para uma matriz 10 x 10.

Eu acho que essa solução é O (n!) Complexidade, então você pode esquecê-la para qualquer coisa um pouco maior do que uma matriz 10 x 10. Tenho a sensação de que deveria haver uma solução O (n³) , mas minha matemática não é boa o suficiente para descobri-la.

Eu acho que pelo menos dá uma referência para outros baterem com métodos mais sofisticados ...

Answer 2

Eu peguei a solução de Allan Cameron e a comparei com uma heurística, Aceitação de Limiar (TA; uma variante do Simulated Annealing). Essencialmente, ele começa com uma submatriz aleatória e, em seguida, altera essa submatriz de forma incremental, por exemplo, trocando índices de linha ou adicionando ou removendo uma coluna.

Uma solução seria codificada como uma lista, fornecendo os índices de linha e coluna. Portanto, para uma matriz de tamanho 5x5, uma solução candidata pode ser

x
## [[1]]
## [1]  TRUE FALSE FALSE  TRUE FALSE
## 
## [[2]]
## [1]  TRUE FALSE  TRUE FALSE FALSE

Essa solução é alterada por meio de uma função de vizinhança nb,. Por exemplo:

nb(x)
## [[1]]
## [1]  TRUE FALSE FALSE  TRUE  TRUE
## 
## [[2]]
## [1]  TRUE FALSE  TRUE  TRUE FALSE
##                       ^^^^^

Dada essa solução, precisaremos de uma função objetivo.

OF <- function(x, M)
    -det(M[x[[1]], x[[2]], drop = FALSE])

Desde a implementação do TA, usarei minimiza, coloquei um sinal de menos na frente do determinante.

Uma função de bairro nbpoderia ser esta (embora certamente pudesse ser melhorada):

nb <- function(x, ...) {
    if (sum(x[[1L]]) > 0L &&
        sum(x[[1L]]) < length(x[[1L]]) &&
        runif(1) > 0.5) {
        rc <- if (runif(1) > 0.5)
                  1 else 2
        select1 <- which( x[[rc]])
        select2 <- which(!x[[rc]])
        size <- min(length(select1), length(select2))
        size <- sample.int(size, 1)
        i <- select1[sample.int(length(select1), size)]
        j <- select2[sample.int(length(select2), size)]
        x[[rc]][i] <- !x[[rc]][i]
        x[[rc]][j] <- !x[[rc]][j]        
    } else {            
        i <- sample.int(length(x[[1L]]), 1)
        if (x[[1L]][i]) {
            select <- which( x[[2L]])
        } else {
            select <- which(!x[[2L]])
        }
        j <- select[sample.int(length(select), 1)]
        x[[1L]][i] <- !x[[1L]][i]
        x[[2L]][j] <- !x[[2L]][j]
    }
    x
}

Essencialmente, nbjogue uma moeda no ar e reorganize os índices de linha ou coluna (ou seja, deixe o tamanho da submatriz inalterado), ou adicione ou remova uma linha e uma coluna.

Por fim, crio uma função auxiliar para criar soluções iniciais aleatórias.

x0 <- function() {
    k <- sample(n, 1)
    x1 <- logical(n)
    x1[sample(n, k)] <- TRUE
    x2 <- sample(x1)
    list(x1, x2)
}

Podemos executar a aceitação de limite. Eu uso uma implementação chamada TAopt, fornecida no NMOFpacote (que mantenho). Para um bom estilo, faço 10 reinicializações e mantenho o melhor resultado.

n <- 5
M <- matrix(rnorm(n*n), n, n)
max_det(M)$indices ## $rows
## [1] 1 2 4
## 
## $columns ## [1] 2 3 5 library("NMOF") restartOpt(TAopt, 10, OF, list(x0 = x0, neighbour = nb, printBar = FALSE, printDetail = FALSE, q = 0.9, nI = 1000, drop0 = TRUE), M = M, best.only = TRUE)$xbest
## [[1]]
## [1]  TRUE  TRUE FALSE  TRUE FALSE
## 
## [[2]]
## [1] FALSE  TRUE  TRUE FALSE  TRUE

Portanto, obtemos as mesmas linhas / colunas. Executei o seguinte pequeno experimento, para tamanhos crescentes de M, de 2 a 20. Cada vez, comparo a solução de TA com a ótima e também registro os tempos (em segundos) que TA e a enumeração completa exigem.

set.seed(134345)
message(format(c("Size",
        "Optimum",
        "TA",
        "Time optimum",
        "Time TA"), width = 13, justify = "right"))
for (i in 2:20) {
    n <- i
    M <- matrix(rnorm(n*n), n, n)
    t.opt <- system.time(opt <- max_det(M)$max_determinant) t.ta <- system.time(ta <- -restartOpt(TAopt, 10, OF, list(x0 = x0, neighbour = nb, printBar = FALSE, printDetail = FALSE, q = 0.9, nI = 1000, drop0 = TRUE), M = M, best.only = TRUE)$OFvalue)

    message(format(i, width = 13),
            format(round(opt, 2),  width = 13),
            format(round(ta, 2),  width = 13),
            format(round(t.opt[[3]],1), width = 13),
            format(round(t.ta[[3]],1), width = 13))
}

Os resultados:

     Size      Optimum           TA Time optimum      Time TA
        2           NA         1.22            0          0.7
        3         1.46         1.46            0          0.6
        4         2.33         2.33            0          0.7
        5        11.75        11.75            0          0.7
        6         9.33         9.33            0          0.7
        7          9.7          9.7            0          0.7
        8       126.38       126.38          0.1          0.7
        9         87.5         87.5          0.3          0.7
       10       198.63       198.63          1.3          0.7
       11      1019.23      1019.23          5.1          0.7
       12     34753.64     34753.64           20          0.7
       13     16122.22     16122.22         80.2          0.7
       14     168943.9     168943.9        325.3          0.7
       15     274669.6     274669.6       1320.8          0.7
       16      5210298      5210298       5215.4          0.7

Portanto, pelo menos até o tamanho 16x16, os dois métodos retornam o mesmo resultado. Mas o TA precisa de um tempo constante de menos de um segundo (as iterações são fixadas em 1000).