我采用了 Allan Cameron 的解决方案,并将其与启发式阈值接受(TA;模拟退火的一种变体)进行了比较。本质上,它从一个随机子矩阵开始,然后逐渐改变这个子矩阵,例如通过交换行索引,或者通过添加或删除一列。
解决方案将被编码为列表,给出行和列索引。因此,对于大小为 5x5 的矩阵,一个候选解决方案可能是
x
## [[1]]
## [1] TRUE FALSE FALSE TRUE FALSE
##
## [[2]]
## [1] TRUE FALSE TRUE FALSE FALSE
这样的解决方案是通过邻域函数来改变的,nb
。例如:
nb(x)
## [[1]]
## [1] TRUE FALSE FALSE TRUE TRUE
##
## [[2]]
## [1] TRUE FALSE TRUE TRUE FALSE
## ^^^^^
给定这样的解决方案,我们将需要一个目标函数。
OF <- function(x, M)
-det(M[x[[1]], x[[2]], drop = FALSE])
自从实施 TA 以来,我将使用最小值,我在行列式前面加上了一个减号。
邻里功能nb
可能是这样的(尽管它肯定可以改进):
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
}
本质上,nb
掷硬币然后重新排列行或列索引(即保持子矩阵的大小不变),或者添加或删除一行和一列。
最后,我创建了一个辅助函数来创建随机初始解决方案。
x0 <- function() {
k <- sample(n, 1)
x1 <- logical(n)
x1[sample(n, k)] <- TRUE
x2 <- sample(x1)
list(x1, x2)
}
我们可以运行阈值接受。TAopt
我使用包中提供的一个名为 的实现NMOF
(我维护)。为了获得良好的风格,我会重新启动 10 次并保持最佳结果。
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
所以我们得到相同的行/列。我进行了以下小实验,将 的大小M
从 2 增加到 20。每次我将 TA 的解决方案与最佳解决方案进行比较,我还记录了 TA 和完整枚举所需的时间(以秒为单位)。
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))
}
结果:
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
因此,至少在大小为 16x16 之前,两种方法都返回相同的结果。但是 TA 需要不到一秒的恒定时间(迭代次数固定为 1000)。