Алгоритм Видемана

Алгоритм Видемана — алгоритм, позволяющий получить решение системы линейных уравнений A x = b , b ≠ 0 {displaystyle Ax=b,b eq 0} над конечным полем K = G F ( q ) {displaystyle K=GF(q)} . Был предложен Дугласом Видеманом (англ. Douglas Wiedemann) в 1986 году. В течение некоторого времени после опубликования статьи, алгоритм не получил большой поддержки и считался пригодным только для получения наилучших оценок сложности. Но позже алгоритмы Видемана были реализованы на компьютере и использовались, например, для поиска разложения многочленов на множители над конечными полями.

История возникновения

Алгоритмы Видемана были представлены общественности в прошлом столетии. В 1986 году в январском выпуске журнала IEEE Transactions on Information Theory была опубликована статья Дугласа Видемана под названием «Решение разреженной системы линейных уравнений над конечным полем» (англ. Solving sparse linear equations over finite fields). В ней были описаны алгоритмы для решения системы линейных уравнений над конечным полем в случае когда матрица системы является разреженной. Причём в статье были рассмотрены случаи с различными разреженными матрицами. Также в статье было опубликовано обоснование алгоритмов и оценка сложности их работы.

Задача

Алгоритм нужен чтобы решить систему линейных уравнений A x = b , b ≠ 0 {displaystyle Ax=b,b eq 0} . Матрица A {displaystyle A} имеет размерность n × n {displaystyle n imes n} и предполагается разреженной, количество ненулевых элементов в ней равно w {displaystyle w} .

Теория

С помощью матрицы A {displaystyle A} определяется невырожденное линейное отображение(которое обозначается также A {displaystyle A} ) на пространстве K n {displaystyle mathrm {K} ^{n}} . Рассматривается пространство S {displaystyle S} , порождённое множеством векторов { ( A i b k ) } i = 0 , 1 , 2 {displaystyle left{{(A^{i}b_{k})} ight}_{i=0,1,2}} и определяется A S = A | S {displaystyle A_{S}=A|_{S}} - линейное отображение S {displaystyle S} на S {displaystyle S} .

Обозначим f z ∈ K n {displaystyle f_{z}in mathrm {K} ^{n}} — минимальный многочлен A S {displaystyle A_{S}} , то есть ненулевой многочлен наименьшей степени, такой, что f ( A s ) {displaystyle f(A_{s})} является нулевым отображением S {displaystyle S} , при чём нормализованный так, что его свободный член равен единице. Отметим, что если g ( z ) ∈ K [ z ] {displaystyle g(z)in mathrm {K} [z]} , то g ( A s ) {displaystyle g(A_{s})} - нулевое отображение тогда и только тогда, когда g ( A ) b = 0 {displaystyle g(A)b=0} . Кроме того, f ( z ) {displaystyle f(z)} делит многочлен d e t ( z l n − A ) {displaystyle det(zl_{n}-A)} , и поэтому d e g f ( z ) ⩽ n {displaystyle degf(z)leqslant n} .

Обозначим d = d e g f ( z ) , f ( z ) = ∑ i = 0 d f [ i ] z i {displaystyle d=degf(z),f(z)=sum _{i=0}^{d}{f[i]z^{i}}} , где f [ i ] ∈ K n {displaystyle f[i]in mathrm {K} ^{n}} - коэффициенты f ( z ) {displaystyle f(z)} . Если можно найти f ( z ) {displaystyle f(z)} , то решение системы A x = b , b ≠ 0 {displaystyle Ax=b,b eq 0} также находится: так как f ( A ) b = 0 {displaystyle f(A)b=0} и f [ 0 ] = 1 {displaystyle f[0]=1} , то

x = − ∑ i = 1 d f [ i ] A i − 1 b {displaystyle x=-sum _{i=1}^{d}{f[i]A^{i-1}b}}

Пусть u {displaystyle u} - какой-либо фиксированный вектор из K n {displaystyle mathrm {K} ^{n}} . Обозначим стандартное билинейное отображение K n {displaystyle mathrm {K} ^{n}} в K {displaystyle mathrm {K} } как ( , ) {displaystyle (,)} , то есть ( ( v 1 , . . . , v n ) , ( w 1 , . . . , w n ) ) = ∑ i = 1 n v i w i {displaystyle ((v_{1},...,v_{n}),(w_{1},...,w_{n}))=sum _{i=1}^{n}{v_{i}w_{i}}} .

