Подтвердить что ты не робот

Подсчитайте, сколько матриц имеет полный ранг для всех подматриц

Я хотел бы подсчитать, сколько m матриц, чьи элементы имеют 1 или -1, обладают тем свойством, что все его подматрицы floor(m/2)+1 by n имеют полный ранг. Мой текущий метод наивен и медленный и находится в следующем коде python/numpy. Он просто выполняет итерации по всем матрицам и проверяет все подматрицы.

import numpy as np
import itertools
from scipy.misc import comb

m = 8
n = 4

rowstochoose = int(np.floor(m/2)+1)

maxnumber = comb(m, rowstochoose, exact = True)

matrix_g=(np.array(x).reshape(m,n) for x in itertools.product([-1,1], repeat = m*n))

nofound = 0
for A in matrix_g:
    count = 0
    for rows in itertools.combinations(range(m), int(rowstochoose)):
       if (np.linalg.matrix_rank(A[list(rows)]) == int(min(n,rowstochoose))):
           count+=1
       else:
           break
    if (count == maxnumber):
         nofound+=1   
print nofound, 2**(m*n)

Есть ли лучший/более быстрый способ сделать это? Я бы хотел сделать этот расчет для n и m до 20, но любые существенные улучшения были бы большими.

Контекст. Я заинтересован в получении точных решений для https://math.stackexchange.com/info/640780/probability-that-every-vector-is-not-orthogonal-to-half-of-the-others.


В качестве точки данных для сравнения реализаций. n,m = 4,4 должен выводить 26880. n,m=5,5 для меня слишком медленно. Для n = 2 и m = 2,3,4,5,6 выходы должны быть 8, 0, 96, 0, 1280.


Текущее состояние 2 февраля 2014:

  • Ответ leewangzhong быстрый, но неверен для m > n. leewangzhong рассматривает, как это исправить.
  • Ответ Hooked не выполняется для m > n.
4b9b3361

Ответ 1

