Программирование, 2023, № 5, стр. 79-86

1. ВВЕДЕНИЕ

Python – современный язык программирования общего назначения, вышедший в последние годы на первое место в различных рейтингах популярности и особенно востребованный в научных приложениях. Сильными сторонами Python считаются легкость начального освоения, быстрота процесса разработки, лаконичность и удобочитаемость текстов программ. В качестве основного недостатка обычно называется относительно низкая производительность программ, обусловленная в основном интерпретируемой природой языка и динамической типизацией. Однако этот недостаток во многом компенсируется доступностью большого разнообразия библиотек, предоставляющих Python-разработчикам доступ к эффективным реализациям различных алгоритмов, а также возможностью подключения собственных программных модулей, написанных на компилируемых языках [1]. В контексте нашей задачи удобной особенностью Python является встроенная в язык поддержка работы с целыми числами неограниченной битовой длины, а также с рациональными числами.

Ранг n × n матрицы над полем можно вычислить, используя полиномиальное число процессоров и выполнив на каждом из них лишь $O(\mathop {\log }\nolimits_2^2 n)$ операций над этим полем [2, 3]. Эффективное распараллеливание возможно не только для вычисления ранга, но также для вычисления обратной матрицы, $LDU$-разложения квадратной матрицы и разложения Холецкого симметричной положительно определенной матрицы [4, 5]. С другой стороны, сложность вычисления ранга [6] и характеристического многочлена [7–9] близка к сложности матричного умножения. Для циркулянтов известны легко проверяемые условия невырожденности [10]. Теоретический интерес представляет задача проверки линейной независимости системы многочленов от нескольких переменных. Иногда такие многочлены могут быть заданы короткими выражениями или арифметическими схемами, хотя матрицы коэффициентов таких многочленов имеют большой размер. Поэтому важны легко проверяемые оценки ранга матрицы. Обычно вычисления выполняются над полем рациональных чисел, но символьные вычисления можно проводить и над полями алгебраических чисел [11, 12].

Точки детерминантального многообразия соответствуют (с точностью до умножения на ненулевое число) матрицам фиксированного размера, ранг которых ограничен сверху. Детерминантальное многообразие может иметь очень большую степень [13, 14]. Поэтому над алгебраически замкнутым полем нельзя проверить значение ранга матрицы посредством неветвящейся программы небольшого размера. С другой стороны, на этих многообразиях лежат линейные подпространства, что иногда можно использовать для оценки ранга [15].

Используя оценки ранга, мы рассмотрим эвристический алгоритм проверки инцидентности данного подпространства и какой-либо вершины единичного куба. Такая проверка эквивалентна поиску (0, 1)-решения для системы уравнений

где через ${{\ell }_{j}}$ обозначены линейные формы [16].

Хотя уже известны эвристические алгоритмы для решения этой задачи при дополнительных ограничениях [17–19], принято считать, что эта задача трудная для вычисления в худшем случае. Поэтому актуально создание новых алгоритмов, у которых вычислительная сложность в среднем ограничена многочленом от длины входа при других ограничениях.

Говоря о задачах распознавания, мы предполагаем три варианта ответа: вход может быть не только принят или отвергнут, но также возможно явное уведомление о неопределенности выбора. При этом ответ должен быть получен за конечное время и без ошибок, а при выполнении легко проверяемого условия уведомление о неопределенности может быть выдано лишь на малой доле входов среди всех входов данной длины [20, 21]. При оценке алгебраической сложности достаточно требовать, чтобы уведомление о неопределенности выдавалось на множестве входов, на которых обращается в нуль некоторый многочлен, отличный от тождественно нулевого [16, 22].

В разделе 2 представлены теоретические результаты. В разделе 3 описана реализация алгоритма на языке Python. В разделе 4 дано краткое заключение.

2. ТЕОРЕТИЧЕСКИЕ РЕЗУЛЬТАТЫ

