Determinante máximo de submatriz

Sep 29 2020

Suponiendo que tenemos una matriz cuadrada M, por ejemplo,

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

Me pregunto si hay una manera eficiente de encontrar la submatriz a tal que su determinante sea el máximo entre todas las submatrices. El tamaño de la matriz debe ser mayor 1x1pero menor o igual que 5x5. Algunos ejemplos de submatriz son los siguientes

> 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

Puedo hacerlo de manera de fuerza bruta, es decir, iterando a través de todas las submatrices posibles, pero creo que debe haber algún enfoque de optimización que pueda tomarlo con más calma.

Prefiero ver soluciones con CVXRpero no estoy seguro de si este problema de optimización se puede formular de manera convexa. ¿Alguien puede ayudar? De lo contrario, también son bienvenidos otros paquetes de optimización.

Respuestas

7 AllanCameron Oct 02 2020 at 15:25

Como han pasado cuatro días sin una respuesta, pensé que podría poner en marcha la pelota con una solución generalizable funcional. Desafortunadamente, cae en la categoría de fuerza bruta, aunque para una matriz de 5 x 5 es bastante rápido, y se completa en aproximadamente 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)
}

El formato de la salida es:

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

El problema, por supuesto, es que esto no se adapta bien a matrices más grandes. Aunque todavía funciona:

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

Me estoy dando más de un segundo para encontrar esta solución para una matriz de 10 x 10.

Creo que esta solución es de complejidad O (n!) , Por lo que puede olvidarse de ella para cualquier cosa incluso un poco más grande que una matriz de 10 x 10. Tengo la sensación de que debería haber una solución O (n³) , pero mis matemáticas no son lo suficientemente buenas para resolverlo.

Supongo que al menos da un punto de referencia para que otros lo superen con métodos más sofisticados ...

4 EnricoSchumann Oct 16 2020 at 22:34

Tomé la solución de Allan Cameron y la comparé con una heurística, Aceptación de umbral (TA; una variante de Recocido simulado). Esencialmente, comienza con una submatriz aleatoria y luego cambia gradualmente esta submatriz, por ejemplo, intercambiando índices de fila o agregando o quitando una columna.

Una solución se codificaría como una lista, dando los índices de fila y columna. Entonces, para una matriz de tamaño 5x5, una solución candidata podría ser

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

Solución tal se cambia a través de una función de vecindad, nb. Por ejemplo:

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

Dada tal solución, necesitaremos una función objetivo.

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

Dado que la implementación de TA que usaré minimiza, he puesto un signo menos delante del determinante.

Una función de vecindario nbpodría ser esta (aunque ciertamente podría mejorarse):

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
}

Esencialmente, nblanza una moneda y luego reorganiza los índices de fila o columna (es decir, deja el tamaño de la submatriz sin cambios), o agrega o quita una fila y una columna.

Finalmente, creo una función auxiliar para crear soluciones iniciales aleatorias.

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

Podemos ejecutar Threshold Accepting. Utilizo una implementación llamada TAopt, proporcionada en el NMOFpaquete (que mantengo). Para tener buen estilo, hago 10 reinicios y mantengo el mejor 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

Entonces obtenemos las mismas filas / columnas. Ejecuté el siguiente pequeño experimento, para aumentar los tamaños de M, de 2 a 20. Cada vez comparo la solución de TA con la óptima y también registro los tiempos (en segundos) que requieren TA y la enumeración completa.

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

Los 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

Entonces, al menos hasta el tamaño 16x16, ambos métodos devuelven el mismo resultado. Pero TA necesita un tiempo constante de menos de un segundo (las iteraciones se fijan en 1000).