Как мне вычислить нулевое пространство прямоугольной разреженной матрицы над GF (2) в C/C++?

Обновленный ответ - спасибо Гарольду: QR-разложение вообще не применимо для вашего случая.

См. например

https://math.stackexchange.com/questions/1346664/how-to-find-orthogonal-vectors-in-gf2

Я ошибочно предположил, что QR здесь применим, но это не теоретически.

Если вас все еще интересуют подробности о QR-алгоритмах, откройте новую тему.

Назовем рассматриваемую матрицу М. Тогда (поправьте меня, если я ошибаюсь):

1. GF(2) подразумевает, что M эквивалентна битовой матрице — каждый элемент может иметь одно из двух значений.

2. Арифметика на GF(2) аналогична целочисленной арифметике с неотрицательными числами, но выполняется по модулю 2, поэтому сложение — побитовое XOR, а умножение — побитовое AND. Неважно, какие именно элементы содержит GF(2) — все они эквивалентны битам.

3. Векторы в GF(2) линейно независимы, пока они не равны, или пока они отличаются хотя бы на один бит, или v_1 + v_2 ≠ 0 (поскольку сложение в GF(2) является логическим XOR).

По определению (правое) нулевое пространство охватывает базисные векторы, которые матрица преобразует в 0. Вектор v окажется в нулевом пространстве, если умножить каждый j-й столбец M на j-й бит v, просуммировать их и получить результат нулевой.

Я вижу как минимум два пути решения этой проблемы.

1. Выполните плотное исключение Гаусса с точки зрения битовых операций, организуйте данные и напишите циклы так, чтобы компилятор все векторизовал и работал с 512-битными типами данных. Вы можете использовать Compiler Explorer на godbolt.org, чтобы легко проверить, что векторизация имеет место, и, например. Используются инструкции AVX512. Конечно, линейный выигрыш в конечном итоге проиграет при квадратичном масштабировании проблемы, но прирост производительности по сравнению с наивной реализацией на основе bool будет огромным и может быть достаточным для ваших нужд. Разреженность добавляет возможную сложность: если матрица не помещается в памяти в плотном представлении, то необходимо разработать подходящее представление, которое обеспечит хорошую работу исключения Гаусса. Необходимо больше знать о матрицах, над которыми вы работаете. Вообще говоря, операции со строками будут выполняться с пропускной способностью памяти, если реализация правильная, порядка 1E10 элементов/с, поэтому 1E3x1E3 M должно обрабатываться не более чем за секунду.

2. Поскольку проблема эквивалентна набору логических уравнений, используйте решатель SAT (решатель логических задач выполнимости) для постепенного создания нулевого пространства. Исходным набором уравнений является M × v = 0 и v ≠ 0, где v — битовый вектор. Запускайте SAT, пока не найдёт какой-нибудь v, назовём его v_i. Затем добавьте ограничение v ≠ v_i и снова запустите SAT, добавляя ограничения на каждой итерации. То есть k-я итерация имеет ограничения v ≠ 0, v ≠ v1, ... v ≠ v(k-1).

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

   Современный SAT превосходно справляется с разреженными задачами с большим количеством логических уравнений, чем с переменными, поэтому я полагаю, что это будет работать очень хорошо - чем разреженнее матрица, тем лучше. Задача должна быть предварительно обработана, чтобы удалить все нулевые столбцы в M, чтобы минимизировать комбинаторный взрыв. Решатели SAT с открытым исходным кодом могут легко справляться с задачами с 1 млн переменных, поэтому для разреженной задачи вы могли бы реально решить 100 000–1 млн столбцов в M и около 10 «единиц» в каждой строке. Таким образом, разреженная матрица 1Mx1M с 10 «единицами» в каждой строке в среднем была бы разумной задачей для обычных SAT-решателей, и я полагаю, что современный уровень техники может работать с матрицами 10Mx10M и выше.

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

   Поскольку вы уже используете Eigen, проблема, по крайней мере, будет соответствовать представлению SparseMatrix с элементами размером в байт, так что это не очень большая проблема с точки зрения SAT.

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