Элементами (0, 1)-матрицы служат лишь нули или единицы. Подматрицей называется матрица, полученная удалением некоторых строк и столбцов. Матрица размера $n \times n$ называется вполне неразложимой, когда в ней нет нулевой подматрицы размера $s \times (n - s)$ ни для какого s.

Рассмотрим квадратные матрицы, имеющие в каждой строке и каждом столбце по два ненулевых элемента. Если такого типа (0, 1)-матрица вполне неразложимая, то ее перманент равен двум [23]. В этом случае значение определителя принадлежит множеству $\{ - 2,0,2\} $. Очевидно, вырождена любая такая (0, 1)-матрица над полем характеристики два. Примером вырожденной матрицы над произвольным полем служит 2 × 2 матрица ранга один, каждый элемент которой равен единице. Другим примером служит (0, 1)-циркулянт

Эта матрица тоже вырожденная, но ее ранг равен трем. Характеристический многочлен этой матрицы равен ${{x}^{4}} - 4{{x}^{3}} + 6{{x}^{2}} - 4x$. Мы дадим нижнюю оценку ранга таких матриц.

Теорема 1. Дана вполне неразложимая $n \times n$ матрица $A$, где в каждой строке и каждом столбце ровно по два ненулевых элемента. Ранг матрицы $A$ ограничен снизу: ${\text{rank}}(A) \geqslant n - 1$.

Доказательство. Условие полной неразложимости не нарушается при перестановке строк или столбцов. Поэтому без ограничения общности полагаем, что на главной диагонали матрицы $A$ все элементы ненулевые. Тогда матрица $A$ равна сумме двух невырожденных матриц $A = D + CP$, где обе матрицы $C$ и $D$ диагональные и невырожденные, а через $P$ обозначена матрица перестановки с нулями на главной диагонали. Поскольку матрица $A$ вполне неразложимая, перестановка $P$ состоит из одного цикла длины $n$.

Для любой матрицы перестановки $Q$ матрица $QA{{Q}^{{ - 1}}}$ получается из матрицы $A$ согласованной перестановкой строк и столбцов. При этом элементы на главной диагонали остаются на главной диагонали. Поэтому для некоторой матрицы перестановки $Q$ в матрице $QA{{Q}^{{ - 1}}}$ отличны от нуля элементы на главной диагонали, непосредственно над главной диагональю и еще один элемент в первом столбце и последней строке. Такая матрица $QA{{Q}^{{ - 1}}}$ имеет вид (для $n > 6$)

где символом $*$ отмечены ненулевые элементы матрицы. Поскольку при перестановке строк или столбцов ранг матрицы не меняется, ранг матрицы $QA{{Q}^{{ - 1}}}$ совпадает с рангом матрицы $A$. Для завершения доказательства теоремы достаточно найти невырожденную $(n - 1) \times (n - 1)$ подматрицу $B$ в матрице $QA{{Q}^{{ - 1}}}$. Выберем подматрицу $B$, удаляя из матрицы $A$ первую строку и первый столбец. Эта матрица верхнетреугольная. Каждый элемент на главной диагонали матрицы $B$ отличен от нуля. Поэтому матрица $B$ невырожденная.                    ◻

В теореме 1 условие полной неразложимости важно. Например, для четного $n$ у блочно-диагональной матрицы размера $n \times n$, на диагонали которой стоят 2 × 2 блоки по четыре единицы в каждом, ранг равен $n{\text{/}}2$.

