РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ В КОММУТАТИВНЫХ ОБЛАСТЯХ ^{∗}

\begin{matrix}c \\ ◯ \end{matrix} Г.И. Малашонок

Malaschonok G.I. Solving systems of linear equations in commutative domains. There are considered fast methods of linear systems solving in commutative domains and in their fields of fractions. There is suggested the realization of these methods for the rings of principal ideals with the help of Hensel lifting techniques. The complexity of the deterministic algorithm to find all the rational solutions for linear systems with coefficient matrix A of size n × m and rank n with p-adic lifting is O(n^{2}m(m − n + 1)(log n + log ∥A∥ + log p)^{2}) bit operations. The randomized algorithm to find all the Diophantine solutions for this system has essentially the same complexity.

∗Вестник ТГУ, Т.5, N.1, 2000, 147–154.

ВВЕДЕНИЕ

Пусть R – коммутативная область с единицей, K – поле частных R. Одной из важных задач вычислительной коммутативной алгебры является задача решения систем линейных уравнений Ax = b, коэффициенты которых берутся из коммутативной области R: A ∈ R^{n×m}, b ∈ R^{n}.

При этом могут ставиться разные задачи: (1^{o}) найти одно решение в K^{m}, (2^{o}) найти все решения в K^{m}, (3^{o}) найти одно решение в R^{m}, (4^{o}) найти все решения R^{m}. Рассмотрим каждую из этих задач для разных областей.

Область R является полем. Основным вычислительным методом здесь является метод исключения Гаусса и его модификации со сложностью O(n^{2}m) операций в поле R. Быстрые методы решения систем линейных уравнений над полем ведут свое начало от работы Штрассена [1] и имеют такую же сложность, как и матричное умножение – O(n^{β−1}m). O(n^{β}) – это число операций, необходимое для умножения квадратных матриц n-го порядка. Для обычного умножения матриц β = 3, для алгоритма Штрассена β = log_{2} 7, а для лучшего на сегодня алгоритма β < 2.376 [2].

Произвольная коммутативная область R с единицей. Оценку алгоритмов будем производить по количеству операций над элементами R. При этом предполагаем, что в R существует алгоритм сокращения, то есть если a = bc и известны a и c, c ≠ 0, то существует алгоритм вычисления b. Основным вычислительным методом здесь является метод Доджсона, иногда называемый методом Гаусса с точным делением, со сложностью O(n^{2}m) операций в R. Впервые такой метод был опубликован в работе Доджсона [3], а его улучшения получены в работах [4], [5], [6]. Быстрый метод решения систем линейных уравнений над коммутативной областью опубликован в работе [7]. Сложность этого метода O(n^{β−1}m) такая же, как и сложность матричного умножения. И теперь дальнейший прогресс в задачах (1^{o}) и (2^{o}) целиком связан с улучшениями алгоритма матричного умножения.

Вероятностный алгоритм для решения задач (3^{o}) и (4^{o}) при самых общих предположениях относительно коммутативной области обсуждается в настоящей работе. Предполагается, что в области R имеется алгоритм, который может установить, является ли конечнопорожденный идеал I = (a_{1}, a_{2}, …, a_{s}) единичным. Если I является единичным, то этот алгоритм вычисляет коэффициенты k_{i} ∈ R в разложении единицы 1 = ∑_{i=1}^{s} k_{i}a_{i}. Например, достаточно иметь алгоритм вычисления в области R базиса Гребнера.

Коммутативная область главных идеалов R. Это основная прикладная область. Лучшим методом для задачи (1^{o}) в случае решения определенных систем в поле рациональных чисел (R = Z, m = n, det A ≠ 0) остается метод Диксона [9], использующий линейный p-адический подъем. Его сложность O(n^{3}(log n + log ∥A∥ + log p)^{2}) бит-операций.

Если применять метод Гаусса с точным делений для R = Z, то для решения задачи (2^{o}) в обычной арифметике потребуется O(mn^{4}(log n+log ∥A∥)^{2}) бит-операций. ∥A∥ обозначает модуль наибольшего из коэффициентов системы. Если применить китайскую теорему об остатках, то сложность можно уменьшить до величины O(mn^{3}(log n+log ∥A∥)^{2}) битопераций [8].

Мой подход к решению задачи (2^{o}) в области главных идеалов состоит в применении pадического подьема, как и в методе Диксона [9]. Предлагаемый в настоящей работе метод имеет сложность O(n^{2}m(m − n + 1)(log m + log ∥A∥ + log p)^{2}) бит-операций.

В работе [10] предлагается вероятностный алгоритм для задачи (3^{o}) – получения одного решения в кольце целых чисел Z и кольце F [x] многочленов над полем F, который опирается на алгоритм Диксона.

Я предлагаю вероятностный алгоритм для задачи (4^{o}), который основан на p-адическом подьеме и сложность которого имеет существенно кубическую зависимость от порядка системы, как и [10], но количество вычислений обратных матриц становится в m − n + 1 раз меньше.

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

