Как выполнить дискретную оптимизацию функций над матрицами?

Question

Я хотел бы оптимизировать все 30 на 30 матриц с элементами, которые являются 0 или 1. Моя целевая функция является определяющим. Одним из способов сделать это будет какой-то стохастический градиентный спуск или имитированный отжиг.

Я смотрел scipy.optimize, но, похоже, он не поддерживает такую ​​оптимизацию, насколько я могу судить. scipy.optimize.basinhopping выглядел очень заманчивым, но, похоже, ему нужны непрерывные переменные.

Есть ли какие-либо инструменты в Python для такого рода общей дискретной оптимизации?

Answer 1

Я не знаю ни одного прямого метода для дискретной оптимизации в scipy. Одним из вариантов является использование simanneal package от ПУМ или GitHub, что позволяет ввести свою собственную функцию перемещения, так что вы можете ограничить его движется в пределах своего домена:

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) 56 миллионов - это максимальное значение, когда размер матрицы равен 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.

EDIT: Чтение более тщательно с этой ссылке:

Максимальные детерминанты {1, -1} матриц с точностью до размера п = 21 приведены в следующей таблице. Размер 22 - самый маленький открытый случай. В таблице D (n) представляет собой максимальный детерминант, деленный на 2n-1. Эквивалентно, D (n) представляет собой максимальный детерминант матрицы {0, 1} размера n-1.

Чтобы таблица могла использоваться для верхних границ, но помните, что они разделены на 2n-1. Также обратите внимание, что 22 является наименьшим открытым случаем, поэтому попытка найти максимум матрицы 30x30 не была выполнена и даже не приближается к тому, чтобы быть выполненным.

2) Причина, по которой код Дэвида Цвиккера дает ответ в 30 миллионов, вероятно, объясняется тем, что он сводит к минимуму. Не максимизация.

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

Посмотрите, как у него есть знак минуса?

3) Кроме того, пространство для решения этой проблемы очень велико. Я вычислил число различных матриц как 2^(30 * 30), т. Е. Порядка 10^270. Поэтому смотреть на каждую матрицу просто невозможно, и даже смотреть на большинство из них тоже.

У меня здесь есть код (адаптированный из кода Дэвида Цвиккера), но я понятия не имею, насколько он близок к фактическому максимуму. Это займет около 45 минут, чтобы выполнить 10 миллионов итераций на моем ПК или всего около 2 минут для 1-миллиметровой итерации. Я получаю максимальное значение около 3,4 миллиарда. Но опять же, я не знаю, насколько это близко к теоретическому максимуму.

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

Это занимает около 40 секунд, чтобы запустить 1000 поколений на моем ноутбуке, который получает меня от минимального значения определителя около -5.7845x10 до -6.41504x10 . Я не очень много играл с мета-параметрами (численность населения, скорость мутации, скорость кроссовера и т. Д.), Поэтому я уверен, что можно сделать намного лучше.

---

Вот значительно улучшенная версия, которая реализует гораздо умнее функцию кроссовера, которая SWAPS блоков строк или столбцов через людей, и использует 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() 

Мой лучший результат до сих пор составляет около -3.23718x10 -3.92366x10 после 1000 поколений, которые занимают около 45 секунд на моей машине.

На основе решения cthonicdaemon связанных в комментариях, максимальный определитель для матрицы 31x31 Адамар должен быть по крайней мере 75960984159088 × 2 ~ = 8.1562x10 (это еще не доказано, что ли раствор является оптимальным). Максимальный определитель для (п-1 х N-1) бинарной матрицы 2 ^1-N раз значение для (NxN) матрицы Адамара, т.е. 8.1562x10 х 2 ^-30 ~ = 7.5961x10 , поэтому генетический алгоритм попадает в порядок величины текущего наиболее известного решения.

Однако функция фитнеса, кажется, здесь плато, и я с трудом ломаю -4x10 . Поскольку это эвристический поиск, нет никакой гарантии, что он в конечном итоге найдет глобальный оптимум.