Далее мы рассматриваем поля характеристики нуль. Рассмотрим задачу об инцидентности данного аффинного подпространства и какой-либо вершины $n$-мерного куба. Метод состоит в попытке построить алгебраическую гиперповерхность малой степени, проходящую через каждую вершину куба, но не пересекающую данное аффинное подпространство. Вычислительная сложность этого метода зависит от степени гиперповерхности. С другой стороны, этот метод эффективен лишь при условии, что размерность аффинного подпространства достаточно мала. Пусть это подпространство определяется системой из $m$ линейных уравнений. Размерность линейного пространства форм степени $d$ от переменных ${{x}_{0}}$, …, ${{x}_{n}}$ равна биномиальному коэффициенту $\left( {\begin{array}{*{20}{c}} {n + d} \\ d \end{array}} \right)$. Поэтому рассматриваемый метод эффективен в среднем, когда число уравнений удовлетворяет неравенству типа m ≥ n – $\sqrt {d(d - 1)n - o(n)} $, где через $d$ обозначена степень искомой гиперповерхности. Мы ограничимся случаем d = 2, когда гиперповерхностью служит квадрика [16].

Пусть $n$-мерный куб вложен в проективное пространство, а однородные координаты вершин принадлежат множеству $\{ 0,1\} $, но ни одна из вершин куба не лежит на бесконечно удаленной гиперплоскости ${{x}_{0}} = 0$. Известно [24], что квадрика, проходящая через каждую вершину этого куба, определяется квадратичной формой типа

Если подпространство задано системой уравнений ${{x}_{j}} = {{\ell }_{j}}({{x}_{0}}, \ldots ,{{x}_{{n - m}}})$, где $n - m < j \leqslant n$, то существование квадрики, проходящей через каждую вершину куба и пересекающей это подпространство лишь на бесконечности, эквивалентно разрешимости уравнения

относительно переменных ${{\lambda }_{1}}$, …, ${{\lambda }_{n}}$. В свою очередь, эта задача эквивалентна разрешимости системы линейных уравнений. По теореме Кронекера–Капелли, эта проверка сводится к сравнению рангов двух матриц. Алгоритм 1 не обеспечивает полиномиальную вычислительную сложность в худшем случае. Но иногда он гораздо эффективнее полного перебора вершин куба.

В алгоритме 1 первые $n - m$ столбцов матрицы $A$ всегда содержат по два ненулевых элемента 1 и –1. В столбце с номером $j > n - m$ элементы равны многочленам второй степени от коэффициентов линейной формы ${{\ell }_{j}}$.

Алгоритм 1. Проверка существования вершины куба, инцидентной данному подпространству.

Input: целые числа $0 < m < n$ и линейные формы ${{\ell }_{j}}({{x}_{0}}, \ldots ,{{x}_{{n - m}}})$, где $n - m < j \leqslant n$.

1: Элементами матрицы $A$ служат коэффициенты формы

где строки матрицы $A$ соответствуют мономам второй степени от переменных ${{x}_{0}}$, …, ${{x}_{{n - m}}}$, а столбцы – переменным ${{\lambda }_{1}}$, …, ${{\lambda }_{n}}$;

2: Расширенная матрица $B$ получается из матрицы $A$ добавлением столбца, в котором единица стоит в строке, соответствующей моному $x_{0}^{2}$, а остальные элементы этого столбца равны нулю;

3: if ${\text{rank}}(A) = {\text{rank}}(B)$ then вход отвергается

4: else уведомление о неопределенности.

Рассмотрим пример, соответствующий прямой на плоскости. Пусть $n = 2$, $m = 1$, линейная форма ${{\ell }_{2}} = {{x}_{0}} + {{x}_{1}}$. Получаем уравнение

Здесь три монома от переменных ${{x}_{0}}$ и ${{x}_{1}}$. Это $x_{0}^{2}$, ${{x}_{0}}{{x}_{1}}$ и $x_{1}^{2}$. Соответствующая матрица $A$ содержит два столбца и три строки:

Расширенная матрица равна

Здесь ${\text{rank}}(A) = 2$ и ${\text{rank}}(B) = 3$. Поэтому искомые ${{\lambda }_{1}}$ и ${{\lambda }_{2}}$ не существуют. Ответ алгоритма 1 состоит в уведомлении о неопределенности.