РЕШЕНИЕ СИСТЕМЫ В ПОЛЕ ЧАСТНЫХ

Пусть R – коммутативная область с единицей, K – поле частных области R, A ∈ R^{n×m}, c ∈ R^{n},

Ax = c  (1)

– система линейных уравнений над R. Пусть rankA = r, S и T – матрицы перестановок, такие, что det A_{0} ≠ 0, где A_{0} верхний левый блок

размера r × r матрицы SAT = ( \begin{matrix}A_{0} & A_{1} \\ A_{2} & A_{3} \end{matrix} )

и пусть Sc = ( \begin{matrix}c_{0} \\ c_{1} \end{matrix} ), c_{0} ∈ R^{r}. M = \{x|x ∈

K^{m}, Ax = c} – множество решений системы (1). Единичную матрицу порядка r будем обозначать I_{r}, а также обозначим E_{i,j} – квадратные матрицы, у которых только один ненулевой элемент – (i, j), и этот элемент равен единице. Будем пользоваться этими обозначениями на протяжении всей статьи.

Нам потребуются некоторые факты из теории систем линейных уравнений.

1. Если rank(A, c) ≠ r, то M = ∅. Если rank(A, c) = r, то M - это гиперплоскость размерности m − r в пространстве размерности m. Она определяется m − r + 1 точками, которые не лежат в одной гиперплоскости меньшей размерности.

2. Если система (1) однородная (c = 0) и x_{1}, x_{2}, …, x_{m−r} ее линейно независимые решения, то M = {∑_{i=1}^{m−r+1} x_{i}u_{i}|u_{i} ∈ K}.

3. Если система (1) неоднородная (c ≠ 0) и x_{1}, x_{2}, …, x_{m−r+1} ее линейно независимые решения, то M = {∑_{i=1}^{m−r+1} x_{i}u_{i}|u_{i} ∈ K, ∑_{i=1}^{m−r+1} u_{i} = 1}.

Определение. Базисным множеством решений однородной системы линейных уравнений (1) называется множество, состоящее из m − r линейно независимых решений системы (1). Базисным множеством решений неоднородной системы линейных уравнений (1) называется множество, состоящее из m−r+1 линейно независимых решений системы (1).

В следующих двух теоремах задача вычисления базисного множества решений системы (1) сводится к нескольким задачам решения определенных систем. В первой теореме рассматриваются однородные системы, во второй теореме – неоднородные системы.

Теорема 1 Пусть (1) – однородная система линейных уравнений и A_{1} = (a_{1}, a_{2}, …, a_{m−r}), a_{j} ∈ R^{r}. Тогда системы

A_{0}x_{j} = −a_{j}, j = 1, …, m − r  (2)

определены, и их решения x_{j} ∈ K^{r} задают базисное множество решений системы (1):

T ( \begin{matrix}x_{j} \\ e_{j} \end{matrix} ), j = 1, …, m − r,  (3)

где e_{j} ∈ R^{m−r} — это столбцы единичной матрицы I_{m−r} = (e_{1}, e_{2}, …, e_{m−r}).

Доказательство. Обозначим y = T^{−1}x. По условию теоремы (A_{0}, A_{1})y = 0. Если искать ре-

шение в виде y = ( \begin{matrix}x_{j} \\ e_{j} \end{matrix} ), то придем к системе

(2).

Следствие 1 Пусть B = A_{0}^{−1}A_{1} = (b_{1}, b_{2},

…, b_{m−r}), b_{j} ∈ K^{r}, тогда T ( \begin{matrix}−b_{1} \\ e_{1} \end{matrix} ),

T ( \begin{matrix}−b_{2} \\ e_{2} \end{matrix} ), …, T ( \begin{matrix}−b_{m−r} \\ e_{m−r} \end{matrix} ) базисное множе-

ство решений однородной системы линейных уравнений (1).

Доказательство. В (3) подставим решения систем (2).

Теорема 2 Пусть (1) – неоднородная система линейных уравнений, P – матрица перестановок, такая, что последний элемент вектора b = P A_{0}^{−1}c_{0} не равен 0. Пусть B = P A_{0}^{−1}A_{1}, J ⊂ {1, …, m−r} – номера столбцов матрицы B с нулевыми элементами в последней строке и U = I_{m−r+1} + ∑_{j∈J} E_{1,j+1}. Пусть P A_{0}P = (A′_{0}, a_{0}), (a_{0}, P A_{1})U = (a_{0}, a_{1}, a_{2}, …, a_{m−r}). Тогда системы

(A′_{0}, a_{j})x_{j} = P c_{0}, j = 0, 1, …, m − r, (4)

определены и их решения x_{j} = ( \begin{matrix}x′ & j \\ ξ_{j} \end{matrix} ), x′_{j} ∈

K^{r−1}, ξ_{j} ∈ K задают базисное множество решений системы (1):

