Wie führt man eine diskrete Optimierung von Funktionen über Matrizen durch?

Question

Ich möchte über alle 30 bis 30 Matrizen mit Einträgen, die 0 oder 1 sind, optimieren. Meine Zielfunktion ist die Determinante. Ein Weg, dies zu tun, wäre eine Art stochastischer Gradientenabfall oder simuliertes Tempern.

Ich schaute auf scipy.optimize, aber es scheint nicht diese Art von Optimierung zu unterstützen, soweit ich das beurteilen kann. scipy.optimize.basinhopping sah sehr verlockend aus, aber es scheint kontinuierliche Variablen zu erfordern.

Gibt es in Python irgendwelche Werkzeuge für diese Art der allgemeinen diskreten Optimierung?

Answer 1

Ich kenne keine direkte Methode für diskrete Optimierung in scipy. Eine Alternative ist die simanneal package von pip oder Github verwenden, mit dem Sie Ihre eigene Move-Funktion einführen können, so dass Sie es einschränken können innerhalb Ihrer Domain bewegt:

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

Ich habe in diese ein wenig aussah.

Ein paar Dinge zuerst: 1) 56 Millionen ist der Maximalwert, wenn die Größe der Matrix 21x21 ist, nicht 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.

Aber das ist auch eine obere Grenze für -1, 1 Matrizen, nicht 1,0.

EDIT: Lesen genauer von diesem Link:

Die maximalen Determinanten von {1, -1} Matrizen Größe bis n = 21 sind in der folgenden Tabelle angegeben. Größe 22 ist der kleinste offene Fall. In der Tabelle repräsentiert D (n) die maximale Determinante dividiert durch 2n-1. Äquivalent stellt D (n) die maximale Determinante einer {0,1} -Matrix der Größe n-1 dar.

So kann diese Tabelle für obere Grenzen verwendet werden, aber denken Sie daran, dass sie durch 2n-1 geteilt sind. Beachten Sie auch, dass 22 der kleinste offene Fall ist. Daher wurde nicht versucht, das Maximum einer 30x30-Matrix zu finden, und ist noch nicht einmal annähernd fertig.

2) Der Grund, warum David Zwickers Code eine Antwort von 30 Millionen gibt, liegt wahrscheinlich an der Tatsache, dass er minimiert. Nicht maximierend.

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

Sehen Sie, wie er dort das Minuszeichen hat?

3) Auch der Lösungsraum für dieses Problem ist sehr groß. Ich berechne die Anzahl der verschiedenen Matrizen zu 2^(30 * 30), d.h. in der Größenordnung von 10^270. Es ist einfach unmöglich, jede Matrix zu betrachten, und selbst die meisten von ihnen sehen sie auch.

Ich habe hier ein bisschen Code (adaptiert aus David Zwickers Code), der läuft, aber ich habe keine Ahnung, wie nah es am tatsächlichen Maximum ist. Es dauert ungefähr 45 Minuten, um 10 Millionen Iterationen auf meinem PC durchzuführen, oder nur ungefähr 2 Minuten für 1 Milliteriterationen. Ich bekomme einen Maximalwert von rund 3,4 Milliarden. Aber ich habe keine Ahnung, wie nah das dem theoretischen Maximum ist.

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

Hilft das?

Answer 3

Ich denke, ein genetic algorithm könnte in diesem Fall recht gut funktionieren.Hier ist ein kurzes Beispiel zusammengeworfen deap verwendet, basierte lose auf ihrem Beispiel 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() 

Es dauert etwa 40 Sekunden 1000 Generationen auf meinem Laptop laufen zu lassen, die mich von einem minimalen Determinante Wert von etwa bekommen -5.7845x10 bis -6,41504x10 . Ich habe mit den Meta-Parametern (Bevölkerungsgröße, Mutationsrate, Crossover-Rate usw.) nicht wirklich viel herumgespielt, also bin ich mir sicher, dass es viel besser geht.

---

Hier ist eine stark verbesserte Version, die eine viel intelligentere Crossover-Funktion implementiert, die Blöcke von Zeilen oder Spalten zwischen Individuen Swaps und verwendet ein cachetools.LRUCache zu gewährleisten, dass jede Mutation Schritt eine neue Konfiguration erzeugt und die Bewertung der überspringen Determinante für Konfigurationen, die bereits versucht haben:

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

bisher mein bestes Ergebnis ist über -3.23718x10 -3.92366x10 nach 1000 Generationen, das dauert etwa 45 Sekunden auf meiner Maschine.

Auf der Basis der Lösung cthonicdaemon in den Kommentaren verknüpfen, muss die maximale Determinante für eine 31x31 Hadamard-Matrix mindestens 75960984159088 × 2 ~ = 8.1562x10 sein (ob diese Lösung bewährt noch nicht ist ist optimal). Die maximale Determinante für eine (n-1 × n-1) binäre Matrix ist der Wert für eine (n × n) Hadamard-Matrix, dh 8,1562 × 10 × 2 ^-30 ~ = 7,5961 × 10 , so kommt der genetische Algorithmus in eine Größenordnung der derzeit bekanntesten Lösung.

Allerdings scheint die Fitness-Funktion hier zu plateau, und ich habe eine harte Zeit brechen -4x10 . Da es sich um eine heuristische Suche handelt, kann nicht garantiert werden, dass sie letztendlich das globale Optimum findet.