Рассмотрим другой пример, соответствующий прямой в пространстве. Пусть $n = 3$, $m = 2$, линейные формы ${{\ell }_{2}} = 2{{x}_{0}} + {{x}_{1}}$ и ${{\ell }_{3}} = 3{{x}_{1}}$. Получаем уравнение ${{\lambda }_{1}}( - {{x}_{0}}{{x}_{1}} + x_{1}^{2})\,\, + $ ${{\lambda }_{2}}(2x_{0}^{2} + 3{{x}_{0}}{{x}_{1}} + x_{1}^{2}) + $ ${{\lambda }_{3}}( - 3{{x}_{0}}{{x}_{1}} + 9x_{1}^{2}) = x_{0}^{2}$. Здесь снова три монома от переменных ${{x}_{0}}$ и ${{x}_{1}}$. Это $x_{0}^{2}$, ${{x}_{0}}{{x}_{1}}$ и $x_{1}^{2}$. Но матрица $A$ содержит три столбца и три строки:

Расширенная матрица равна

Матрица $A$ невырожденная, ранги матриц $A$ и $B$ совпадают. Алгоритм 1 отвергает вход.

Следующая теорема обеспечивает достаточное условие типа $m \geqslant n - \sqrt {2n - o(n)} $, при котором алгоритм 1 отвергает большую долю входов при фиксированных значениях $n$ и $m$. Доказательство основано на построении базиса специального вида в пространстве квадратичных форм.

Теорема 2. Даны положительные целые числа $n$ и $m < n$, для которых выполнено неравенство $2n \geqslant (n - m + 1)(n - m + 2)$. Для почти каждого набора из $m$ линейных форм ${{\ell }_{j}}({{x}_{0}}, \ldots ,{{x}_{{n - m}}})$, где $n - m < j \leqslant n$, найдутся значения коэффициентов ${{\lambda }_{1}}$, …, ${{\lambda }_{n}}$, для которых выполнено равенство квадратичных форм

Если для некоторого $\varepsilon > 0$ коэффициенты линейных форм ${{\ell }_{j}}$ независимо и равномерно распределены на множестве из чисел, то вероятность отсутствия искомых значений ${{\lambda }_{1}}$, …, ${{\lambda }_{n}}$ не превышает $\varepsilon $.

Доказательство. Рассмотрим матрицу $A$ из алгоритма 1. По условию теоремы, в матрице $A$ число строк $r = (n - m + 1)(n - m + 2){\text{/}}2$ не превышает числа столбцов $n$. В этом случае достаточным условием существования искомых ${{\lambda }_{1}}$, …, ${{\lambda }_{n}}$ служит равенство ${\text{rank}}(A) = r$. Каждый элемент матрицы $A$ равен либо константе, либо многочлену второй степени от коэффициентов линейных форм ${{\ell }_{j}}$. Следовательно, минор порядка $r$ матрицы $A$ равен многочлену от коэффициентов линейных форм ${{\ell }_{j}}$, степень которого не превышает $2r$.

По лемме Шварца–Зиппеля [25], если этот многочлен не равен нулю тождественно, то вероятность обращения в нуль не превышает .

Осталось показать, что для некоторого набора форм ${{\ell }_{j}}$ угловой минор порядка $r$ матрицы $A$ отличен от нуля. При этом достаточно рассмотреть набор форм, коэффициенты которых не обязательно принадлежат указанному в условии множеству. Положим ${{\ell }_{{n - m + k}}} = {{x}_{0}} + {{x}_{k}}$ для $1 \leqslant k \leqslant n$ – m. Обозначим через $f(i,k)$ номер пары индексов $1 \leqslant i < k \leqslant n - m$, принимающий значения от 1 до $(n - m)(n - m - 1){\text{/}}2$. Положим ${{\ell }_{j}} = {{x}_{i}} + {{x}_{k}}$ для $j = 2(n - m) + f(i,k)$. При $j = r$ полагаем ${{\ell }_{r}}\, = \,2{{x}_{0}}$.

