Determinante massimo della sotto-matrice

Sep 29 2020

Supponendo di avere una matrice quadrata M, ad esempio,

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

Mi chiedo se esista un modo efficiente per trovare la sotto-matrice in modo tale che il suo determinante sia il massimo tra tutte le sotto-matrici. La dimensione della matrice deve essere maggiore 1x1ma minore o uguale a 5x5. Alcuni esempi di sub-matrici sono come di seguito

> 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 farlo in modo brute-force, cioè iterando attraverso tutte le possibili sotto-matrici, ma credo che ci debba essere un approccio di ottimizzazione che possa renderlo più facile.

Preferisco vedere soluzioni con CVXRma non sono sicuro che questo problema di ottimizzazione possa essere formulato in modo convesso. Qualcuno può aiutare? Altrimenti, sono ben accetti anche altri pacchetti di ottimizzazione!

Risposte

7 AllanCameron Oct 02 2020 at 15:25

Dato che sono passati quattro giorni senza una risposta, ho pensato di far partire la palla con una soluzione generalizzabile funzionante. Sfortunatamente, rientra nella categoria della forza bruta, anche se per una matrice 5 x 5 è abbastanza veloce, completando in circa 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)
}

Il formato dell'output è:

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

Il problema ovviamente è che questo non si adatta bene a matrici più grandi. Anche se funziona ancora:

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

Sto ottenendo più di un secondo per trovare questa soluzione per una matrice 10 x 10.

Penso che questa soluzione sia O (n!) Complessità, quindi puoi dimenticartene per qualsiasi cosa anche un po 'più grande di una matrice 10 x 10. Ho la sensazione che dovrebbe esserci una soluzione O (n³) , ma la mia matematica non è abbastanza buona per capirla.

Immagino che almeno fornisca un punto di riferimento per gli altri da battere con metodi più sofisticati ...

4 EnricoSchumann Oct 16 2020 at 22:34

Ho preso la soluzione di Allan Cameron e l'ho confrontata con un'euristica, Threshold Accepting (TA; una variante di Simulated Annealing). Essenzialmente, inizia con una sottomatrice casuale e quindi modifica in modo incrementale questa sottomatrice, ad esempio scambiando gli indici di riga o aggiungendo o rimuovendo una colonna.

Una soluzione sarebbe codificata come un elenco, fornendo gli indici di riga e di colonna. Quindi, per una matrice di dimensioni 5x5, una soluzione candidata potrebbe essere

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

Tale soluzione viene modificata attraverso una funzione di quartiere, nb. Per esempio:

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

Data una tale soluzione, avremo bisogno di una funzione obiettivo.

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

Poiché l'implementazione dell'AT che userò riduce al minimo, ho messo un segno meno davanti al determinante.

Una funzione di quartiere nbpotrebbe essere questa (anche se potrebbe certamente essere migliorata):

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
}

Essenzialmente, nblancia una moneta e poi riorganizza gli indici di riga o colonna (cioè lascia invariata la dimensione della sottomatrice), oppure aggiungi o rimuovi una riga e una colonna.

Infine, creo una funzione di supporto per creare soluzioni iniziali casuali.

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

Possiamo eseguire l'accettazione della soglia. Uso un'implementazione chiamata TAopt, fornita nel NMOFpacchetto (che mantengo). Per un buon stile, faccio 10 ripartenze e mantengo il miglior risultato.

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

Quindi otteniamo le stesse righe / colonne. Ho eseguito il seguente piccolo esperimento, per dimensioni crescenti di M, da 2 a 20. Ogni volta confronto la soluzione di TA con quella ottimale, e registro anche i tempi (in secondi) che TA e l'enumerazione completa richiedono.

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))
}

I risultati:

     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

Quindi, almeno fino alla dimensione 16x16, entrambi i metodi restituiscono lo stesso risultato. Ma TA ha bisogno di un tempo costante inferiore a un secondo (le iterazioni sono fissate a 1000).