T V ( \begin{matrix}x′ & j \\ ξ_{j}e′_{j} \end{matrix} ) j = 0, 1, …, m − r,  (5)

где e′_{j} ∈ R^{m−r+1} — это столбцы единичной матрицы I_{m−r+1} = (e′_{0}, e′_{1}, …, e′_{m−r}), Q = diag(P, I_{m−r}), W = diag(I_{r−1}, U), V = QW.

Доказательство. Обозначим y = V^{−1}T^{−1}x. По условию теоремы P (A_{0}, A_{1})V y = P c_{0} и P (A_{0}, A_{1})V = (A′_{0}, a_{0}, a_{1}, …, a_{m−r}).

Если искать решение в виде y = ( \begin{matrix}x′ & j \\ ξ_{j}e′_{j} \end{matrix} ), то

придем к системам (4).

Покажем, что системы (4) определены и ξ_{j} ≠ 0. Для этого умножим их слева на P A_{0}^{−1}P. Так как P = P^{−1}, то получим P A_{0}^{−1}P (A′_{0}, a_{j})x_{j} = P A_{0}^{−1}c_{0}, j = 0, 1, …, m − r.

Так как b = P A_{0}^{−1}c_{0} и P A_{0}^{−1}P (A′_{0}, a_{0}) = I_{r}, то первая из систем (4) принимает вид I_{r}x_{0} = b. Обозначим I_{r} = (I′, e), e ∈ R^{r}, тогда P A_{0}^{−1}P A′_{0} = I′, P A_{0}^{−1}P a_{0} = e.

Так как по определению P A_{0}^{−1}P (a_{0}, …, a_{m−r}) = P A_{0}^{−1}P (a_{0}, P A_{1})U = (e, B)U, то обозначим d_{j} = P A_{0}^{−1}P a_{j}, j = 1, …, m − r столбцы матрицы (e, B)U. Так как U = I_{m−r+1} + ∑_{j∈J} E_{1,j+1} и J – номера столбцов матрицы B с нулевыми элементами в последней строке, то последние элементы всех столбцов d_{j} матрицы (e, B)U не равны нулю. Системы (4) принимают вид

(I′, d_{j})x_{j} = b, j = 0, 1, …, m − r, (6)

Algorithm RationalBasis[A, c]

Input: A ∈ R^{n×m}, c ∈ R^{n}

Output: Решения системы Ax = c или NIL if c = 0 then

{B, S, T, r} := DiagonalF ormH(A); Comment : B = (b_{1}, b_{2}, …, b_{m−r}) Comment : I_{m−r} = (e_{1}, e_{2}, …, e_{m−r}).

return (T ( \begin{matrix}−b_{1} \\ e_{1} \end{matrix} ), …, T ( \begin{matrix}−b_{m−r} \\ e_{m−r} \end{matrix} ))

else {B, b, S, T, P, r} := DiagonalF orm(A, c); Comment : B = (b_{1}, b_{2}, …, b_{m−r}),

Comment : b = ( \begin{matrix}b′ \\ β \end{matrix} ), b_{j} = ( \begin{matrix}b′ & j \\ β_{j} \end{matrix} )

Comment : Q = diag(P, I_{m−r}),

Comment : I_{m−r+1} = (e′_{0}, e′_{1}, …, e′_{m−r}). if r < rank(A, c) then return NIL; else J := ∅

for j := 1 to m − r do

if β_{j} = 0 then J := J ∪ {j}; for j := 1 to m − r do

if j ∉ J

then ξ_{j} := β/β_{j}, f_{j} := e′_{j} else ξ_{j} := β, f_{j} := e′_{j} + e′_{0}

return (TQ( \begin{matrix}b′ \\ βe′_{0} \end{matrix} ),

TQ( \begin{matrix}b′ − ξ_{j}b′_{j} \\ ξ_{j}f_{j} \end{matrix} ), j = 1, …, m − r)

Рис. 1. Алгоритм "Рациональный базис"

при этом det(I′, d_{j}) ≠ 0. Так как у вектора b последний элемент не равен нулю, то решения

систем x_{j} = ( \begin{matrix}x′ & j \\ ξ_{j} \end{matrix} ) таковы, что ξ_{j} ≠ 0. По-

этому векторы ξ_{j}e′_{j}, j = 1, …, m − r линейно независимы, следовательно, и векторы (5) линейно независимы.

Следствие 2 Пусть B = (b_{1}, b_{2}, …, b_{m−r}),

b = ( \begin{matrix}b′ \\ β \end{matrix} ), b_{j} = ( \begin{matrix}b′ & j \\ β_{j} \end{matrix} ). ξ_{j} = β/β_{j}, f_{j} = e′_{j}

для j ∉ J. ξ_{j} = β, f_{j} = e′_{j} + e′_{0} для j ∈ J. Тогда

T Q( \begin{matrix}b′ \\ βe′_{0} \end{matrix} ), T Q( \begin{matrix}b′ − ξ_{j}b′_{j} \\ ξ_{j}f_{j} \end{matrix} ), (7)

