如何在矩阵上执行函数的离散优化？
Posted
技术标签:
【中文标题】如何在矩阵上执行函数的离散优化？【英文标题】：How to perform discrete optimization of functions over matrices? 【发布时间】：2015-09-28 06:33:55 【问题描述】：我想优化所有 30 x 30 矩阵，其条目为 0 或 1。我的目标函数是行列式。一种方法是使用某种随机梯度下降或模拟退火。
我查看了scipy.optimize，但据我所知，它似乎不支持这种优化。 scipy.optimize.basinhopping 看起来很诱人，但似乎需要连续变量。
Python 中是否有任何工具可用于这种一般的离散优化？
【问题讨论】：
@ali_m 我相信现在已弃用。 @ali_m OP的问题似乎很明确。所以，如果你能弄清楚如何让盆地跳跃为它工作，那么请将它作为答案发布。据我所知，出于重要原因不推荐使用 anneal，因为它效果不佳，因此最好不要推荐它。 我想知道模拟退火的适用性，它假设有一个“附近”配置的概念会产生“附近”品质因数——本质上，你可能一开始可能会跳很多，但很快或稍后您将微调解决方案。我不知道行列式是否是这样的——我可以想象“附近”的矩阵有非常不同的行列式。不只是 SA，所有本地搜索方法都会有这个问题。另外，如果我没记错的话，会有矩阵的等价类具有相同的行列式——尝试搜索等价类。 知道有大量关于maximising determinants of 0-1 matrices的研究可能很有用。 我知道它并没有真正回答关于如何做到这一点的问题，但如果你只是在寻找最大行列式，你可以使用 this solution for -1, 1 矩阵和this transformation从-1、1到0、1的问题。 【参考方案1】：我不知道scipy 中有任何直接的离散优化方法。一种替代方法是使用 pip 或 github 中的 simanneal package，它允许您引入自己的移动功能，以便您可以将其限制为在您的域内移动：
import random
import numpy as np
import simanneal
class BinaryAnnealer(simanneal.Annealer):
def move(self):
# choose a random entry in the matrix
i = random.randrange(self.state.size)
# flip the entry 0 <=> 1
self.state.flat[i] = 1 - self.state.flat[i]
def energy(self):
# evaluate the function to minimize
return -np.linalg.det(self.state)
matrix = np.zeros((5, 5))
opt = BinaryAnnealer(matrix)
print(opt.anneal())
【讨论】：
谢谢。我用n = 20 进行了尝试，它几乎立即达到了大约 3000 万，然后停止改进。根据oeis.org/A003432，真正的最大值为 56,640,625。有没有办法让它探索更多的空间。
github 页面为您提供了一些有关如何影响模拟退火的信息。特别是，您应该使用初始温度和最终温度，以便算法探索足够多的能量景观。
谢谢我尝试了各种选项。我什至试过 opt.Tmax = 500000 opt.Tmin = 10000 opt.steps = 5000000 。它总是卡在 3000 万左右。
嗯，恐怕这与您的成本函数结构有关，即能源格局。也许可以利用行列式的结构，但我对此知之甚少。查看集群更新也可能会有所帮助，即每次迭代翻转多个条目。【参考方案2】：
我已经研究了一下。
首先要注意几点：1) 5600 万是当矩阵大小为 21x21 而不是 30x30:https://en.wikipedia.org/wiki/Hadamard%27s_maximal_determinant_problem#Connection_of_the_maximal_determinant_problems_for_.7B1.2C.C2.A0.E2.88.921.7D_and_.7B0.2C.C2.A01.7D_matrices 时的最大值。
但这也是 -1, 1 矩阵的上限，而不是 1,0。
编辑：从该链接更仔细地阅读：
下表给出了大小为 n = 21 的 1, -1 矩阵的最大行列式。 22 号是最小的开箱。在表中，D(n) 表示最大行列式除以 2n-1。等效地，D(n) 表示大小为 n-1 的 0, 1 矩阵的最大行列式。
因此该表可用于上限，但请记住它们除以 2n−1。另请注意，22 是最小的开箱，因此尚未尝试找到 30x30 矩阵的最大值，甚至还没有接近完成。
2) David Zwicker 的代码给出 3000 万个答案的原因可能是因为他正在最小化。没有最大化。
return -np.linalg.det(self.state)
看看他是怎么得到减号的？
3) 另外，这个问题的解空间很大。我计算出不同矩阵的数量为 2^(30*30)，即 10^270 的数量级。因此，查看每个矩阵根本不可能，甚至查看大多数矩阵也是如此。
我这里有一些运行的代码（改编自 David Zwicker 的代码），但我不知道它与实际最大值有多接近。在我的 PC 上进行 1000 万次迭代大约需要 45 分钟，或者 1 次迭代只需大约 2 分钟。我得到的最大值约为 34 亿。但同样，我不知道这与理论最大值有多接近。
import numpy as np
import random
import time
MATRIX_SIZE = 30
def Main():
startTime = time.time()
mat = np.zeros((MATRIX_SIZE, MATRIX_SIZE), dtype = int)
for i in range(MATRIX_SIZE):
for j in range(MATRIX_SIZE):
mat[i,j] = random.randrange(2)
print("Starting matrix:\n", mat)
maxDeterminant = 0
for i in range(1000000):
# choose a random entry in the matrix
x = random.randrange(MATRIX_SIZE)
y = random.randrange(MATRIX_SIZE)
mat[x,y] = 1 - mat[x,y]
#print(mat)
detValue = np.linalg.det(mat)
if detValue > maxDeterminant:
maxDeterminant = detValue
timeTakenStr = "\nTotal time to complete: " + str(round(time.time() - startTime, 4)) + " seconds"
print(timeTakenStr )
print(maxDeterminant)
Main()
这有帮助吗？
【讨论】：
我在能量函数中有-det(mat)的原因是因为模拟退火算法做了最小化。因此，我计算了行列式的最大值。不过，尚不清楚模拟退火是否是最大化行列式的正确方法。
谢谢。您的代码没有执行随机爬山，而只是执行随机游走，我说得对吗？
是的，这只是随机游走。所以会非常依赖起始位置。因此，可以沿着这些思路完成的其他事情是在计算行列式之前每次更改一个以上的条目（比如 10 个），或者甚至只是在每次迭代中提出一个新的随机矩阵。【参考方案3】：
我认为genetic algorithm 在这种情况下可能工作得很好。这是一个使用deap 的简单示例，大致基于他们的示例here：
import numpy as np
import deap
from deap import algorithms, base, tools
import imp
class GeneticDetMinimizer(object):
def __init__(self, N=30, popsize=500):
# we want the creator module to be local to this instance, since
# creator.create() directly adds new classes to the module's globals()
# (yuck!)
cr = imp.load_module('cr', *imp.find_module('creator', deap.__path__))
self._cr = cr
self._cr.create("FitnessMin", base.Fitness, weights=(-1.0,))
self._cr.create("Individual", np.ndarray, fitness=self._cr.FitnessMin)
self._tb = base.Toolbox()
# an 'individual' consists of an (N^2,) flat numpy array of 0s and 1s
self.N = N
self.indiv_size = N * N
self._tb.register("attr_bool", np.random.random_integers, 0, 1)
self._tb.register("individual", tools.initRepeat, self._cr.Individual,
self._tb.attr_bool, n=self.indiv_size)
# the 'population' consists of a list of such individuals
self._tb.register("population", tools.initRepeat, list,
self._tb.individual)
self._tb.register("evaluate", self.fitness)
self._tb.register("mate", self.crossover)
self._tb.register("mutate", tools.mutFlipBit, indpb=0.025)
self._tb.register("select", tools.selTournament, tournsize=3)
# create an initial population, and initialize a hall-of-fame to store
# the best individual
self.pop = self._tb.population(n=popsize)
self.hof = tools.HallOfFame(1, similar=np.array_equal)
# print summary statistics for the population on each iteration
self.stats = tools.Statistics(lambda ind: ind.fitness.values)
self.stats.register("avg", np.mean)
self.stats.register("std", np.std)
self.stats.register("min", np.min)
self.stats.register("max", np.max)
def fitness(self, individual):
"""
assigns a fitness value to each individual, based on the determinant
"""
return np.linalg.det(individual.reshape(self.N, self.N)),
def crossover(self, ind1, ind2):
"""
randomly swaps a subset of array values between two individuals
"""
size = self.indiv_size
cx1 = np.random.random_integers(0, size - 2)
cx2 = np.random.random_integers(cx1, size - 1)
ind1[cx1:cx2], ind2[cx1:cx2] = (
ind2[cx1:cx2].copy(), ind1[cx1:cx2].copy())
return ind1, ind2
def run(self, ngen=int(1E6), mutation_rate=0.3, crossover_rate=0.7):
np.random.seed(seed)
pop, log = algorithms.eaSimple(self.pop, self._tb,
cxpb=crossover_rate,
mutpb=mutation_rate,
ngen=ngen,
stats=self.stats,
halloffame=self.hof)
self.log = log
return self.hof[0].reshape(self.N, self.N), log
if __name__ == "__main__":
np.random.seed(0)
gd = GeneticDetMinimizer()
best, log = gd.run()
在我的笔记本电脑上运行 1000 代大约需要 40 秒，这使我从大约 -5.7845x108 到 -6.41504x1011 的最小行列式值。我对元参数（种群规模、突变率、交叉率等）并没有真正玩弄太多，所以我相信它可以做得更好。
这是一个大大改进的版本，它实现了一个更智能的交叉函数，可以在个体之间交换行或列块，并使用cachetools.LRUCache 来保证每个突变步骤产生一个新的配置，并跳过对行列式的评估已经尝试过的配置：
import numpy as np
import deap
from deap import algorithms, base, tools
import imp
from cachetools import LRUCache
# used to control the size of the cache so that it doesn't exceed system memory
MAX_MEM_BYTES = 11E9
class GeneticDetMinimizer(object):
def __init__(self, N=30, popsize=500, cachesize=None, seed=0):
# an 'individual' consists of an (N^2,) flat numpy array of 0s and 1s
self.N = N
self.indiv_size = N * N
if cachesize is None:
cachesize = int(np.ceil(8 * MAX_MEM_BYTES / self.indiv_size))
self._gen = np.random.RandomState(seed)
# we want the creator module to be local to this instance, since
# creator.create() directly adds new classes to the module's globals()
# (yuck!)
cr = imp.load_module('cr', *imp.find_module('creator', deap.__path__))
self._cr = cr
self._cr.create("FitnessMin", base.Fitness, weights=(-1.0,))
self._cr.create("Individual", np.ndarray, fitness=self._cr.FitnessMin)
self._tb = base.Toolbox()
self._tb.register("attr_bool", self.random_bool)
self._tb.register("individual", tools.initRepeat, self._cr.Individual,
self._tb.attr_bool, n=self.indiv_size)
# the 'population' consists of a list of such individuals
self._tb.register("population", tools.initRepeat, list,
self._tb.individual)
self._tb.register("evaluate", self.fitness)
self._tb.register("mate", self.crossover)
self._tb.register("mutate", self.mutate, rate=0.002)
self._tb.register("select", tools.selTournament, tournsize=3)
# create an initial population, and initialize a hall-of-fame to store
# the best individual
self.pop = self._tb.population(n=popsize)
self.hof = tools.HallOfFame(1, similar=np.array_equal)
# print summary statistics for the population on each iteration
self.stats = tools.Statistics(lambda ind: ind.fitness.values)
self.stats.register("avg", np.mean)
self.stats.register("std", np.std)
self.stats.register("min", np.min)
self.stats.register("max", np.max)
# keep track of configurations that have already been visited
self.tabu = LRUCache(cachesize)
def random_bool(self, *args):
return self._gen.rand(*args) < 0.5
def mutate(self, ind, rate=1E-3):
"""
mutate an individual by bit-flipping one or more randomly chosen
elements
"""
# ensure that each mutation always introduces a novel configuration
while np.packbits(ind.astype(np.uint8)).tostring() in self.tabu:
n_flip = self._gen.binomial(self.indiv_size, rate)
if not n_flip:
continue
idx = self._gen.random_integers(0, self.indiv_size - 1, n_flip)
ind[idx] = ~ind[idx]
return ind,
def fitness(self, individual):
"""
assigns a fitness value to each individual, based on the determinant
"""
h = np.packbits(individual.astype(np.uint8)).tostring()
# look up the fitness for this configuration if it has already been
# encountered
if h not in self.tabu:
fitness = np.linalg.det(individual.reshape(self.N, self.N))
self.tabu.update(h: fitness)
else:
fitness = self.tabu[h]
return fitness,
def crossover(self, ind1, ind2):
"""
randomly swaps a block of rows or columns between two individuals
"""
cx1 = self._gen.random_integers(0, self.N - 2)
cx2 = self._gen.random_integers(cx1, self.N - 1)
ind1.shape = ind2.shape = self.N, self.N
if self._gen.rand() < 0.5:
# row swap
ind1[cx1:cx2, :], ind2[cx1:cx2, :] = (
ind2[cx1:cx2, :].copy(), ind1[cx1:cx2, :].copy())
else:
# column swap
ind1[:, cx1:cx2], ind2[:, cx1:cx2] = (
ind2[:, cx1:cx2].copy(), ind1[:, cx1:cx2].copy())
ind1.shape = ind2.shape = self.indiv_size,
return ind1, ind2
def run(self, ngen=int(1E6), mutation_rate=0.3, crossover_rate=0.7):
pop, log = algorithms.eaSimple(self.pop, self._tb,
cxpb=crossover_rate,
mutpb=mutation_rate,
ngen=ngen,
stats=self.stats,
halloffame=self.hof)
self.log = log
return self.hof[0].reshape(self.N, self.N), log
if __name__ == "__main__":
np.random.seed(0)
gd = GeneticDetMinimizer(0)
best, log = gd.run()
到目前为止，我的最好成绩是在 10000 1000 代之后 -3.23718x1013 -3.92366x1013 ，在我的机器上大约需要 45 秒。
根据 cmets 中链接到的解决方案 cthonicdaemon，31x31 Hadamard 矩阵的最大行列式必须至少为 75960984159088×230 ~= 8.1562x1022（尚未证明该解决方案是否最佳）。 (n-1 x n-1) 二元矩阵的最大行列式是 (nxn) Hadamard 矩阵的值的 21-n 倍，即 8.1562x1022 x 2-30 ~= 7.5961x1013，因此遗传算法在当前最佳解决方案的一个数量级范围内。
但是，适应度函数似乎在此处趋于平稳，我很难突破 -4x1013。由于这是一种启发式搜索，因此无法保证最终会找到全局最优值。
【讨论】：
以上是关于如何在矩阵上执行函数的离散优化？的主要内容，如果未能解决你的问题，请参考以下文章