Déterminant maximum de la sous-matrice
En supposant que nous ayons une matrice carrée M, par exemple,
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
Je me demande s'il existe un moyen efficace de trouver la sous-matrice a telle que son déterminant soit le maximum parmi toutes les sous-matrices. La taille de la matrice doit être supérieure 1x1mais inférieure ou égale à 5x5. Certains exemples de sous-matrice sont comme ci-dessous
> 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
Je peux le faire de manière brutale, c'est-à-dire en itérant à travers toutes les sous-matrices possibles, mais je crois qu'il doit y avoir une approche d'optimisation qui peut faciliter les choses.
Je préfère voir des solutions avec CVXRmais je ne sais pas si ce problème d'optimisation peut être formulé de manière convexe. Quelqu'un peut-il aider? Sinon, d'autres packages d'optimisation sont également les bienvenus!
Réponses
Comme cela fait quatre jours sans réponse, j'ai pensé que je lancerais le bal avec une solution généralisable fonctionnelle. Malheureusement, il entre dans la catégorie de la force brute, bien que pour une matrice 5 x 5, il soit assez rapide, se terminant en environ 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)
}
Le format de la sortie est:
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
Le problème, bien sûr, est que cela ne s'adapte pas bien à des matrices plus grandes. Bien que cela fonctionne toujours:
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
Je reçois plus d'une seconde pour trouver cette solution pour une matrice 10 x 10.
Je pense que cette solution est de complexité O (n!) , Vous pouvez donc l'oublier pour tout ce qui est même un peu plus grand qu'une matrice 10 x 10. J'ai le sentiment qu'il devrait y avoir une solution O (n³) , mais mes calculs ne sont pas assez bons pour le comprendre.
Je suppose que cela donne au moins un point de repère pour que les autres puissent battre avec des méthodes plus sophistiquées ...
J'ai pris la solution d'Allan Cameron et l'ai comparée à une heuristique, Threshold Accepting (TA; une variante de Simulated Annealing). Essentiellement, il commence par une sous-matrice aléatoire, puis modifie progressivement cette sous-matrice, par exemple en échangeant des indices de ligne, ou en ajoutant ou en supprimant une colonne.
Une solution serait codée sous forme de liste, donnant les indices de ligne et de colonne. Donc, pour une matrice de taille 5x5, une solution candidate pourrait être
x
## [[1]]
## [1] TRUE FALSE FALSE TRUE FALSE
##
## [[2]]
## [1] TRUE FALSE TRUE FALSE FALSE
Une telle solution est modifiée par une fonction de voisinage, nb. Par exemple:
nb(x)
## [[1]]
## [1] TRUE FALSE FALSE TRUE TRUE
##
## [[2]]
## [1] TRUE FALSE TRUE TRUE FALSE
## ^^^^^
Étant donné une telle solution, nous aurons besoin d'une fonction objective.
OF <- function(x, M)
-det(M[x[[1]], x[[2]], drop = FALSE])
Depuis la mise en œuvre de TA que j'utiliserai minimise, j'ai mis un moins devant le déterminant.
Une fonction de quartier nbpourrait être la suivante (bien qu'elle puisse certainement être améliorée):
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
}
Essentiellement, nbretourne une pièce, puis réorganise les indices de ligne ou de colonne (c'est-à-dire laisse la taille de la sous-matrice inchangée), ou ajoute ou supprime une ligne et une colonne.
Enfin, je crée une fonction d'assistance pour créer des solutions initiales aléatoires.
x0 <- function() {
k <- sample(n, 1)
x1 <- logical(n)
x1[sample(n, k)] <- TRUE
x2 <- sample(x1)
list(x1, x2)
}
Nous pouvons exécuter l'acceptation de seuil. J'utilise une implémentation appelée TAopt, fournie dans le NMOFpackage (que je maintiens). Pour du bon style, je fais 10 redémarrages et garde le meilleur résultat.
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
Nous obtenons donc les mêmes lignes / colonnes. J'ai effectué la petite expérience suivante, pour des tailles croissantes de M, de 2 à 20. Chaque fois, je compare la solution de TA avec la solution optimale, et j'enregistre également les temps (en secondes) requis par TA et l'énumération complète.
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))
}
Les résultats:
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
Ainsi, au moins jusqu'à la taille 16x16, les deux méthodes renvoient le même résultat. Mais TA a besoin d'un temps constant de moins d'une seconde (les itérations sont fixées à 1000).