Так как f ( A ) b = 0 {displaystyle f(A)b=0} , то последовательность

( u , A i b ) , i = 0 , 1 , 2 , . . . , {displaystyle (u,A^{i}b),i=0,1,2,...,}

удовлетворяет линейному рекуррентному соотношению, характеристический многочлен которого равен f ( z ) {displaystyle f(z)} . Пусть f u ( z ) {displaystyle f_{u}(z)} - характеристический многочлен самого короткого рекуррентного соотношения. Тогда f ( z ) | f u ( z ) {displaystyle f(z)|f_{u}(z)} . Действительно, если разделить с остатком

f ( z ) = q ( z ) f u ( z ) + r ( z ) {displaystyle f(z)=q(z)f_{u}(z)+r(z)} , d e g r ( z ) < d e g   f u ( z ) {displaystyle degr(z)<deg~f_{u}(z)}

то из равенств

0 = ( u , f ( A ) b ) = ( u , q ( A ) f u ( A ) b ) + ( u , r ( A ) b ) {displaystyle 0=(u,f(A)b)=(u,q(A)f_{u}(A)b)+(u,r(A)b)} ,

( u , f u ( A ) A j b ) = 0 , j = 0 , 1 , 2 , . . . , {displaystyle (u,f_{u}(A)A^{j}b)=0,j=0,1,2,...,}

и минимальности f u ( z ) {displaystyle f_{u}(z)} будет следовать, что r ( z ) = 0 {displaystyle r(z)=0} . Поскольку свободный член f ( z ) {displaystyle f(z)} равен единице, то можно принять, что свободный член f u ( z ) {displaystyle f_{u}(z)} равен единице.

Минимальный многочлен f u ( z ) {displaystyle f_{u}(z)} для последовательности ( u , A i b ) , i = 0 , 1 , 2 , . . . , {displaystyle (u,A^{i}b),i=0,1,2,...,} может быть получен с помощью алгоритма Берлекэмпа-Месси по первым её 2 n {displaystyle 2n} членам. Существуют два метода решения исходной системы.

Первый метод. Выбирается случайный вектор u ( z ) {displaystyle u(z)} . Строится f u ( z ) {displaystyle f_{u}(z)} и в предположении, что f ( z ) = f u ( z ) {displaystyle f(z)=f_{u}(z)} , находится x {displaystyle x} по формуле

x = − ∑ i = 1 d f [ i ] A i − 1 b {displaystyle x=-sum _{i=1}^{d}{f[i]A^{i-1}b}}

Этим путём с достаточно высокой вероятностью можно найти решение.

Второй метод. Пусть b 0 = b , f 1 ( z ) = f u 1 ( z ) {displaystyle b_{0}=b,f_{1}(z)=f_{u_{1}}(z)} для некоторого вектора u 1 {displaystyle u_{1}} . Если вектор b 1 = f 1 ( A ) b 0 {displaystyle b_{1}=f_{1}(A)b_{0}} равен 0, то находится x {displaystyle x} по формуле

x = − ∑ i = 1 d f [ i ] A i − 1 b {displaystyle x=-sum _{i=1}^{d}{f[i]A^{i-1}b}} (так как тогда f 1 ( z ) = f ( z ) {displaystyle f_{1}(z)=f(z)} ).

Если же b 1 ≠ 0 {displaystyle b_{1} eq 0} , то повторяется процедура, то есть выбирается случайный вектор u 2 {displaystyle u_{2}} и строится минимальный многочлен f 2 ( z ) = f u 2 ( z ) {displaystyle f_{2}(z)=f_{u_{2}}(z)} для последовательности ( u 2 , A i b 1 ) {displaystyle (u_{2},A_{i}b_{1})} . Если b 2 = f 2 ( A ) b 1 = 0 {displaystyle b_{2}=f_{2}(A)b_{1}=0} , то f ( z ) = f 1 ( z ) f 2 ( z ) {displaystyle f(z)=f_{1}(z)f_{2}(z)} и можно найти решение x по формуле

x = − ∑ i = 1 d f [ i ] A i − 1 b {displaystyle x=-sum _{i=1}^{d}{f[i]A^{i-1}b}} ,

иначе выбирается u 3 {displaystyle u_{3}} и так далее.

