如何对矩阵进行离散优化？

Question

我想优化所有30×30矩阵的条目是0或1。我的目标函数是行列式。一种方法是随机梯度下降或模拟退火。在

我看了^{}，但据我所知，它似乎不支持这种优化。^{}看起来很诱人，但似乎需要连续变量。在

在Python中有没有用于这种一般离散优化的工具？在

Answer 1

我不知道在scipy中有任何直接的离散优化方法。另一种选择是使用来自pip或github的^{} 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())

Answer 2

我认为在这种情况下，genetic algorithm可能会很好地工作。下面是一个使用^{}组合起来的快速示例，大致基于他们的示例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.7845x10⁸到-6.41504x10¹¹。我对元参数（种群大小、突变率、交叉率等）并没有做太多的研究，所以我相信可以做得更好。在

---

这是一个大大改进的版本，它实现了一个更智能的交叉函数，它可以在个体之间交换行或列的块，并使用^{}来保证每一个变异步骤都会产生一个新的配置，并跳过对已经尝试过的配置的行列式的评估：

^{pr2}$

到目前为止，我的最佳成绩是-3.23718x10¹³-3.92366x10¹³1000代，这在我的机器上大约需要45秒。在

根据注释中链接的解，31x31哈达玛矩阵的最大行列式必须至少为75960984159088×2³⁰~=8.1562x10²²（尚未证明该解是否最优）。（n-1×n-1）二元矩阵的最大行列式为2^1-n乘以（n×n）Hadamard矩阵的值，即8.1562x10²²x2^-30~=7.5961x10¹³，因此遗传算法在当前已知解的数量级内。在

然而，健身功能似乎在这里趋于平稳，我很难打破-4x10¹³。由于这是一种启发式搜索，因此不能保证它最终会找到全局最优。在

Answer 3

我已经调查过了。在

首先要注意的是：1）当矩阵的大小为21x21，而不是30x30时，5600万是最大值：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。在

编辑：请仔细阅读以下链接：

The maximal determinants of {1, −1} matrices up to size n = 21 are given in the following table. Size 22 is the smallest open case. In the table, D(n) represents the maximal determinant divided by 2n−1. Equivalently, D(n) represents the maximal determinant of a {0, 1} matrix of size n−1.

所以这个表可以用于上界，但要记住它们被2n−1除。还要注意的是22是最小的开放情况，所以试图找到30x30矩阵的最大值还没有完成，甚至还没有接近完成。在

2）大卫·茨威克的代码给出3000万的答案可能是因为他在最小化。不是最大化。在

return -np.linalg.det(self.state)

你看他是怎么得到负号的吗？在

3）而且这个问题的解空间很大。我计算出不同矩阵的数量为2^（30*30），即10^270的顺序。所以看每一个矩阵都是不可能的，甚至大多数矩阵也是如此。在

我不知道怎么修改代码。在我的电脑上进行1000万次迭代大约需要45分钟，或者100万次迭代只需要大约2分钟。我得到的最大值是34亿左右。但我还是不知道这离理论上的最大值有多近。在

^{pr2}$

这有帮助吗？在