Для индексов $1 \leqslant i < k \leqslant n - m$ моному ${{x}_{i}}{{x}_{k}}$ соответствует строка, где отличен от нуля лишь элемент в столбце $j = 2(n - m) + f(i,k)$. Этот элемент равен двум. С другой стороны, в столбцах с номерами $1 \leqslant j \leqslant n - m$ по два ненулевых элемента: $ - 1$ в строке, соответствующей моному ${{x}_{0}}{{x}_{j}}$, и $1$ в строке, соответствующей моному $x_{j}^{2}$. Переставим строки в матрице $A$ в соответствии с таким порядком мономов: ${{x}_{0}}{{x}_{1}}$, …, ${{x}_{0}}{{x}_{{n - m}}}$, $x_{1}^{2}$, $x_{2}^{2}$, …, $x_{{n - m}}^{2}$, мономы ${{x}_{i}}{{x}_{k}}$ для $1 \leqslant i < k \leqslant n - m$ и последний $x_{0}^{2}$. Угловая подматрица порядка $r$ примет вид

где через $I$ обозначена единичная матрица порядка $n - m$, через $I{\kern 1pt} '$ – единичная матрица порядка $(n - m)(n - m - 1){\text{/}}2$, а через 0 – нулевая матрица. Эта подматрица невырожденная, поскольку элементарные преобразования строк дают треугольную матрицу с ненулевыми элементами на главной диагонали

Следовательно, ${\text{rank}}(A) = r$. ◻

Заметим, что в доказательстве теоремы 2 используется семейство разреженных матриц. Поэтому разреженность матрицы $A$ не служит препятствием для успешного применения обсуждаемого алгоритма. Более того, для разреженных матриц меньше сложность вычисления ранга, что позволяет работать в больших размерностях. В доказательстве теоремы 2 рассмотрен частный случай, когда матрица $A$ имеет полный ранг. Поэтому вычисление ранга можно заменить на проверку невырожденности матрицы $A$. Однако это было бы неоправданным ограничением применимости алгоритма 1.

3. РЕАЛИЗАЦИЯ

Обозначим через $L$ матрицу, составленную из коэффициентов линейных форм ${{\ell }_{j}}$ в алгоритме 1 и теореме 2. Первая строка матрицы $L$ соответствует форме ${{\ell }_{{n - m + 1}}}$, а последняя строка – форме ${{\ell }_{n}}$. Например, для двух линейных форм ${{\ell }_{2}}\, = \,2{{x}_{0}}\, + \,{{x}_{1}}$ и ${{\ell }_{3}} = 3{{x}_{1}}$ эта матрица равна

Сопоставление матрице $L$ матрицы $B$ требует меньше времени, чем вычисление рангов матриц $A$ и $B$. Такое преобразование реализует функция compose_b(), которая получает на вход матрицу $L$ и возвращает матрицу $B$, где обе матрицы представлены массивами NumPy. Порядок строк в $B$ соответствует лексикографическому мономиальному упорядочению, а первая строка матрицы – моному $x_{0}^{2}$. Матрица $A$ легко получается удалением последнего столбца из матрицы $B$. Функция compose_b() реализована на языке Python с использованием библиотеки NumPy.

При использовании целых чисел фиксированной битовой длины, NumPy позволяет значительно повысить эффективность работы с матрицами (как по скорости, так и по памяти) по сравнению со встроенными типами данных языка Python, например списками чисел [26]. При работе с другими типами данных, включая целые неограниченной длины и рациональные числа, эти преимущества в основном утрачиваются. Однако в любом случае NumPy позволяет кратко записывать стандартные операции над матрицами и их частями, например, как в нашем случае, выделение подматрицы, поэлементное сложение и умножение строк и т. п.