Докажем, что если сделано k {displaystyle k} итераций, то f 1 ( z ) . . . f k ( z ) {displaystyle f_{1}(z)...f_{k}(z)} делит f ( z ) {displaystyle f(z)} . Выше было показано, что f − 1 ( z ) | f ( z ) {displaystyle f-1(z)|f(z)} . Далее, если предположить что f 1 ( z ) . . . f k − 1 ( z ) {displaystyle f_{1}(z)...f_{k-1}(z)} делит f ( z ) {displaystyle f(z)} , то поскольку f k ( z ) {displaystyle f_{k}(z)} - минимальный многочлен для последовательности { ( u k , A i b k − 1 ) } i = { ( u k , f k − 1 ( A ) . . . f 1 ( A j ) b ) } j {displaystyle left{(u_{k},A^{i}b_{k-1}) ight}_{i}=left{(u_{k},f_{k-1}(A)...f_{1}(A^{j})b) ight}_{j}} , а многочлен f ( z ) f 1 ( z ) . . . f k − 1 ( z ) {displaystyle {frac {f(z)}{f_{1}(z)...f_{k-1}(z)}}} её аннулирует, то f k ( z ) | f ( z ) f 1 ( z ) . . . f k − 1 ( z ) {displaystyle {f_{k}(z)}|{frac {f(z)}{f_{1}(z)...f_{k-1}(z)}}} , что и требовалось доказать.

Теперь очевидно, что если b x = f k ( A ) . . . f 1 ( A ) b = 0 {displaystyle b_{x}=f_{k}(A)...f_{1}(A)b=0} , то f ( x ) = f 1 ( x ) . . . f k ( x ) {displaystyle f(x)=f_{1}(x)...f_{k}(x)} . То есть, как только будет найден нулевой вектор b k = f k ( A ) b k − 1 {displaystyle b_{k}=f_{k}(A)b_{k-1}} , то можно найти решение исходной системы по формуле

x = − ∑ i = 1 d f [ i ] A i − 1 b {displaystyle x=-sum _{i=1}^{d}{f[i]A^{i-1}b}} .

Алгоритм 1

В оригинальной статье алгоритм имеет такое название. На его основе строится детерминированный алгоритм, который в оригинальной статье называется алгоритм 2.

Описание алгоритма

1 этап. Приравнивается b 0 = b , k = 0 , y 0 = 0 , d 0 = 0 {displaystyle b_{0}=b,k=0,y_{0}=0,d_{0}=0} .

2 этап. Если b k = 0 {displaystyle b_{k}=0} , то решение равно x = − y k {displaystyle x=-y_{k}} , и алгоритм прекращает работу.

3 этап. Выбирается случайный вектор u k + 1 ∈ K n , u k + 1 ≠ 0 {displaystyle u_{k+1}in mathrm {K} ^{n},u_{k+1} eq 0} .

4 этап. Вычислить первые 2 ( n − d k ) {displaystyle 2(n-d_{k})} членов последовательности { ( u k + 1 , A i b k ) } i = 0 ∞ {displaystyle left{{(u_{k+1},A^{i}b_{k})} ight}_{i=0}^{infty }} .

5 этап. Вычислить минимальный многочлен f k + 1 ( z ) {displaystyle f_{k+1}(z)} последовательности из 4-го этапа, причём нормализовать его так, чтобы его свободный член равнялся единице. Это можно осуществить с помощью алгоритма Берлекэмпа-Месси.

6 этап. Присвоить

y k + 1 = y k + f ^ k + 1 ( A ) b k {displaystyle y_{k+1}=y_{k}+{hat {f}}_{k+1}(A)b_{k}} , где f ^ ( z ) = f ( z ) − f ( 0 ) z {displaystyle {hat {f}}(z)={frac {f(z)-f(0)}{z}}}

b k + 1 = b 0 + A y k + 1 {displaystyle b_{k+1}=b_{0}+Ay_{k+1}} ,

d k + 1 = d k + d e g   f k + 1 ( z ) {displaystyle d_{k+1}=d_{k}+deg~f_{k+1}(z)} .

7 этап. Присвоить k = k + 1 {displaystyle k=k+1} и вернуться на второй этап.

Обоснование корректности алгоритма с помощью метода математической индукции

