Come eseguire l'ottimizzazione discreta delle funzioni su matrici?

Question

Mi piacerebbe ottimizzare tutte le matrici 30 per 30 con voci 0 o 1. La mia funzione obiettivo è determinante. Un modo per farlo sarebbe una sorta di discesa gradiente stocastica o ricottura simulata.

Ho guardato a scipy.optimize ma non sembra supportare questo tipo di ottimizzazione per quanto posso dire. scipy.optimize.basinhopping sembrava molto allettante ma sembra richiedere variabili continue.

Esistono strumenti in Python per questo tipo di ottimizzazione generale discreta?

Answer 1

Non conosco alcun metodo diretto per l'ottimizzazione discreta in scipy. Un'alternativa è usare la simanneal package da pip o github, che consente di introdurre la propria funzione mossa, in modo tale che si può limitare ad muove all'interno del dominio:

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

Ho guardato in questo un po '.

Un paio di cose prima: 1) 56 milioni è il valore massimo quando la dimensione della matrice è 21x21, non 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.

Ma questo è anche un limite superiore su -1, 1 matrici, non 1,0.

EDIT: lettura più attenzione da quel collegamento:

I determinanti massime di {1, -1} matrici fino alla taglia n = 21 sono riportati nella tabella seguente. La taglia 22 è la custodia aperta più piccola. Nella tabella, D (n) rappresenta il determinante massimale diviso per 2n-1. Equivalentemente, D (n) rappresenta il determinante massimale di una matrice {0, 1} di dimensione n-1.

Quindi questa tabella può essere utilizzata per i limiti superiori, ma ricorda che sono divisi per 2n-1. Si noti inoltre che 22 è il caso più piccolo aperto, quindi non è stato fatto il tentativo di trovare il massimo di una matrice 30x30, e non è ancora vicino ad essere ancora fatto.

2) Il motivo per cui il codice di David Zwicker fornisce una risposta di 30 milioni è probabilmente dovuto al fatto che sta riducendo al minimo. Non massimizzare.

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

Vedere come ha il segno meno lì?

3) Inoltre, lo spazio di soluzione per questo problema è molto grande. Calcolo il numero di matrici diverse da 2^(30 * 30), vale a dire nell'ordine di 10^270. Quindi guardare ogni matrice è semplicemente impossibile, e anche guardare la maggior parte di loro è troppo.

Ho un po 'di codice qui (adattato dal codice di David Zwicker) che gira, ma non ho idea di quanto sia vicino al massimo effettivo. Ci vogliono circa 45 minuti per fare 10 milioni di iterazioni sul mio PC, o solo circa 2 minuti per 1 millilitro di iterazioni. Ottengo un valore massimo di circa 3,4 miliardi. Ma ancora, non ho idea di quanto questo sia vicino al massimo teorico.

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

Questo aiuto?

Answer 3

Penso che un genetic algorithm potrebbe funzionare abbastanza bene in questo caso.Ecco un rapido esempio gettato insieme usando deap, liberamente ispirato loro esempio 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() 

ci vogliono circa 40 secondi per eseguire 1000 generazioni sul mio computer portatile, che mi viene da un valore minimo di circa determinante -5.7845x10 a -6.41504x10 . Non ho davvero giocato molto con i meta-parametri (dimensione della popolazione, tasso di mutazione, tasso di crossover ecc.), Quindi sono sicuro che è possibile fare molto meglio.

---

Ecco una versione migliorata che implementa una funzione crossover molto più intelligente che scambia blocchi di righe o colonne tra gli individui, e utilizza un cachetools.LRUCache per garantire che ogni passo mutazione produce una nuova configurazione, e per saltare valutazione della determinante per le configurazioni già provato:

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

Il mio miglior punteggio finora è di circa -3.23718x10 -3.92366x10 dopo 1000 generazioni, che richiede circa 45 secondi sulla mia macchina.

basato sulla soluzione cthonicdaemon legata a nei commenti, il determinante massima per una matrice 31x31 Hadamard deve essere almeno 75.960.984,159088 millions × 2 ~ = 8.1562x10 (non è ancora dimostrato se tale soluzione è ottimale). Il determinante massima per un (n-1 x n-1) matrice binaria è 2 ^1-n volte il valore di un (nxn) matrice di Hadamard, cioè 8.1562x10 x 2 ^-30 ~ = 7.5961x10 , quindi l'algoritmo genetico raggiunge un ordine di grandezza dell'attuale soluzione più nota.

Tuttavia, la funzione di fitness sembra altopiano qui intorno, e sto attraversando un periodo difficile: -4x10 . Poiché si tratta di una ricerca euristica, non vi è alcuna garanzia che alla fine troverà l'optimum globale.