Через $s = n - m$ обозначена размерность подпространства, для которого проверяется инцидентность вершине $n$-мерного куба. Через height обозначено число мономов второй степени от переменных ${{x}_{0}}$, …, ${{x}_{s}}$, равное $(s + 1)(s + 2){\text{/}}2$.

Для оценки времени работы функции compose_b() при n = height генерировалась $(n - s) \times (s + 1)$ матрица $L$, элементами которой служат случайно выбранные целые числа из отрезка $[ - {{10}^{9}}{{,10}^{9}}]$. Время, затраченное на вычисление матрицы $B$, показано в таблице 1. В третьем столбце таблицы 1 указаны результаты для 64-битовых целых чисел, а в последнем столбце – для целых чисел неограниченной битовой длины, вычисления над которыми реализованы в языке Python. Мы использовали версии Python 3.10.4 и NumPy 1.22.4. Вычисления проведены на персональном компьютере на базе процессора Intel^® Core i5-3570 и с оперативной памятью 16 гигабайтов. Поэтому для больших матриц вычисления с числами неограниченной битовой длины не были проведены из-за переполнения оперативной памяти, хотя для 64-битовых чисел вычисления продолжены вплоть до значений $s = 300$.

Одновременное вычисление рангов матриц $A$ и $B$ по данной на вход матрице $B$ также реализовано на языке Python. Для вычислений с рациональными числами используется класс Fraction из стандартной библиотеки Python.

Вычисление ранга основано на приведении матрицы к ступенчатому виду. Поскольку матрица $A$ служит подматрицей в $B$, для обеих матриц $A$ и $B$ такое преобразование можно сделать одновременно. Это позволяет примерно вдвое уменьшить время работы по сравнению с независимым вычислением рангов двух матриц.

Чтобы оценить время для вычисления функции rank_ab(), при n = height генерировались $(n - s) \times (s + 1)$ матрицы $L$, элементами которых служат случайно выбранные целые числа из отрезка $[ - N,N]$. Время вычисления для разных значений границы диапазона $N \in {{\{ 10}^{3}}{{,10}^{6}}{{,10}^{9}}\} $ показано в таблице 2. Чем больше размер матрицы, тем быстрее возрастает время работы при увеличении средней битовой длины элементов матрицы.

Отметим, что реализованный в библиотеке NumPy метод вычисления ранга, основанный на сингулярном разложении матрицы, в этой задаче применять небезопасно из-за возможных ошибок, связанных с использованием чисел с плавающей запятой. Например, при значениях $k \geqslant 16$ вычисляемый в NumPy 1.22.4 таким неточным способом ранг диагональной 2 × 2 матрицы с элементами 1 и 10^k на главной диагонали оказывается равным единице. Поэтому для вычисления ранга мы использовали новую реализацию известного алгоритма Гаусса приведения матрицы к ступенчатому виду, которая позволяет работать с числами неограниченной битовой длины.

В итоге реализацией алгоритма 1 служит функция may_have_01solution(), которая получает на вход матрицу $L$. Значение равно True, когда алгоритм 1 дал бы неопределенный ответ. Значение False означает, что алгоритм 1 отвергает вход. Входом может служить матрица с рациональными элементами.

Листинги обсуждаемых функций и примеры их использования доступны по адресу http://lab6.iitp.ru/-/qualg

4. ЗАКЛЮЧЕНИЕ

Реализован на языке Python эвристический алгоритм полиномиального времени для распознавания подпространств достаточно малой размерности, инцидентных какой-либо из вершин единичного многомерного куба, хотя в худшем случае эта задача считается вычислительно трудной. Алгоритм либо корректно отвергает вход, либо сообщает об отсутствии очевидного препятствия к инцидентности. В частности, полученные ограничения на сложность этой задачи могут быть полезны для оценки надежности криптографических протоколов, основанных на поиске (0, 1)-решения системы уравнений [19, 27].

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

Авторы благодарны анонимному рецензенту за сделанные замечания, позволившие улучшить статью.