f ^ ( z ) = f ( z ) − f ( 0 ) z {displaystyle {hat {f}}(z)={frac {f(z)-f(0)}{z}}} соответствует правой части формулы x = − ∑ i = 1 d f [ i ] A i − 1 b {displaystyle x=-sum _{i=1}^{d}{f[i]A^{i-1}b}} без знака минус. При k = 0 {displaystyle k=0} выбирается u 1 {displaystyle u_{1}} , рассматривается 2 n {displaystyle 2n} членов последовательности { ( u 1 , A i b 0 ) } i = 0 , 1 , 2 {displaystyle left{{(u_{1},A^{i}b_{0})} ight}_{i=0,1,2}} и находится f 1 ( x ) {displaystyle f_{1}(x)} по алгоритму Берлекэмпа-Месси. Тогда y 1 = f ^ 1 ( A ) b , b 1 = b 0 + A y 1 = b + A f 1 ( A ) − 1 A b = f 1 ( A ) b , d 1 = d e g f 1 ( z ) {displaystyle y_{1}={hat {f}}_{1}(A)b,b_{1}=b_{0}+Ay_{1}=b+A{frac {f_{1}(A)-1}{A}}b=f_{1}(A)b,d_{1}=degf_{1}(z)} .

Пусть после k {displaystyle k} проходов алгоритма выполнены равенства

y k = f k ( A ) . . . f 1 ( A ) − 1 A b {displaystyle y_{k}={frac {f_{k}(A)...f_{1}(A)-1}{A}}b}

b k = f k ( A ) . . . f 1 ( A ) b {displaystyle b_{k}=f_{k}(A)...f_{1}(A)b}

Тогда после k + 1 {displaystyle k+1} прохода

y k + 1 = y k + f ^ k + 1 ( A ) b k = f k ( A ) . . . f 1 ( A ) − 1 A b + f k + 1 ( A ) − 1 A f k ( A ) . . . f 1 ( A ) b = f k + 1 ( A ) . . . f 1 ( A ) − 1 A b {displaystyle y_{k+1}=y_{k}+{hat {f}}_{k+1}(A)b_{k}={frac {f_{k}(A)...f_{1}(A)-1}{A}}b+{frac {f_{k+1}(A)-1}{A}}f_{k}(A)...f_{1}(A)b={frac {f_{k+1}(A)...f_{1}(A)-1}{A}}b} ,

b k + 1 = b k + A f k + 1 ( A ) . . . f 1 ( A ) − 1 A b = f k + 1 ( A ) . . . f 1 ( A ) b {displaystyle b_{k+1}=b_{k}+A{frac {f_{k+1}(A)...f_{1}(A)-1}{A}}b=f_{k+1}(A)...f_{1}(A)b}

То есть формулы для y k {displaystyle y_{k}} и b k {displaystyle b_{k}} сохраняются. Теперь корректность алгоритма следует из раздела теория.

Детерминированный алгоритм

Описание алгоритма

1 этап. Найти значение A i b , i = 0 , 1 , . . . , 2 n − 1 {displaystyle A^{i}b,i=0,1,...,2n-1} .

2 этап. Приравнять нулю k {displaystyle k} , а g 0 ( z ) {displaystyle g_{0}(z)} единице.

3 этап. Присвоить u k + 1 = ( 0 , . . . , 0 , 1 , 0 , . . . , 0 ) {displaystyle u_{k+1}=(0,...,0,1,0,...,0)} (единица находится на k + 1 {displaystyle k+1} месте).

4 этап. Найти последовательность ( u k + 1 , A i b ) , i = 0 , 1 , . . . , 2 n − 1 {displaystyle (u_{k+1},A^{i}b),i=0,1,...,2n-1} при помощи первого этапа.

5 этап. Найти последовательность ( u k + 1 , g k ( A ) A i b ) , i = 0 , 1 , . . . , 2 n − 1 − d e g   g k ( z ) {displaystyle (u_{k+1},g_{k}(A)A^{i}b),i=0,1,...,2n-1-deg~g_{k}(z)} , можно использовать дискретное преобразование Фурье.

6 этап. Найти минимальный многочлен f k + 1 ( z ) {displaystyle f_{k+1}(z)} со свободным членом равным единице с помощью алгоритма Берлекэмпа-Месси.

7 этап. Присвоить g k + 1 ( z ) = f k + 1 ( z ) g k ( z ) {displaystyle g_{k+1}(z)=f_{k+1}(z)g_{k}(z)} .

8 этап. Увеличить k {displaystyle k} на единицу. Если d e g   g k ( z ) < n {displaystyle deg~g_{k}(z)<n} и k < n {displaystyle k<n} , то возвратиться на 3 этап.

