如何在矩阵上执行函数的离散优化？

Question

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

我查看了scipy.optimize，但据我所知，它似乎不支持这种优化。 scipy.optimize.basinhopping 看起来很诱人，但似乎需要连续变量。

Python 中是否有任何工具可用于这种一般的离散优化？

Answer 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())

Answer 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()

这有帮助吗？

Answer 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.7845x10⁸ 到 -6.41504x10¹¹ 的最小行列式值。我对元参数（种群规模、突变率、交叉率等）并没有真正玩弄太多，所以我相信它可以做得更好。

---

这是一个大大改进的版本，它实现了一个更智能的交叉函数，可以在个体之间交换行或列块，并使用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.23718x10¹³ -3.92366x10¹³ ，在我的机器上大约需要 45 秒。

根据 cmets 中链接到的解决方案 cthonicdaemon，31x31 Hadamard 矩阵的最大行列式必须至少为 75960984159088×2³⁰ ~= 8.1562x10²²（尚未证明该解决方案是否最佳）。 (n-1 x n-1) 二元矩阵的最大行列式是 (nxn) Hadamard 矩阵的值的 2^1-n 倍，即 8.1562x10²² x 2^-30 ~= 7.5961x10¹³，因此遗传算法在当前最佳解决方案的一个数量级范围内。

但是，适应度函数似乎在此处趋于平稳，我很难突破 -4x10¹³。由于这是一种启发式搜索，因此无法保证最终会找到全局最优值。