(Теперь частичное решение для n = m//2 + 1 и запрошенный код.)

Пусть k: = m//2 + 1

Это несколько эквивалентно заданию: "Сколько наборов m n-мерных векторов {-1,1} не имеют линейно-зависимых множеств размера min (k, n)?"

Для этих матриц мы знаем или можем предположить:

  • Первая запись каждого вектора равна 1 (если нет, умножьте все на -1). Это уменьшает счетчик в 2 ** м.
  • Все векторы в списке различны (если нет, любая подматрица с двумя идентичными векторами имеет не полный ранг). Это устраняет много. Выбирают (2 ** m, n) матрицы различных векторов.
  • Список векторов сортируется лексикографически (ранг не зависит от перестановок). Поэтому мы действительно думаем о наборах векторов, а не о списках. Это уменьшает счетчик в метрах! (потому что нам требуется отчетливость).

При этом мы имеем решение для n = 4, m = 8. Существует только восемь разных векторов с тем свойством, что первая запись положительна. Существует только одна комбинация (отсортированный список) из 8 различных векторов из 8 различных векторов.

array([[ 1,  1,  1,  1],
       [ 1,  1,  1, -1],
       [ 1,  1, -1,  1],
       [ 1,  1, -1, -1],
       [ 1, -1,  1,  1],
       [ 1, -1,  1, -1],
       [ 1, -1, -1,  1],
       [ 1, -1, -1, -1]], dtype=int8)

100 комбинаций размера-4 из этого списка имеют ранг 3. Таким образом, есть 0 матриц с свойством.


Для более общего решения:

Обратите внимание, что существуют векторы 2**(n-1) с первой координатой -1 и choose(2**(n-1),m) для проверки. Для n = 8 и m = 8 существует 128 векторов и 1.4297027e + 12 матриц. Это может помочь ответить: "Для я = 1,..., k, сколько комбинаций имеет ранг i?"

В качестве альтернативы, "Какие матрицы (с приведенными выше предположениями) имеют меньший ранг?" И я думаю, что ответ в точности, Достаточное условие: "Два столбца кратно друг другу". У меня такое чувство, что это правда, и я тестировал это для всех матриц 4x4, 5x5 и 6x6. (Must've прикрутил тесты). Поскольку первый столбец был выбран однородным и так как все однородные векторы кратно друг другу, любая подматрица размера k с однородным столбцом, отличным от первого столбца, будет иметь ранг меньше k.

Это не является необходимым условием. Следующая матрица является сингулярной (первая плюс четвертая равна трети плюс вторая).

array([[ 1,  1,  1,  1,  1],
       [ 1,  1,  1,  1, -1],
       [ 1,  1, -1, -1,  1],
       [ 1,  1, -1, -1, -1],
       [ 1, -1,  1, -1,  1]], dtype=int8)

Поскольку существует только два возможных значения (-1 и 1), все матрицы mxn, где m>2, k := m//2+1, n = k и с первым столбцом -1 имеют мажоритарный член в каждом столбце (т.е. по меньшей мере k членов одинаковы). Таким образом, для n = k ответ равен 0.


Для n <= 8 здесь код для генерации векторов.

from numpy import unpackbits, arange, uint8, int8

#all distinct n-length vectors from -1,1 with first entry -1
def nvectors(n):
    if n > 8:
        raise ValueError #is that the right error?
    return -1 + 2 * (
        #explode binary numbers to arrays of 8 zeroes and ones
        unpackbits(arange(2**(n-1),dtype=uint8)) #unpackbits only takes uint
            .reshape((-1,8)) #unpackbits flattens, so we need to shape it to 8 bits
            [:,-n:] #only take the last n bytes
            .view(int8) #need signed
    )

Матричный генератор:

#generate all length-m matrices that are combinations of distinct n-vectors
def matrix_g(n,m):
    return (array(mat) for mat in combinations(nvectors(n),m))

Ниже приведена функция проверки того, что все подматрицы длины maxrank имеют полный ранг. Он останавливается, если есть меньше, чем maxrank, вместо проверки всех комбинаций.

rankof = np.linalg.matrix_rank
#all submatrices of at least half size have maxrank
#(we only need to check the maxrank-sized matrices)
def halfrank(matrix,maxrank):
return all(rankof(submatr) == maxrank for submatr in combinations(matrix,maxrank))

Сгенерируем все матрицы, имеющие все полуматрицы с полным рангом   def nicematrices (m, n):       maxrank = min (m//2 + 1, n)       return (matr для matr в matrix_g (n, m), если halfrank (matr, maxrank))

Объединяя все это:

import numpy as np
from numpy import unpackbits, arange, uint8, int8, array
from itertools import combinations

#all distinct n-length vectors from -1,1 with first entry -1
def nvectors(n):
    if n > 8:
        raise ValueError #is that the right error?
    if n==0:
        return array([])
    return -1 + 2 * (
        #explode binary numbers to arrays of 8 zeroes and ones
        unpackbits(arange(2**(n-1),dtype=uint8)) #unpackbits only takes uint
            .reshape((-1,8)) #unpackbits flattens, so we need to shape it to 8 bits
            [:,-n:] #only take the last n bytes
            .view(int8) #need signed
    )

#generate all length-m matrices that are combinations of distinct n-vectors
def matrix_g(n,m):
    return (array(mat) for mat in combinations(nvectors(n),m))

rankof = np.linalg.matrix_rank
#all submatrices of at least half size have maxrank
#(we only need to check the maxrank-sized matrices)
def halfrank(matrix,maxrank):
    return all(rankof(submatr) == maxrank for submatr in combinations(matrix,maxrank))

#generate all matrices that have all half-matrices with full rank
def nicematrices(m,n):
    maxrank = min(m//2+1,n)
    return (matr for matr in matrix_g(n,m) if halfrank(matr,maxrank))

#returns (number of nice matrices, number of all matrices)
def count_nicematrices(m,n):
    from math import factorial
    return (len(list(nicematrices(m,n)))*factorial(m)*2**m, 2**(m*n))

for i in range(0,6):
    print (i, count_nicematrices(i,i))

count_nicematrices(5,5) занимает около 15 секунд для меня, подавляющее большинство из которых выполняется функцией matrix_rank.

Ответ 2

Так как никто еще не ответил, вот ответ без кода. Полезные симметрии, которые я вижу, следующие.

  • Умножьте строку на -1.
  • Умножьте столбец на -1.
  • Перенаправить строки.
  • Перенять столбцы.

Я бы атаковал эту проблему исчерпывающе генерировал неизоморфы, отфильтровывал их и суммировал размеры своих орбит. nauty будет весьма полезна для первого и третьего шага. Предполагая, что большинство матриц имеют несколько симметрий (несомненно, отличное предположение для n больших, но не очевидное априори, насколько велико), я ожидал бы, что 8x8 выполнимо, 9x9 будет пограничным, а 10x10 - вне досягаемости.

Расширенный псевдокод:

  • Создайте один представитель каждой орбиты (m - 1) матрицами (n - 1) 0-1, на которые действует группа, сгенерированная перестановками строк и столбцов, вместе с размером орбиты (= (m - 1)! (n - 1)!/размер группы автоморфизмов). Возможно, автор документа, связанный с Тимом, хотел бы поделиться своим кодексом; в противном случае см. ниже.

  • Для каждой матрицы замените записи x на (-1) ^ x. Добавьте одну строку и один столбец из 1s. Умножьте размер его орбиты на 2 ^ (m + n - 1). Это касается симметрий смены знака.

  • Отфильтруйте матрицы и суммируйте размеры орбит тех, которые остаются. Вы можете немного сократить время вычислений, выполнив сам Gram-Schmidt, чтобы при попытке всех комбинаций в лексикографическом порядке появилась возможность повторно использовать частичные результаты для общих префиксов.

Без изоморфного перечисления:

Шаблон McKay можно использовать для генерации представителей для (m + 1) матрицами n 0-1 из представителей для m с помощью n 0-1 матриц, способным к поиску глубины. С каждой матрицей m на n 0-1 свяжите двудольный граф с m черными вершинами, n белыми вершинами и соответствующим ребром для каждой 1 записи. Сделайте следующее для каждого представителя m на n.

  • Для каждого вектора длины-n постройте граф для (m + 1) с помощью n-матрицы, состоящей из представителя вместе с новым вектором, и запустите nauty, чтобы получить каноническую метку и вершинные орбит.

  • Отфильтруйте возможности, где вершина, соответствующая новому вектору, находится на другой орбите от черной вершины с наименьшим числом.

  • Отфильтруйте возможности с помощью дубликатов канонических меток.

nauty также вычисляет порядки групп автоморфизмов.

Ответ 3

Вам нужно будет переосмыслить эту проблему с математической точки зрения. При этом даже с грубой силой есть несколько программных трюков, которые вы можете использовать для ускорения процесса (так как SO - сайт программирования). Маленькие трюки, такие как пересчет int(min(n,rowstochoose)) и itertools.combinations(range(m), int(rowstochoose)), могут сэкономить несколько процентов, но реальный выигрыш возникает из воспоминаний. Другие упомянули об этом, но я подумал, что было бы полезно иметь полный, рабочий, пример кода:

import numpy as np
from scipy.misc import comb
import itertools, hashlib

m,n = 4,4

rowstochoose = int(np.floor(m/2)+1)
maxnumber    = comb(m, rowstochoose, exact = True)
combo_itr    = (x for x in itertools.product([-1,1], repeat = m*n))
matrix_itr   = (np.array(x,dtype=np.int8).reshape((n,m)) for x in combo_itr)

sub_shapes = map(list,(itertools.combinations(range(m), int(rowstochoose))))
required_rank = int(min(n,rowstochoose))

memo = {}

no_found = 0
for A in matrix_itr:
    check = True
    for s in sub_shapes:
        view  = A[s].view(np.int8)
        h    = hashlib.sha1(view).hexdigest()

        if h not in memo: 
            memo[h] =  np.linalg.matrix_rank(view)

        if memo[h] != required_rank:
            check = False
            break
    if check: no_found+=1   

print no_found, 2**(m*n)

Это дает прирост скорости почти 10 раз для случая 4x4 - вы увидите существенные улучшения для больших матриц, если вы хотите подождать достаточно долго. Это возможно для больших матриц, где стоимость ранга пропорционально более дорогая, что вы можете заказать матрицы раньше времени на хеширование:

idx  = np.lexsort(view.T)
h    = hashlib.sha1(view[idx]).hexdigest()

Для случая 4x4 это делает его немного хуже, но я ожидаю, что он изменится для случая 5x5.

Ответ 4

Алгоритм 1 - запоминание малых

Я бы использовал запоминание уже проверенных меньших матриц.

Вы можете просто записать в двоичном формате (0 для -1, 1 для 1) все меньшие матрицы. Кстати, вы можете напрямую проверять матрицы диапазонов (0 и 1) вместо (-1 и 1) - это одно и то же. Назовем эти кодирующие ИЗОБРАЖЕНИЯ. Используя длинные типы, вы можете иметь матрицы до 64 ячеек, поэтому до 8x8. Это быстро. Используя String, вы можете иметь столько, сколько вам нужно. Действительно, 8x8 более чем достаточно - в памяти 8 ГБ мы можем разместить 1G longs. это около 2 ^ 30, поэтому вы можете вспомнить матрицы размером до 25-28 элементов.

Для каждого размера у вас будет набор изображений:

для 2x2: 1001, 0110, 1000, 0100, 0010, 0001, 0111, 1011, 1101, 1110.

Итак, у вас будет archive= массив из NxN, каждый элемент которого будет упорядоченным списком двоичных изображений хороших матриц. - (для размера матрицы MxN, где M >= N, соответствующее место в архиве будет иметь координаты M, N. Если M

  • Когда вы проверяете новую большую матрицу, разделите ее на маленькие
  • Для каждой малой матрицы T
    • Если подходящее место в архиве для размера T не имеет списка, создайте его и заполните изображения всех полноразмерных матриц размера T и изображений заказов. Если вы потеряли память, остановите процесс заполнения архива.
    • Если T может быть в архиве, в соответствии с размером:
      • Сделать изображение T
      • Посмотрите на изображение (t) в списке - если оно в нем, это нормально, если нет, следует отбросить большую матрицу.
    • Если T слишком велико для архива, проверьте его, как вы это делаете.

Алгоритм 2 - увеличение размеров

Другая возможность - создать более крупные матрицы, добавив части к меньшим, уже найденным. Вы должны решить, до какого размера матрицы будут расти.

Когда вы найдете "правильную" матрицу размера MxN, попробуйте добавить строку вверх. Новые матрицы следует проверять только для подматриц, которые включают новую строку. То же самое с новым столбцом.

Вы должны установить точный алгоритм, из которого производятся размеры. Таким образом, вы можете свести к минимуму количество запоминаемых матриц. Я думал об этой последовательности:

  • Начните с матриц 2x2.
  • продолжить с 3x2
  • 4x2, 3x3
  • 5x2, 4x3
  • 6x2, 5x3, 4x4
  • 7x2, 6x3, 5x4
  • ...

Итак, вы можете помнить только (M + N)/2-1 матрицы для поиска среди размеров MxN.

Если каждый раз, когда мы можем создать новый размер из двух старых, мы получаем более квадратные, мы могли бы также значительно избавиться от запоминания матриц: для "длинных" матриц как 7x2 нам нужно помнить и проверять только последние строка 1x2. Для матриц 6x3 мы должны помнить их заглушку 2x3 и т.д.

Кроме того, вам не нужно запоминать самые большие матрицы - вы не будете использовать их для дальнейшего подсчета.

Снова используйте "изображения" для запоминания матрицы.