j=1,…,m-r,

– базисное множество решений неоднородной системы линейных уравнений (1).

Доказательство. В (5) подставим решения систем (6). При этом будем учитывать, что, по

построению d_{j} = ( \begin{matrix}b′ & j \\ δ_{j} \end{matrix} ), δ_{j} = 1 при j ∈ J и

δ_{j} = β_{j} при j ∉ J. Затем умножим на матрицу W.

Cледствия теорем 1 и 2 позволяют записать алгоритм построения базисного множества решений системы линейных уравнений для произвольной коммутативной области R.

Алгоритм приведен на рис. 1. В этом алгоритме полагается следующее.

R – коммутативная область с единицей, в которой cуществует алгоритм сокращения (вычисляющий сомножитель b, по известным a и c, a = bc, c ≠ 0), (1) – система линейных уравнений над R.

DiagonalF orm(A, c) – алгоритм диагонализации матрицы (A, c) для неоднородной системы (1). Он вычисляет {B, b, T, S, P, r}, при этом r = rankA, а матрицы B, b, T, S, P описаны выше.

DiagonalF ormH(A) – алгоритм диагонализации матрицы A, который вычисляет {B, T, S, r} для однородной системы (1).

Для этих алгоритмов можно взять любой из методов [3], [4], [5], [6] со сложность O(n^{2}m) или метод [7] со сложностью O(n^{β−1}m).

РЕШЕНИЕ СИСТЕМЫ С ПОМОЩЬЮ ПОДЪЕМА ПО ГЕНЗЕЛЮ

Теоремы 1 и 2 сводят задачу построения базисного множества решений системы (1) к задаче решения нескольких определенных систем. Для вычисления базисного множества решений неоднородной системы необходимо решить m − r + 1 систем с матрицей коэффициентов размера r × r (для однородных систем нужно решить m − r систем). Для решения определенных систем можно воспользоваться p-адическим подъемом.

Напомним общую схему p-адических методов. Выбирается подходящий простой элемент p кольца R. Для решения системы необходимо, чтобы элемент p не делил определитель матрицы коэффициентов. Например, в кольце Z этот выбор может оказаться неудачным с вероятностью не более 1/p. Если проверка покажет, что решение неверное, то выбирается другой простой элемент p. Кольцо вычетов по простому модулю p является полем, и решение в этом поле можно искать, например, методом Гаусса. С помощью неравенства Адамара находится верхняя оценка для числителей и знаменателей решения системы, а по ним находится оценка для величины p^{k} – границы подъема. Затем решение поднимается по модулю p до p^{k}, и затем реконструируется рациональное решение. Для решения одной определенной системы используем алгоритм, приведенный

Algorithm RationalBasisH[A]

Input: A ∈ R^{n×m}

Output: Решения системы Ax = 0.

p := F irstP rime;

repeat Initialization: {Ã_{0}^{−1}, S, T, r} := DiagonalF ormH[mod_{p} A]; if m = r then return {∅};

( \begin{matrix}A_{0} & A_{1} \\ A_{2} & A_{3} \end{matrix} ) := SAT;

(A′_{0}, a_{0}) := A_{0};

(a_{1}, …, a_{m−r}) := A_{1}; j := 1; T est := T RUE; repeat \underline{Lifting:} l := 1; x := 0; c^{∗} := −a_{j}; repeat

x̃ := mod_{p}Ã_{0}^{−1}c^{∗}; c^{∗} := (c^{∗} − A_{0}x̃)/p; x := x + lx̃; l := lp; until ∥l∥ ≥ L_{j}^{H}[A];

( \begin{matrix}x′ & j \\ ξ_{j} \end{matrix} ) :=

RationalReconstruction[x, N_{j}^{H}, D]; \underline{Test:} Comment : (e_{1}, …, e_{m−r}) = I_{m−r}

x_{j} := T ( \begin{matrix}x′ & j \\ ξ_{j}e_{j} \end{matrix} );

if Ax_{j} ≠ 0 then T est := F ALSE; j := j + 1;

until j > m − r or T est = F ALSE; p := NextP rime[p]; until T est = T RUE; return {x_{1}, …, x_{m−r}}.

Рис. 2. Алгоритм "Вычисление рационального базиса однородной системы линейных уравнений с помощью p-адического подъема".

Диксоном [9] и использующий линейный pадический подъем.

Алгоритмы построения базисного множества решений системы (1), использующие линейный p-адический подъем, приведены на рис. 2 и 3. Алгоритм для решения однородных систем приведен на рис. 2, алгоритм для решения неоднородных систем — на рис. 3.

В алгоритмах используются следующие функции.

Функция ∥ ∥ имеет следующее значение: ∥α∥ = |α| для α ∈ Z, ∥α∥ = deg α для α ∈ F [x], для матрицы A ∥A∥ = max_{α∈A} ∥α∥, для столбца k матрицы A ∥a_{k}∥ = max_{i} ∥a_{i,k}∥.

