i optimize on 30 30 matrices entries 0 or 1. objective function determinant. 1 way sort of stochastic gradient descent or simulated annealing.
i looked @ scipy.optimize doesn't seem support sort of optimization far can tell. scipy.optimize.basinhopping looked tempting seems require continuous variables.
are there tools in python sort of general discrete optimization?
i think genetic algorithm might work quite in case. here's quick example thrown using deap, based loosely on example here:
import numpy np import deap deap import algorithms, base, tools import imp class geneticdetminimizer(object): def __init__(self, n=30, popsize=500): # want creator module local instance, since # creator.create() directly adds new classes 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() # 'individual' consists of (n^2,) flat numpy array of 0s , 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) # 'population' consists of 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 initial population, , initialize hall-of-fame store # best individual self.pop = self._tb.population(n=popsize) self.hof = tools.halloffame(1, similar=np.array_equal) # print summary statistics 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 fitness value each individual, based on determinant """ return np.linalg.det(individual.reshape(self.n, self.n)), def crossover(self, ind1, ind2): """ randomly swaps subset of array values between 2 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() it takes 40 seconds run 1000 generations on laptop, gets me minimum determinant value of -5.7845x108 -6.41504x1011. haven't played around meta-parameters (population size, mutation rate, crossover rate etc.), i'm sure it's possible lot better.
here's improved version implements smarter crossover function swaps blocks of rows or columns across individuals, , uses cachetools.lrucache guarantee each mutation step produces novel configuration, , skip evaluation of determinant configurations have been tried:
import numpy np import deap deap import algorithms, base, tools import imp cachetools import lrucache # used control size of cache doesn't exceed system memory max_mem_bytes = 11e9 class geneticdetminimizer(object): def __init__(self, n=30, popsize=500, cachesize=none, seed=0): # 'individual' consists of (n^2,) flat numpy array of 0s , 1s self.n = n self.indiv_size = n * n if cachesize none: cachesize = int(np.ceil(8 * max_mem_bytes / self.indiv_size)) self._gen = np.random.randomstate(seed) # want creator module local instance, since # creator.create() directly adds new classes 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) # 'population' consists of 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 initial population, , initialize hall-of-fame store # best individual self.pop = self._tb.population(n=popsize) self.hof = tools.halloffame(1, similar=np.array_equal) # print summary statistics 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 have 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 individual bit-flipping 1 or more randomly chosen elements """ # ensure each mutation introduces 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 fitness value each individual, based on determinant """ h = np.packbits(individual.astype(np.uint8)).tostring() # fitness configuration if has 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 block of rows or columns between 2 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() my best score far -3.23718x1013 -3.92366x1013 after 10000 1000 generations, takes 45 seconds on machine.
based on solution cthonicdaemon linked in comments, maximum determinant 31x31 hadamard matrix must @ least 75960984159088×230 ~= 8.1562x1022 (it's not yet proven whether solution optimal). maximum determinant (n-1 x n-1) binary matrix 21-n times value (n x n) hadamard matrix, i.e. 8.1562x1022 x 2-30 ~= 7.5961x1013, genetic algorithm gets within order of magnitude of current best known solution.
however, fitness function seems plateau around here, , i'm having hard time breaking -4x1013. since it's heuristic search there no guarantee find global optimum.
Comments
Post a Comment