Об инцидентности подпространства и какой-либо из вершин единичного куба
Реализован на языке Python эвристический алгоритм полиномиального времени для распознавания подпространств достаточно малой размерности, инцидентных какой-либо из вершин единичного многомерного куба (в худшем случае эта задача считается вычислительно трудной) [1, 2].
Алгоритм принимает на вход матрицу \(L\) коэффициентов линейных форм \(\ell_j\), входящих в систему уравнений $$\left\{\begin{array}{rcl} x_0&=&1\\ x_j&=&\ell_j(x_0,\dots,x_s), j>s\\ x_i&\in&\{0,1\}\\ \end{array}\right.$$ и возвращает одно из двух значений:
False, когда установлено, что подпространство не инцидентно никакой из вершин единичного куба иTrue, когда не исключена такая возможность.
(Ответ True не гарантирует, что подпространство инцидентно какой-либо вершине.)
Строки матрицы \(L\) соответствуют линейным формам. Например, для двух линейных форм \(\ell_2=2x_0+x_1\) и \(\ell_3=3x_1\) эта матрица равна $$L=\left(\begin{array}{cc} 2&1\\ 0&3\\ \end{array}\right)$$ (Матрица \(L\) содержит не все коэффициенты исходной системы линейных уравнений, а только изменяемую часть.)
Пример использования
from fractions import Fraction
import numpy as np
from qualg import may_have_01solution
inputs = [
np.array([[1, 1]]),
np.array([[2, 1],
[0, 3]]),
]
inputs.append(inputs[0] * 10**20)
inputs.append(inputs[1] * Fraction(1, 7))
for matrix_ell in inputs:
print("Input:")
print(matrix_ell)
print("Result:")
if may_have_01solution(matrix_ell):
print("A (0,1)-solution *may* exist.")
else:
print("There is *no* (0,1)-solution.")
print()
Реализация: qualg.py
from fractions import Fraction
import numpy as np
def may_have_01solution(matrix_ell):
matrix_b = compose_b(matrix_ell)
rank_a, rank_b = rank_ab(matrix_b)
return rank_a != rank_b
def compose_b(matrix_ell):
s = matrix_ell.shape[1] - 1
n = matrix_ell.shape[0] + s
height = (s + 1) * (s + 2) // 2
dtype = matrix_ell.dtype
matrix_b = np.zeros(shape=(height, n + 1), dtype=dtype)
for i in range(s):
matrix_b[i + 1, i] = -1
matrix_b[(i + 1) * (2 * s - i + 2) // 2, i] = 1
looo = np.array([1] + s * [0], dtype=dtype)
triu_indices = np.triu_indices(s + 1)
for j in range(n - s):
matrix_q = matrix_ell[j] * (matrix_ell[j] - looo).reshape(-1, 1)
matrix_q += np.tril(matrix_q, -1).T
matrix_b[:, s + j] = matrix_q[triu_indices]
matrix_b[0, -1] = 1
return matrix_b
def rank_ab(matrix_b):
matrix = matrix_b + Fraction()
height, width = matrix.shape
a_rank = b_rank = 0
for column_index in range(width):
if b_rank == height:
break
submatrix = matrix[b_rank:, column_index:]
first_column = submatrix[:, 0]
nonzero_indices = first_column.nonzero()[0]
if nonzero_indices.size == 0:
continue
if first_column[0] == 0:
first_nonzero = nonzero_indices[0]
first_row_copy = submatrix[0].copy()
submatrix[0] = submatrix[first_nonzero]
submatrix[first_nonzero] = first_row_copy
for row in range(1, height - b_rank):
if first_column[row] != 0:
ratio = first_column[row] / first_column[0]
submatrix[row] -= ratio * submatrix[0]
a_rank += column_index < width - 1
b_rank += 1
return a_rank, b_rank
Время работы
Ниже приведены иллюстрации, дающие общее представление о времени работы алгоритма. Время оценивалось на матрицах \(L\) со случайными целочисленными элементами, равномерно распределёнными на отрезке \([-10^9,10^9]\) (если не указано другое). Вычисления производились с использованием персонального компьютера (процессор Core i5-3570, 16 ГБ оперативной памяти).
Время работы алгоритма почти полностью определяется временем точного вычисления ранга матриц.
При наличии более быстрой реализации вычисления ранга, можно решать задачи значительно большей размерности.
Время работы алгоритма существенно зависит, в числе прочего, от максимальной величины коэффициентов на входе.
Публикации
- А.В. Селиверстов. Двоичные решения для больших систем линейных уравнений. Прикладная дискретная математика, 2021, № 52, стр. 5–15. DOI: 10.17223/20710410/52/1
- О.А. Зверков, А.В. Селиверстов. Эффективные нижние границы для ранга матрицы и приложения. Программирование, 2023, № 5, стр. 79–86. DOI: 10.31857/S0132347423020176