Из неравенства Адамара следуют оценки для числителя и знаменателя решения системы.

Для кольца Z получим:

Algorithm RationalBasis[A, c]

Input: A ∈ R^{n×m}, c ∈ R^{n} Output: Решения системы Ax = c или NIL p := F irstP rime; i := 1; repeat

Initialization: {Ã_{0}^{−1}, B̃, b̃, S, T, P, r, r_{c}} :=

DiagonalF orm[mod_{p}(A, c)]; if r < r_{c} then return NIL; Comment : B̃ = (b̃_{1}, b̃_{2}, …, b̃_{m−r})

Comment : b̃_{j} = ( b̃′_{j} ; β̃_{j} ) Q = diag(P, I_{m−r})

Comment : I_{m−r+1} = (e′_{0}, e′_{1}, …, e′_{m−r}).

( \begin{matrix}A_{0} & A_{1} \\ A_{2} & A_{3} \end{matrix} ) := SAT; ( \begin{matrix}c_{0} \\ c_{1} \end{matrix} ) := Sc;

ĉ := Pc_{0}; (A′_{0}, a_{0}) := PA_{0}P; (a_{1}, …, a_{m−r}) := PA_{1}; J := ∅; for j := 1 to m − r do

if β_{j} = 0 then J := J ∪ {j}; j := 0; Test = TRUE; repeat

if j ∈ J then A_{j} := (A′_{0}, a_{j} + a_{0}) else A_{j} := (A′_{0}, a_{j}) \underline{Lifting:}

Ã_{j}^{−1} := ( \begin{matrix}I_{r−1} & −b̃′_{j}β̃_{j}^{−1} \\ 0 & β̃_{j}^{−1} \end{matrix} )Ã_{0}^{−1}

l := 1; x := 0; c^{∗} := ĉ; repeat

x̃ := mod_{p}Ã_{j}^{−1}c^{∗}; c^{∗} := (c^{∗} − A_{j}x̃)/p; x := x + lx̃; l := lp;

until ∥l∥ ≥ L_{j}[A, c]

( \begin{matrix}x′ & j \\ ξ_{j} \end{matrix} ) :=

RationalReconstruction[x, N_{j}, D_{j}]; \underline{Test:}

if j = 0 then x_{0} := T( \begin{matrix}P ( \begin{matrix}x′ & 0 \\ ξ_{0} \end{matrix} ) \\ 0 \end{matrix} )

else if j ∈ J

then x_{j} := T( \begin{matrix}P ( \begin{matrix}x′ & j \\ ξ_{j} \end{matrix} ) \\ ξ_{j}e_{j} \end{matrix} )

else x_{j} := T( \begin{matrix}P ( \begin{matrix}x′ & j \\ 0 \end{matrix} ) \\ ξ_{j}e_{j} \end{matrix} );

if Ax_{j} ≠ c then Test := FALSE; j := j + 1;

until j > m − r or Test = FALSE i := i + 1; p := NextPrime[p];

until Test = TRUE or i = NumberOfPrimes if Test = FALSE then return NIL else return {x_{0}, x_{1}, …, x_{m−r}}.

Рис. 3. Алгоритм "Вычисление рационального базиса неоднородной системы линейных уравнений с помощью p-адического подъема". γ = ∏_{k=1}^{n−1} ∥a_{k}∥, λ = min_{k=1,...,n−1} ∥a_{k}∥, μ_{j} = min(λ, ∥a_{j}∥),

D_{j} = n^{n/2}γ∥a_{j}∥ — оценка для знаменателя решения, N_{j} = n^{n/2}∥a_{j}∥γ∥c∥μ_{j}^{−1} — оценка для числителя, L_{j}[A, c] = 2N_{j}D_{j} – оценка для границы p-адического подъема.

Для F [x] получим: D_{j} = ∥a_{j}∥ + ∑_{k=1}^{n−1} ∥a_{k}∥ — оценка для знаменателя, N_{j} = D_{j} + ∥c∥ + min_{k=1,...,n−1,j} ∥a_{k}∥ — оценка для числителя, L_{j}[A, c] = N_{j} + D_{j} — оценка для границы pадического подъема.

NextP rime[p] — генератор простых элементов кольца R, который возвращает следующий простой элемент, после p. F irstP rime и MinP rime обозначают первый простой элемент и наименьший по норме простой элемент.

NumberOfP rimes — количество простых

элементов ⌊( max_{j=n,n+1,...,m} D_{j})∥MinP rime∥^{−1}⌋.

Пусть R̃ = R/pR фактор-кольцо R по модулю pR.

