Maximale Determinante der Submatrix
Angenommen, wir haben eine quadratische Matrix M, z.
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
Ich frage mich, ob es einen effizienten Weg gibt, eine Submatrix so zu finden, dass ihre Determinante das Maximum unter allen Submatrizen ist. Die Größe der Matrix sollte größer als, 1x1aber kleiner oder gleich sein 5x5. Einige Submatrix-Beispiele sind wie folgt
> 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
Ich kann es auf Brute-Force-Weise tun, dh durch alle möglichen Untermatrizen iterieren, aber ich glaube, es muss einen Optimierungsansatz geben, der es einfacher machen kann.
Ich bevorzuge Lösungen mit CVXR, bin mir aber nicht sicher, ob dieses Optimierungsproblem konvex formuliert werden kann. Kann jemand helfen? Ansonsten sind auch andere Optimierungspakete willkommen!
Antworten
Da es vier Tage ohne Antwort waren, dachte ich, ich würde den Ball mit einer funktionierenden verallgemeinerbaren Lösung ins Rollen bringen. Leider fällt es in die Kategorie der Brute Force, obwohl es für eine 5 x 5-Matrix ziemlich schnell ist und in etwa 5 ms abgeschlossen ist:
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)
}
Das Format der Ausgabe lautet:
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
Das Problem ist natürlich, dass dies nicht gut auf größere Matrizen skaliert werden kann. Obwohl es immer noch funktioniert:
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
Ich brauche mehr als eine Sekunde, um diese Lösung für eine 10 x 10-Matrix zu finden.
Ich denke, diese Lösung ist O (n!) Komplexität, so dass Sie sie für alles vergessen können, was noch etwas größer als eine 10 x 10-Matrix ist. Ich habe das Gefühl, es sollte eine O (n³) -Lösung geben, aber meine Mathematik ist nicht gut genug, um es herauszufinden.
Ich denke, das gibt zumindest einen Maßstab, den andere mit ausgefeilteren Methoden schlagen können ...
Ich nahm Allan Camerons Lösung und verglich sie mit einer Heuristik, Threshold Accepting (TA; eine Variante des simulierten Temperns). Im Wesentlichen beginnt es mit einer zufälligen Submatrix und ändert diese Submatrix dann schrittweise, z. B. durch Austauschen von Zeilenindizes oder durch Hinzufügen oder Entfernen einer Spalte.
Eine Lösung würde als Liste mit den Zeilen- und Spaltenindizes codiert. Für eine Matrix der Größe 5x5 könnte also eine mögliche Lösung sein
x
## [[1]]
## [1] TRUE FALSE FALSE TRUE FALSE
##
## [[2]]
## [1] TRUE FALSE TRUE FALSE FALSE
Eine solche Lösung wird durch eine Nachbarschaftsfunktion geändert nb. Beispielsweise:
nb(x)
## [[1]]
## [1] TRUE FALSE FALSE TRUE TRUE
##
## [[2]]
## [1] TRUE FALSE TRUE TRUE FALSE
## ^^^^^
Bei einer solchen Lösung brauchen wir eine objektive Funktion.
OF <- function(x, M)
-det(M[x[[1]], x[[2]], drop = FALSE])
Da die Implementierung von TA, die ich verwenden werde, minimiert wird, habe ich der Determinante ein Minus vorgesetzt.
Eine Nachbarschaftsfunktion nbkönnte dies sein (obwohl sie sicherlich verbessert werden könnte):
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
}
nbWirft im Wesentlichen eine Münze und ordnet dann entweder die Zeilen- oder Spaltenindizes neu an (dh lässt die Größe der Submatrix unverändert) oder fügt eine Zeile und eine Spalte hinzu oder entfernt sie.
Schließlich erstelle ich eine Hilfsfunktion, um zufällige Anfangslösungen zu erstellen.
x0 <- function() {
k <- sample(n, 1)
x1 <- logical(n)
x1[sample(n, k)] <- TRUE
x2 <- sample(x1)
list(x1, x2)
}
Wir können Threshold Accepting ausführen. Ich verwende eine Implementierung TAopt, die im NMOFPaket enthalten ist (das ich pflege). Für einen guten Stil mache ich 10 Neustarts und behalte das beste Ergebnis.
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
Wir erhalten also die gleichen Zeilen / Spalten. Ich habe das folgende kleine Experiment durchgeführt, um die Größe von M2 auf 20 zu erhöhen . Jedes Mal, wenn ich die Lösung von TA mit der optimalen vergleiche, zeichne ich auch die Zeiten (in Sekunden) auf, die TA und die vollständige Aufzählung erfordern.
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))
}
Die Ergebnisse:
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
So, at least until size 16x16, both methods return the same result. But TA needs a constant time of less than one second (iterations are fixed at 1000).