9 этап. Для многочлена f ( z ) = g k ( z ) {displaystyle f(z)=g_{k}(z)} с помощью найденных на первом этапе значений A i b {displaystyle A^{i}b} отыскать решение x {displaystyle x} системы с помощью формулы

x = − ∑ i = 1 d f [ i ] A i − 1 b {displaystyle x=-sum _{i=1}^{d}{f[i]A^{i-1}b}} .

Обоснование корректности алгоритма

Обратим внимание, что фактически алгоритм работает также, как и алгоритм 1, только векторы u k {displaystyle u_{k}} выбираются не случайно, а идёт перебор единичных векторов ( 0 , . . . , 0 , 1 , 0 , . . . , 0 ) {displaystyle (0,...,0,1,0,...,0)} . Очевидно, что g k ( z ) = f k ( z ) . . . f 1 ( z ) {displaystyle g_{k}(z)=f_{k}(z)...f_{1}(z)} , где f k ( z ) {displaystyle f_{k}(z)} — минимальный многочлен для последовательности

( u k , f k − 1 ( A ) . . . f 1 ( A ) A i b ) , i = 0 , 1 , . . . , 2 n − 1 − d e g ( f k − 1 . . . f 1 ) {displaystyle (u_{k},f_{k-1}(A)...f_{1}(A)A^{i}b),i=0,1,...,2n-1-deg(f_{k-1}...f_{1})} .

Алгоритм закончил работу при некотором значении параметра k {displaystyle k} . Предположим, что k < n , d e g   g k ( z ) = n {displaystyle k<n,deg~g_{k}(z)=n} . Так как d e g   f ( z ) ⩽ n {displaystyle deg~f(z)leqslant n} и g k ( z ) | f ( z ) {displaystyle g_{k}(z)|f(z)} , то g k ( z ) = f ( z ) {displaystyle g_{k}(z)=f(z)} . Отсюда следует, что на этапе 9 решение исходной системы точно будет найдено.

Теперь рассмотрим случай k = n {displaystyle k=n} . Поскольку был совершён перебор всех единичных векторов u 1 , . . . , u n {displaystyle u_{1},...,u_{n}} , то вектор g n ( A ) b {displaystyle g_{n}(A)b} ортогонален u 1 , . . . , u n {displaystyle u_{1},...,u_{n}} . Следовательно, g n ( A ) b = 0 {displaystyle g_{n}(A)b=0} . Так как g n ( z ) | f ( z ) {displaystyle g_{n}(z)|f(z)} и f ( z ) {displaystyle f(z)} -минимальный многочлен, то g k ( z ) = f ( z ) {displaystyle g_{k}(z)=f(z)} . Поэтому и в данном случае подтверждена корректность работы алгоритма.

Оценка сложности алгоритма

Для детерминированного алгоритма Видеманом была получена следующая оценка сложности: O ( n ( w + n   l o g ( n ) l o g ( l o g ( n ) ) ) ) {displaystyle O(n(w+n~log(n)log(log(n))))} . Полученная оценка сложности является наилучшей среди известных. Благодаря алгоритму Видемана возможно улучшение оценки сложности в других алгоритмах, использующих методы решения линейных систем.

Аналогичные алгоритмы

Где может пригодится решение системы линейных уравнений над конечным полем? Потребность в их решении возникает при использовании алгоритмов факторизации и при решении задач дискретного логарифмирования, использующих факторные базы. Существует большое количество алгоритмов для получения решения системы линейных уравнений над конечными полями. Помимо алгоритмов Видемана можно использовать гауссово и структурное гауссово исключения, алгоритм Ланцоша, метод сопряжённых градиентов . Также известны алгоритмы основанные на быстром умножении матриц, например на алгоритмах Штрассена и Копперсмита-Винограда. Свои алгоритмы были предложены Коновальцевым и Бриллхартом.

В общем случае (матрица системы не является разреженной) в последнее время чаще используется алгоритм Ланцоша(вероятно, вместе со структурированным гауссовым исключением для получения более плотной матрицы подсистемы). Но в случае разреженной матрицы эффективнее всего использовать алгоритмы Видемана, так как оценки их сложности являются наилучшими из известных. Не сразу алгоритмы Видемана получили признание, но позже всё-таки были реализованы на компьютере. Алгоритмы использовались, например, для разложения многочленов на множители над конечными полями.

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