DiagonalF orm[mod_{p}(A, c)] алгоритм диагонализации в кольце R̃. Результатом вычислений будет {Ã_{0}^{−1}, B̃, b̃, S, T, P, r, r_{c}}. Здесь все обозначения имеют прежний смысл. Знак тильды обозначает матрицы над факторкольцом R̃. r и r_{c} — это ранги матриц mod_{p}A и mod_{p}(A, c). Так как R̃ является полем, то для вычислений можно воспользоваться методом Гаусса или любым другим методом решения систем линейных уравнений над полем.

DiagonalF ormH[mod_{p} A] алгоритм диагонализации в R̃, вычисляет {Ã_{0}^{−1}, S, T, r}.

RationalReconstraction[x̃, N_{i}, D_{i}] – восстановление вектора решения в поле частных кольца R по решению x̃ в поле R/p^{k}R когда числитель и знаменатель ограничены числами N_{i} и D_{i}. Алгоритм описан в работе [11].

РЕШЕНИЕ СИСТЕМЫ В КОММУТАТИВНОЙ ОБЛАСТИ

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

Пусть α — непустое конечное подмножество в области R. Пересечение главных идеалов, порожденных элементами множества α, является главным идеалом ∪_{p∈α}(p). Порождающий элемент этого идеала обозначим (α), т. е. ((α)) = ∪_{p∈α}(p).

Пусть K – поле частных области R, x – вектор из пространства K^{m}, χ̃ ⊂ R – множество знаменателей компонент вектора x.

Определение. Знаменателем вектора x будем называть χ = НОК (χ̃).

Вектор x будем записывать в виде произведения x = xχ^{−1}, x ∈ R^{m}, χ ∈ R. Знаменатель вектора x будем обозначать DEN(x).

Пусть M = {x|x ∈ K^{m}, Ax = b} – множество всех решений системы (1) в K^{m}. Обозначим M_{D} = M ∪ R^{m} – множество решений системы, лежащих в модуле R^{m}. Будем называть M_{D} множеством диофантовых решений.

Для x = (x_{1}, x_{2}, …, x_{s}) и y = (y_{1}, y_{2}, …, y_{s}) будем обозначать ⟨x, y⟩ = ∑_{i=1}^{s} x_{i}y_{i}.

Теорема 3 Пусть {x_{i} = x_{i}χ_{i}^{−1}|i = 1, 2, …, h} – базисное множество решений неоднородной системы линейных уравнений (1) и пусть x̂ = (x_{1}, x_{2}, …, x_{h}), χ = (χ_{1}, χ_{2}, …, χ_{h}). Тогда

M = \{\frac{⟨x̂, q⟩}{⟨χ, q⟩}|q ∈ R^{h}, ⟨χ, q⟩ ≠ 0\}

Доказательство: Пусть u = (u_{1}, u_{2}, …, u_{h}) ∈ K^{h} такой, что ∑_{i=1}^{h} u_{i} = 1. Обозначим s_{i} = u_{i}χ_{i}^{−1}, g – знаменателей всех чисел {s_{i}|i = 1, 2, …, h}, q_{i} = gs_{i}, q_{i} ∈ R, i = 1, 2, …, h, q = (q_{1}, q_{2}, …, q_{h}).

Пусть x̂ = (x_{1}, x_{2}, …, x_{h}) и z = ⟨x̂, u⟩, следовательно, z ∈ M. Тогда получим

z = ⟨x̂, u⟩ = ∑_{i=1}^{h} x_{i}χ_{i}^{−1}u_{i} = ∑_{i=1}^{h} x_{i}s_{i} = ⟨x̂, q⟩g^{−1},

(7)

1 = ∑_{i=1}^{h} u_{i} = ∑_{i=1}^{h} χ_{i}χ_{i}^{−1}u_{i} = ∑_{i=1}^{h} χ_{i}s_{i} = ⟨χ, q⟩g^{−1}.

(8)

Поэтому

\frac{⟨x̂, q⟩}{⟨χ, q⟩} = z ∈ M.

Наоборот, если q ∈ R^{h}, ⟨χ, q⟩ = g ≠ 0 и z = \frac{⟨x̂,q⟩}{⟨χ,q⟩}, то, обозначая u_{i} = χ_{i}q_{i}g^{−1}, получим z =

⟨x, u⟩ и ∑_{i=1}^{h} u_{i} = 1. Поэтому z ∈ M.

Следствие 3 Множество всех знаменателей {DEN(x)|Ax = c, x ∈ K^{m}} решений системы Ax = c вместе с нулевым элементом образует идеал I_{A} в кольце R.

Следствие 4 Для того чтобы система (1) имела диофантовы решения, необходимо и достаточно, чтобы I_{A} = R.

Следствие 5 Пусть идеал, порожденный знаменателями базисного множества решений

x_{i} = x_{i}χ_{i}^{−1}, i = 1, 2, …, h системы (1), является единичным, χ = (χ_{1}, χ_{2}, …, χ_{h}). Тогда найдется такой вектор q = (q_{1}, q_{2}, …, q_{h}) ∈ R^{h}, что ⟨χ, q⟩ = 1 и z = ⟨x, q⟩ – диофантово решение системы (1). Если при этом q_{s} ≠ 0, то x_{1}, …, x_{s−1}, z, x_{s+1}, …, x_{h} – базисное множество решений этой системы.

Следствие 6 Cистема (1) имеет диофантовы решения, если идеал, порожденный знаменателями базисного множества решений, является единичным.

Следствие 7 Пусть идеал, порожденный знаменателями базисного множества решений x_{i} = x_{i}χ_{i}^{−1}, i = 1, 2, …, h системы (1), является единичным, χ = (χ_{1}, χ_{2}, …, χ_{h}). Тогда найдется такой вектор q = (q_{1}, q_{2}, …, q_{h}) ∈ R^{h}, что ⟨χ, q⟩ = 1 и z = ⟨x, q⟩ – диофантово решение системы (1). Если при этом q_{s} ≠ 0, то x_{1}, …, x_{s−1}, z, x_{s+1}, …, x_{h} – базисное множество решений этой системы.

Следствие 8 Пусть x_{i} = x_{i}χ_{i}^{−1}, i = 1, 2, …, h базисное множество решений системы (1) и χ_{1} = 1. Тогда x_{1}, x_{i} − x_{1}(χ_{i} − 1), i = 2, 3, …, h – линейно независимые диофантовы решения этой системы.

Доказательство: Так как все z_{i} = x_{i} − x_{1}(χ_{i} − 1), i = 2, 3, …, h принадлежат R^{m} и являются линейно независимыми вместе с x_{1}, остается показать, что z_{i} – решения системы (1). Так как χ_{1} = 1, то

z_{i} = x_{i} − x_{1}(χ_{i} − 1) = \frac{x_{1}(−χ_{i} + 1) + x_{i} · 1}{χ_{1}(−χ_{i} + 1) + χ_{i} · 1}.

Следовательно, по теореме, z_{i}, i = 2, 3, …, h, – решения системы Ax = b.

Определение. Будем называть диофантовым базисом решений системы Ax = – базисное множество решений этой системы, которое целиком лежит в области R.

Другими словами – это m + 1 − r линейно независимых решений неоднородной (m − r – для однородной) системы, принадлежащих R^{m}.

Следствия последней теоремы позволяют записать алгоритм вычисления диофантова базиса решений системы. Алгоритм приведен на рис. 4.

В этом алгоритме используются следующие функции.

Функция RationalBasis[A, c] вычисляет базисное множество {x_{1}χ_{1}^{−1}, x_{2}χ_{2}^{−1}, …, x_{h}χ_{h}^{−1}}

Algorithm DiophantineBasis[A, c]

Input: A ∈ R^{n×m}, c ∈ R^{n} Output: Диофантовы решения или NIL {h, x_{1}χ_{1}^{−1}, x_{2}χ_{2}^{−1}, …, x_{h}χ_{h}^{−1}} :=

RationalBasis[A, c]; if h = 0 then return NIL; (gcd, u_{1}, u_{2}, …, u_{h}) :=

ExdendedGCD[(χ_{1}, χ_{2}, …, χ_{h})]; if gcd ≠ 1 then return NIL else

k := 1; while u_{k} = 0 do k := k + 1; z_{k} := ∑_{i=1}^{h} x_{i}u_{i}; for i = 1, 2, …, k − 1, k + 1 …, h do

z_{i} := x_{i} − x_{k}(χ_{i} − 1);

return {z_{1}, z_{2}, …, z_{h}}.

Рис. 4. Алгоритм "Диофантов базис".

решений системы Ax = c и количество h векторов в нем.

Функция ExdendedGCD[(χ_{1}, χ_{2}, …, χ_{h})] по заданным h числам χ_{1}, χ_{2}, …, χ_{h} вычисляет порождающий элемент gcd идеала (χ_{1}, χ_{2}, …, χ_{h}), если этот идеал является главным, а также множители u_{1}, u_{2}, …, u_{h} ∈ R, такие, что gcd = ∑_{i=1}^{h} u_{i}χ_{i}. Если этот идеал не является главным, то gcd принимает значение NIL.

Если результатом является NIL, то необходимо выполнять следующую итерацию. Для следующей итерации нужно выбрать случайным образом новые матрицы перестановок S и T и вычислить новый рациональный базис. Алгоритм ExdendedGCD нужно применять к знаменателям всех рациональных решений, полученным ранее. Как доказано в [10], ожидаемое число случайных рациональных решений, которое необходимо для получения одного диофантова решения, будет O(log n + log log ∥(A, c)∥). Если rankA = n, то за одну итерацию, будет получено m − n + 1 рациональных решений. Поэтому ожидаемое количество итераций будет O((m − n + 1)^{−1}(log n + log log ∥(A, c)∥)).

ОЦЕНКА СЛОЖНОСТИ АЛГОРИТМОВ

Сложность вычисления диофантового базиса определяется количеством операций в алгоритмах 1 и 4. Сложность первого из них составляет O(mn^{β−1}), а второго O(m(m−r)+C_{G}) операций в кольце R. C_{G} — количество операций, необходимых для разложения единицы по m − r образующим единичного идеала в кольце R. Для такого разложения единицы можно воспользоваться, например, алгоритмом вычисления базиса Гребнера. Оценка сложности такого алгоритма в общем случае выходит за рамки настоящей статьи. Отметим только, что для евклидовых колец C_{G} — это количество операций в алгоритме Евклида, который вычисляет наибольший общий делитель для m−r чисел.

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

Известны оценки сложности алгоритма Диксона [10]. Для целых чисел Z сложность алгоритма Диксона для решения системы (1) при n = m и det A ≠ 0 ограничена числом D_{Z} = O(n^{3}(log n + log ∥A∥ + log p)^{2} + n(log ∥∥)^{2}) бит-операций. Для кольца многочленов F [x] над полем F сложность алгоритма Диксона ограничена числом D_{F [x]} = O(n^{3}(deg A + deg p)^{2} + n(deg b)^{2}) операций в поле F.

Для вычисления базиса, состоящего из m−r решений потребуется не более, чем C_{Z} = (m − r)D_{Z} и C_{F [x]} = (m − r)D_{F [x]} операций, соответственно для каждого из этих случаев.

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

Количесто операций при вычислении диофантова базиса определяется главным образом количеством операций необходимых для вычисления рационального базиса и количеством итераций. Ожидаемое количество итераций равно O((m − n + 1)^{−1}(log n + log log ∥(A, c)∥)), а количество операций будет соответственно C_{Z} и C_{F [x]} для Z и F [x] соответственно. Следовательно, при m − n << n количество операций для вычисления диофантова базиса от размеров системы зависит существенно по кубическому закону.

По сравнению с алгоритмом [10], среднее число вычислений обратной матрицы теперь в m − r раз меньше.

ЗАКЛЮЧЕНИЕ

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

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

Если использовать арифметику с s-битными словами, то метод Гаусса для решения системы с матрицей размера n × m ранга n потребует O(n^{2}ms^{2}) бит-операций.

Приведенный алгоритм построения рационального базиса решений требует для R = Z примерно O((m − n + 1)n^{2}mw^{2}) бит-операций. Где w = log n + log ∥A∥ + log p. Их отношение равно O((m − n + 1)(w/s)^{2}. В реальных задачах log n+log ∥A∥ меньше s, а p обычно выбирается так, что log p меньше s. Поэтому w/s – величина меньше 1. А искомое отношение – это число (m − n + 1).

Отсюда видно, что при m − n << n алгоритм вычисления рационального решения имеет такую же сложность, что и метод Гаусса для приближенной арифметики при использовании чисел с плавающей точкой. Однако с ростом разности m − n отношение растет линейно, как m − n.

ЛИТЕРАТУРА

1. Strassen. V. Gaussian Elimination is not optimal // Numerische Mathematik. 1969. V. 13. P. 354–356.

2. Coppersmith D. and Winograd S. Matrix multiplication via arithmetic progressions // Journal of Symbolic Computation. 1990. V.

9. P. 251–280.

3. Dodgson C.L. Condensation of determinants, being a new and brief method for computing their arithmetic values // Proc. Royal Soc. Lond. 1866. V. A.15. P. 150–155.

4. Bareiss E.H. Sylvester’s Identity and Multistep Integer-Preserving Gaussian Elimination // Math. Comp. 1968. V. 22 (103). P. 565–578.

5. Malaschonok G.I. Solution of a system of linear equations in an intergal domain // USSR J. of computational Mathematics and Mathematical Physics. 1983. V. 23. P. 1497-1500.

6. Malaschonok G.I. Algorithms for the solution of systems of linear equations in commutative rings. // Effective Methods in Algebraic Geometry / Ed. by T. Mora and C. Traverso, Progress in Mathematics 94, Birkhauser, 1991. P. 289–298.

7. Malaschonok G.I. Recursive Method for the Solution of systems of Linear Equations // Computational Mathematics ( Sydow A. Ed, Proceedings of the 15th IMACS World Congress, Vol. I, Berlin, August 1997 ). Wissenschaft & Technik Verlag. Berlin, 1997. P. 475–480.

8. Bareiss E.H. Computational solutions of matrix problems over an integral domain // J. Inst. Maths Applics. 1972. V. 10. P. 68– 104.

9. Dixon J. Exact solution of linear equations using p-adic expansions // Numer. Math. 1982. V. 40. P. 137–141.

10. Mulders T. and Storjohann A. Diophantine Linear System Solving // Proceedings of ISSAC’99: ACM International Symposium on Symbolic and Algebraic Computation. July 1999, Vancouver, Canada.

11. Collins G.E. and Encarnacion M.J. Efficient rational number reconstruction. Journal of Symbolic Computation. 1995. V.

20. P. 287–297.

Поступила в редакцию 1 февраля 2000 г.