Effective Matrix Methods in Commutative Domains ^{⋆}

Gennadi I. Malaschonok

Tambov State University, 392622 Tambov, Russia

                      malaschonok@math-iu.tambov.su

Abstract. Effective matrix methods for solving standard linear algebra problems in a commutative domains are discussed. Two of them are new. There are a methods for computing adjoined matrices and solving system of linear equations in a commutative domains.

1 Introduction

Let R be a commutative domain with identity, K the field of quotients of R. We assume that R is equipped with an algorithm allowing exact division. This means that if two elements a and b of R are given (a being different from zero) such that b = ac with c ∈ R , then this algorithm can exhibit the exact quotient c. Let R^{n×m} denote the set of n × m matrices with entries in R .

This paper is devoted to the review of effective matrix methods in the domain R for a solution of standard linear algebra problems. They are (1) multiplicating two matrices, (2) solving linear systems in K, (3) solving linear systems in R, (4) computing the adjoint matrix, (5) computing the matrix determinant,

(6) computing the characteristic polynomial of a matrix.

We shall estimate algorithms according to the total number of multiplication and division operations in the ring R.

(1). Multiplication of two matrices. Let O(n^{β}) be the number of multiplication operations, that are necessary for multiplication of square matrices of the order n. For the standard multiplication of matrices β = 3, for Strassen [21] algorithm β = log 7, and for the best on today algorithm β < 2.376 [9].

(2). Solving linear systems in K. Let Ax = b be a systems of linear equations whose coefficients belong to the commutative domain R: A ∈ R^{n×m}, b ∈ R^{n}, x ∈ K^{m}. The main method here is the so-called Gauss method with exact divisions with complexity O(n^{2}m) operations in R. First this method was published in the paper of Dodgson [11], and later it was developed in the works [5], [22], [14], [17], [18]. We adduce the asymptotic complexity of these methods: 1.5n^{2}m for Dodgson method [11], 1.5n^{2}m for Bareiss method [5], n^{3} for forward

⋆ This paper was published in the book Formal Power Series and Algebraic Combinatorics, (Ed. by D.Krob, A.A.Mikhalev, A.V.Mikhalev), Springer, 2000, 506-517. No part of this materials may be reproduced, stored in retrieval system, or transmitted, in any form without prior permission of the copyright owner.

and back-up procedures [14], 2/3n^{2}m for one-pass method [17], 7/12n^{2}m for generalized method [18]. A fast method of solving systems of linear equations over commutative domain is published in the article [19]. The complexity of this method is O(n^{β−1}m), the same as the complexity of matrix multiplication.

(3). Solving linear systems in R. Let Ax = b be a system of linear equations A ∈ R^{n×m}, b ∈ R^{n}. The problem is to find all the solutions x of this system in the module R^{m}. The particular cases of this problem are discussed in [15], [16]. A randomized algorithm for finding all the solutions in R is discassed in the section 3 of this paper. It is supposed that there exists an algorithm that is able to ascertain whether the finitely generated ideal I = (a_{1}, a_{2}, …, a_{s}) is unit or not. If I is unit, then there are calculated the coefficients k_{i} ∈ R in the expansion of the unit 1 = ∑_{i=1}^{s} k_{i}a_{i}. It is possible to take in the capacity of such algorithm the algorithm of the Gröbner bases computation in R.

(4). Computing the adjoint matrix. The best known method of computing the adjoint matrix in an arbitary commutative ring has the complexity O(n^{3}\sqrt{n} logn loglogn) operations of addition, substraction and multiplication [13]. If in a commutative ring the exact division is possible, then the best method has the complexity O(n^{3}) operations of multiplication and exact division [16], [17]. In this work in section 2 we suggest the method with the complexity, equal to the complexity of matrix multiplication, i.e. O(n^{β}).

(5). Computing the matrix determinant. The intermediate result of each algorithm [11], [5], [14], [17], [18], [19] for solving systems in K is the computation of the matrix determinant. So the asymptotic complexity of determinant computation for these methods is 1.5n^{3}, 1.5n^{3}, n^{3}, 2/3n^{3}, 7/12n^{3} respectively. For the methods [19] the complexity of determinant computation is O(n^{β}). The best method of computing the determinant of a matrix without divisions was published by Kaltofen [13]. Complexity of this method is n^{3}\sqrt{n} log n log log n.

(6). Computing the characteristic polynomial of a matrix. In the case of an arbitrary commutative ring, the best algorithms are the Chistov one [7], and the Improved Berkowitz Algorithm [1] with size O(n^{β+1} log n) . In the paper [2] there are described two new efficient methods with O(n^{3}) ring operations (addition, subtraction, multiplication and exact division). The first one is the Quasi-triangular method (with asymptotic multiplicative complexity 5/3n^{3}) and the second one is the Tri-diagonal method (with asymptotic multiplicative complexity 3n^{3}). As in the case of Hessenberg’s method [12], they proceed by reducing the given matrix A to a particular upper quasi-triangular (Hessenberg) form, similar to A.

Commutative domain of principal ideals. This is the basic application field. In section 4 we discuss the problem of solving linear systems in the principal ideals domain R and in the field of fractions K Let Ax = c be a system of linear equations, A ∈ R^{n×m}, c ∈ R^{n}.

Solving linear system in the field of fractions. The best known on today method for solving determined system Ax = c in the field of fractions in the case when R = Z, m = n, det A ≠ 0, is the Dixon method [10], which uses the linear p-adic lifting. Its complexity is O(n^{3}(log n + log ∥A∥ + log p)^{2}) bit operations.

If m > n, R = Z and Gauss method with exact divisions is used then solving system Ax = c in usual arithmetic needs O(mn^{4}(log n + log ∥A∥)^{2}) bit operations. ∥A∥ denotes the absolute value of the greatest coefficient of the system. Using the Chinese remaindering method may reduce the complexity up to O(mn^{3}(log n + log ∥A∥)^{2}) bit operations [6].

My approach to this problem uses p-adic lifting like in [10]. The complexity of the algorithm in the case of the ring Z is O((m−n+1)n^{2}m(log m+log ∥A∥+ log p)^{2}) bit operations.

Solving linear system in the principal ideal domain. In [20] there is given the randomized algorithm for finding one solution of a system in the domain R, in the cases when R = Z and R = F [x] (F[x] — ring of polynomials over a field). This method is based on the Dixon algorithm.

I suggest a randomized algorithm for finding all the solutions of a system in a commutative domain of principal ideals. It is based on using p-adic lifting. Its complexity is essentially cubic in the dimension of system like [20], but the number of matrix inversions is now m − n times less.

2 Adjoint Matrix Computation

2.1 Introduction

Let A = ( \begin{matrix}A & C \\ B & D \end{matrix} ) be the invertible matrix and A — its invertible block. It is

possible to factorize its inverse matrix A^{−1}:

( \begin{matrix}I & −A^{−1}C \\ 0 & I \end{matrix} ) ( \begin{matrix}I & 0 \\ 0 & (D − BA^{−1}C)^{−1} \end{matrix} ) ( \begin{matrix}I & 0 \\ −B & I \end{matrix} ) ( \begin{matrix}A^{−1} & 0 \\ 0 & I \end{matrix} ) .  (2.1)

Let A be a matrix of the order n = 2^{p}. If a block inversion by formula (2.1) is possible for the blocks up to the second order, then the computation of inverse matrix needs 2^{p−1} second order block inversions and 6 · 2^{p−k−1} multiplications of blocks the order 2^{k} (k = 1, 2, …, p − 1).

Using the Strassen [21] algorithm of matrix multiplication for such matrix inversion we need (21n^{β} − 6n)/5 multiplicative operations where β = log_{2} 7.

The similar method is proposed in this section.

2.2 The Theorems of Adjoint Matrix Factorization

Let R be a commutative ring, A = (a_{i,j})—a square matrix of an order n over the ring R. We denote by A_{i,j}^{k} its square submatrix of the order k, obtained by the bordering of upper left block of order k − 1 by the row i and the column j, (i, j > k). Its determinant is denoted by a_{i,j}^{k} = detA_{i,j}^{k}. Denote the corner minor of the order k by δ_{k} = a_{k,k}^{k}. The determinant of the matrix, obtained from the upper left block A_{k,k}^{k} of order k by the replacement of the column i by the column j is denoted by δ_{k(i,j)}, (1 ≤ i ≤ k, k < j ≤ n).

Consider the matrices

A_{t}^{(s)} = (a_{i,j}^{s})_{j=s,...,t}^{i=s,...,t} and G_{s}^{(t)} = (δ_{t(i,j)})_{j=t+1,...,n}^{i=s,...,t}  (2.2)

of the order (t − s + 1) × (t − s + 1) and (t − s + 1) × (n − t), correspondingly.

With the preceding notation the determinant Sylvester identity [3] may be written in the following way:

detA_{t}^{(s)} = δ_{s−1}^{t−s}δ_{t}  (2.3)

where 1 < s < t ≤ n.

Let us prove the two main theorems of the adjoint matrix factorization.

Theorem 1. Let a square matrix A of the order n over the ring R be decomposed into the blocks

A = ( \begin{matrix}A & C \\ B & D \end{matrix} ) ,

A is the square block of order s, (1 < s < n), whose determinant δ_{s} is not zero or zero divisor in R. Then the adjoint matrix A^{∗} can be factorized

( \begin{matrix}δ_{s}^{−1}δ_{n}I & −δ_{s}^{−1}F C \\ 0 & I \end{matrix} ) ( \begin{matrix}I & 0 \\ 0 & G \end{matrix} ) ( \begin{matrix}I & 0 \\ −B & δ_{s}I \end{matrix} ) ( \begin{matrix}F & 0 \\ 0 & I \end{matrix} ) ,  (2.4)

where F = A^{∗}, G = δ_{s}^{−n+s+1}A_{n}^{(s+1)∗}, I is the unit matrix and the following identity takes place:

A_{n}^{(s+1)} = δ_{s}D − BF C.  (2.5)

Theorem 2. Let the square matrix A_{n}^{(s+1)} of order n − s, (s > 0, n − s > 2), over the ring R be decomposed into the blocks

A_{n}^{(s+1)} = ( \begin{matrix}A & C \\ B & D \end{matrix} ) ,

where A is the square block of the order t − s, (s < t < n), δ_{s} and δ_{t} are not zeros or zero divisors in R. Then the matrix δ_{s}^{−n+s+1}A_{n}^{(s+1)∗} can be factorized:

( \begin{matrix}δ_{t}^{−1}δ_{n}I & −δ_{t}^{−1}FC \\ 0 & I \end{matrix} ) ( \begin{matrix}I & 0 \\ 0 & δ_{s}^{−1}G \end{matrix} ) ( \begin{matrix}I & 0 \\ −B & δ_{t}I \end{matrix} ) ( \begin{matrix}F & 0 \\ 0 & I \end{matrix} ) ,  (2.7)

where F = δ_{s}^{−t+s+1}A_{t}^{(s+1)∗}, G = δ_{t}^{−n+t+1}A_{n}^{(t+1)∗}, I is the unit matrix and the following identity is true:

A_{n}^{(t+1)} = δ_{s}^{−1}(δ_{t}D − BFC)  (2.8)

.

Proof. Calculate the products of matrix A_{n}^{(s+1)} by the factors of the matrix (2.7) step by step from the right to the left:

A_{n}^{(s+1)} → ( \begin{matrix}δ_{t}I & FC \\ B & D \end{matrix} ) → ( \begin{matrix}δ_{t}I & FC \\ 0 & δ_{t}D − BFC \end{matrix} ) → ( \begin{matrix}δ_{t}I & FC \\ 0 & δ_{n}I \end{matrix} ) → δ_{n}I.

It is necessary to prove the identity (2.8) and the following identities:

FA = δ_{t}I  (2.9)
δ_{t}^{−n+t+1}A_{n}^{(t+1)∗}A_{n}^{(t+1)} = δ_{n}I  (2.10)

As A = A_{t}^{(s+1)}, the equality (2.9) follows from the determinant Sylvester identity detA_{t}^{(s+1)} = δ_{s}^{t−s−1}δ_{t}.

The identity (2.10) follows from the determinant Sylvester identity

detA_{n}^{(t+1)} = δ_{t}^{n−t−1}δ_{n}.

Let us prove the identity (2.8). Denote by a_{l,j}^{(s+1)∗} the cofactor of the element

(l, j) in the matrix A_{t}^{(s+1)}. From the determinant Sylvester identity we obtain:

δ_{s}^{t−s−1}δ_{t(i,j)} = ∑_{l=s+1}^{t} a_{l,i}^{(s+1)∗}a_{l,j}^{s+1}.

Since F = δ_{s}^{−t+s+1}A_{t}^{(s+1)∗} and C is the block of the matrix A_{n}^{(s+1)}, the last equality for the elements implies the matrix identity G_{t}^{(s)} = FC. We decompose the determinant of the matrix A_{i,j}^{s+1} according to the last column and obtain

a_{i,j}^{s+1} = a_{i,j}δ_{s} − ∑_{l=1}^{s} a_{l,j}σ_{s(l,j)},  (2.11)

where σ_{s(l,j)} is the determinant of the matrix, resulted from A_{s,s}^{s} by the replacement of the row l by the row j. Let σ = (σ_{s(1,i)}, σ_{s(2,i)}, …, σ_{s(s,i)}, 0, 0, …, 0), α = (a_{1,j}, a_{2,j}, …, a_{t,j}), β = (a_{i,1}^{s+1}, a_{i,2}^{s+1}, …, a_{i,t}^{s+1}) denote the rows with t elements. Then according to (2.11) we obtain the matrix identity

( \begin{matrix}I & 0 \\ −σ & δ_{s} \end{matrix} ) · A_{i,j}^{t+1} = ( \begin{matrix}A_{t,t}^{t} & α^{T} \\ β & a_{i,j}^{s+1} \end{matrix} ) .

Correspondingly, we write the following determinant identity, where the determinant of the matrix on the right is decomposed according to the last row:

δ_{s}a_{i,j}^{t+1} = δ_{t}a_{i,j}^{s+1} − ∑_{l=1}^{t} a_{i,l}^{s+1}δ_{s(l,j)}.

In the matrix form it may be written as δ_{s}A_{n}^{(t+1)} = δ_{t}D − BG_{t}^{(s)}. Taking into account G_{t}^{(s)} = FC we get the identity (2.9).

2.3 The Estimate of the Complexity

The dimension of the upper left block A in the process of the factorization of the matrix may be taken arbitrarily. Consider the case, when the dimension of the block A is a degree of 2. We call such decomposition of the adjoin matrix the binary factorization.

Let M (n) be the complexity of the multiplication of two matrices of the order n and its asymptotical estimate is αn^{β}.

Then the complexity of the adjoint matrix calculation for the matrix of the order n = 2^{p} by means of binary factorization is C (n) = ∑_{k=0}^{p−2} 6 · 2^{k}M (2^{n−k−1}). We neglect the complexity of multiplication of a matrix by a scalar, i.e. the terms of the order n^{2}.

Therefore, the asymptotical estimate of the complexity of the adjoint matrix calculation is 3αn^{β}/(1 − 2^{1−β}).

Finally, for the relation of the complexities of the adjoint matrix calculation and of the matrix multiplication we obtain the asymptotical estimate k(β) = lim_{n→∞} \frac{C (n)}{M (n)} = \frac{3}{1−2^{1−β}} . For example we have k(3) = 4 for classical multiplication, and k(log_{2} 7) = 4.2 for Strassen’s multiplication.

3 Linear System Solving in Commutative Domains

3.1 Solving Systems in a Field of Fractions

Let R be a commutative domain with an identity, K be a field of fractions of R, A ∈ R^{n×m}, rankA = r, c ∈ R^{n},

Ax = c  (3.1)

be a system of linear equations. Let S and T be permutation matrices, which transpose linearly independent rows and columns of the matrix A to the upper left corner. We obtain in this corner a square matrix of size r × r, denote it by A_{0} (det A_{0} ≠ 0). The matrices SAT and Sc may be written in a block form:

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}.

Denote by M = {x | x ∈ K^{m}, Ax = c} the set of all the solutions of the system (3.1).

We denote by I_{r} the identity matrix of order r, E_{i,j}—square matrices which have only one nonzero element—(i, j), that equals 1.

We need some facts from the theory of linear equations.

1. If rank(A, c) ≠ r, then M = ∅. If rank(A, c) = r, then M is a hyperplane of dimension m − r in a space of dimension m. It is defined by m − r + 1 points, which do not belong to one hyperplane of lower dimension.

2. If the system (3.1) is homogeneous (c = 0) and x_{1}, x_{2}, …, x_{m−r} are its linearly independent solutions then M = {∑_{i=1}^{m−r+1} x_{i}u_{i} | u_{i} ∈ K}.

3. If the system (3.1) is nonhomogeneous (c ≠ 0) and x_{1}, x_{2}, …, x_{m−r+1} are its linearly independent solutions then M = {∑_{i=1}^{m−r+1} x_{i}u_{i} | u_{i} ∈ K, ∑_{i=1}^{m−r+1} u_{i} = 1}.

Definition. A basis set of solutions of a homogeneous system of linear equations (3.1) is a set which consists from m − r linearly independent solutions of the system (3.1). A basis set of solutions of a nonhomogeneous system of linear equations (3.1) is a set which consists from m − r + 1 linearly independent solutions of the system (3.1).

The next two theorems reduce the problem of getting the basis solutions of (3.1) to several problems of solving determined systems. The first theorem considers homogeneous systems, the second — nonhomogeneous systems.

Theorem 3. Let (3.1) be a homogeneous system of linear equations and A_{1} = (a_{1}, a_{2}, …, a_{m−r}), a_{j} ∈ R^{r}. Then the systems

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

are determined and their solutions x_{j} ∈ K^{r} define the basis set of solutions of (3.1):

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

where e_{j} ∈ R^{m−r} are the columns of the identity matrix I_{m−r} = (e_{1}, e_{2}, …, e_{m−r}).

Proof. Denote y = T^{−1}x. By the condition we have (A_{0}, A_{1})y = 0. We search

for the solution in the form y = ( \begin{matrix}x_{j} \\ e_{j} \end{matrix} ), and get the system (3.2).

Corollary 1. Let it be B = A_{0}^{−1}A_{1} = (b_{1}, b_{2}, …, b_{m−r}), b_{j} ∈ K^{r}. Then

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} ) is the basis set of solutions of (3.1).

Theorem 4. Let (3.1) be a nonhomogeneous system of linear equations, P a permutation matrix such that the last element of the vector b = P A_{0}^{−1}c_{0} is not equal to 0. Let it be B = P A_{0}^{−1}A_{1}, J ⊂ {1, …, m − r} be the numbers of the columns of the matrix B with zero elements in the last row. Let it be U = I_{m−r+1} + ∑_{j∈J} E_{1,j+1}, W = diag(I_{r−1}, U), Q = diag(P, I_{m−r}), V = QW, I_{m−r+1} = (e′_{0}, e′_{1}, …, e′_{m−r}), where e′_{j} ∈ R^{m−r+1} are the columns of the unit matrix, and (A′_{0}, a_{0}, a_{1}, a_{2}, …, a_{m−r})=P (A_{0}, A_{1})V where A′_{0} ∈ R^{r×r−1}, a_{j} ∈ R^{r}. Then the systems

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

are determined. The solutions of these systems x_{j} = ( \begin{matrix}x′ & j \\ ξ_{j} \end{matrix} ), x′_{j} ∈ K^{r−1}, ξ_{j} ∈ K

define the basis set of solutions of (3.1):

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

Proof. Denote y = V^{−1}T^{−1}x. By the condition we have P (A_{0}, A_{1})V y = P c_{0} and P (A_{0}, A_{1})V = (A′_{0}, a_{0}, a_{1}, …, a_{m−r}). If we search for the solution in the

form y = ( \begin{matrix}x′ & j \\ ξ_{j}e′_{j} \end{matrix} ), then we obtain the systems (3.4).

Let us show that the systems (3.4) are determined and ξ_{j} ≠ 0. Multiply them by P A_{0}^{−1}P from the left. Since P = P^{−1}, we get P A_{0}^{−1}P (A′_{0}, a_{j})x_{j} = P A_{0}^{−1}c_{0}, j = 0, 1, …, m − r.

As b = P A_{0}^{−1}c_{0} and P A_{0}^{−1}P (A′_{0}, a_{0}) = I_{r}, the first of the systems (3.4) gets the form I_{r}x_{0} = b. Denote I_{r} = (I′, e), e ∈ R^{r}. We see that P A_{0}^{−1}P A′_{0} = I′, P A_{0}^{−1}P a_{0} = e.

Since P A_{0}^{−1}P (a_{0}, …, a_{m−r}) = P A_{0}^{−1}P (a_{0}, P A_{1})U = (e, B)U, we get d_{j} = P A_{0}^{−1}P a_{j}, j = 1, …, m − r are the columns of the matrix (e, B)U. As U = I_{m−r+1} + ∑_{j∈J} E_{1,j+1} and J are the numbers of the columns of the matrix B with zero elements in the last row, the last elements of all columns d_{j} of the matrix (e, B)U do not equal zero. The systems (3.4) obtain the form

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

and det(I′, d_{j}) ≠ 0. Since the last element of the vector b is nonzero, solutions

x_{j} = ( \begin{matrix}x′ & j \\ ξ_{j} \end{matrix} ) of the systems (3.6) are such that ξ_{j} ≠ 0. So the vectors ξ_{j}e′_{j},

j = 1, …, m−r, are linearly independent, therefore the vectors (3.5) are linearly independent.

Corollary 2. Let it be 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} for j ∉ J; ξ_{j} = β, f_{j} = e′_{j} + e′_{0} for j ∈ J. Then

T Q ( \begin{matrix}b′ \\ βe′_{0} \end{matrix} ) , T Q ( \begin{matrix}b′ − ξ_{j}b′_{j} \\ ξ_{j}f_{j} \end{matrix} ) , j = 1, …, m − r,  (3.7)

is the basis set of solutions of (3.1).

Proof. Substitute solutions of (3.6) into (3.5). We take into account that accord-

ing to the construction, d_{j} = ( \begin{matrix}b′ & j \\ δ_{j} \end{matrix} ), δ_{j} = 1 for j ∈ J and δ_{j} = β_{j} for j ∉ J.

Then we multiply by the matrix W.

Corollaries 1 and 2 allow to present the algorithm to construct the basis set of solutions of a system of linear equations for an arbitrary commutative domain.

3.2 System Solving in a Domain

Now we consider a linear system solving in a commutative domain. Note that it is not a problem for homogeneous systems, since any solution in a field of fractions, been multiplied by a suitable factor, gives the solution in the domain. So further we shall consider only nonhomogeneous systems.

Let α be a nonempty finite subset of R. The intersection of principal ideals generated by elements of the set α, is the principal ideal ∪_{p∈α}(p). We denote by LCM(α) the generator of this ideal, i.e. (LCM(α)) = ∪_{p∈α}(p).

Let K be a field of fractions of R, x be a vector of the space K^{m}, α_{x} ⊂ R be a set of denominators of components of x.

Definition. A denominator of a vector x is χ = LCM(α_{x}). The vector x will be written as a product x = xχ^{−1}, x ∈ R^{m}, χ ∈ R.

The denominator of a vector x is denoted by DEN(x).

Let M = {x | x ∈ K^{m}, Ax = b} be the set of all solutions of (3.1) in K^{m}. Denote by M_{D} = M ∪ R^{m} the set of solutions lying in the module R^{m}. We call M_{D} the set of Diophantine solutions.

For x=(x_{1}, x_{2}, …, x_{s}) and y =(y_{1}, y_{2}, …, y_{s}) we denote ⟨x, y⟩ = ∑_{i=1}^{s} x_{i}y_{i}.

Theorem 5. Let {x_{i} = x_{i}χ_{i}^{−1} | i = 1, 2, …, h} be a basis set of solutions of the nonhomogeneous system (3.1), x̂ = (x_{1}, x_{2}, …, x_{h}), χ = (χ_{1}, χ_{2}, …, χ_{h}). Then

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

Proof. Let u = (u_{1}, u_{2}, …, u_{h}) ∈ K^{h} be such that ∑_{i=1}^{s} u_{i} = 1. Further we use the notations: s_{i} = u_{i}χ_{i}^{−1}, g is LCM of denominators of all numbers {s_{i} | i = 1, 2, …, h}, q_{i} = gs_{i} ∈ R, i = 1, 2, …, h, q = (q_{1}, q_{2}, …, q_{h}).

Let it be x̂ = (x_{1}, x_{2}, …, x_{h}), z = ⟨x̂, u⟩, ∑_{i=1}^{s} u_{i} = 1, and so z ∈ M. Then we obtain

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

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

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

Conversely, if q ∈ R^{h}, ⟨χ, q⟩ = g ≠ 0 and z = \frac{⟨x̂,q⟩}{⟨χ,q⟩}, then, denoting u_{i} = χ_{i}q_{i}g^{−1}, obtain z = ⟨x̂, u⟩ and ∑_{i=1}^{s} u_{i} = 1. Therefore z ∈ M.

Corollary 3. The set I_{A}={DEN(x) | Ax = c, x ∈ K^{m}} ∪ 0 is an ideal in R.

Corollary 4. M_{D} ≠ ∅ if and only if I_{A} = R.

Corollary 5. A system (3.1) has Diophantine solutions if the ideal generated by the denominators of a basis set of solutions is unit.

Corollary 6. Let the ideal, generated by the denominators of basis solutions x_{i} = x_{i}χ_{i}^{−1}, i = 1, 2, …, h of the system (3.1) be unit, χ = (χ_{1}, χ_{2}, …, χ_{h}). Then there exists a nonzero vector q = (q_{1}, q_{2}, …, q_{h}) ∈ R^{h}, such that ⟨χ, q⟩ = 1 and ⟨x, q⟩ is Diophantine solution of (3.1). If in addition q_{s} ≠ 0, then

x_{1}, …, x_{s−1}, ⟨x, q⟩, x_{s+1}, …, x_{h}

is a basis set of solutions for this system.

Definition. We call a Diophantine basis of solutions for a system Ax = b a basis set of solutions for this system, that wholly belongs to R^{m}.

In other words a Diophantine basis consists of the m+1−r linearly independent solutions of nonhomogeneous (m − r for homogeneous) system, that belong to R^{m}.

Corollary 7. Let x_{i} = x_{i}χ_{i}^{−1}, i = 1, 2, …, h be a basis set of solutions of (1) and χ_{1} = 1. Then the set x_{1}, x_{i} − x_{1}(χ_{i} − 1), i = 2, 3, …, h, is a Diophantine basis of solutions for this system.

Proof. : Since all z_{i} = x_{i}−x_{1}(χ_{i}−1), i = 2, 3, …, h belong to R^{m} and are linearly independent together with x_{1}, it remains to show that z_{i} are the solutions of (1). As χ_{1} = 1, we have

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

Therefore, by Theorem 3, z_{i}, i = 2, 3, …, h are the solutions of the system Ax = b.

Corollaries of Theorem 5 allow to present a rendomized algorithm for computing a Diophantine basis of a system solutions.

If the ideal generated by the denominators of a basis set of solutions is unit, than we compute a Diophantine basis. Else we must choose a new permutation matrices S and T and compute a new rational basis, and so on. So we get some kind of iterative Diophantine solve. As it is proved in [20] an expected number of rational solutions that are necessary for getting Diophantine solution is O(log n + log log ∥(A, c)∥) for the ring Z. One of the question here is how to get such evaluation for another rings.

As we get m − n + 1 rational solutions in one iteration for rankA = n, so an expected number of iterations for the ring Z is O((m − n + 1)^{−1}(log n + log log ∥(A, c)∥)).

4 Solving Linear System in Principal Ideal Domain

4.1 Solving Linear System in the Field of Fractions

Theorems 3 and 4 reduce the problem of getting the basis solutions of the system (3.1) to several problems of solving determined systems. Computing a basis set of solutions of nonhomogeneous system needs to solve m − r + 1 systems with the matrix of coefficients of size r × r (the case of homogeneous systems needs to solve m − r systems). To solve a determined system one may use the p-adic lifting.

Recall a general scheme of p-adic methods. The suitable prime element p of the ring R is chosen. The element p must not divide the determinant of the coefficients matrix. For example, in the ring Z this choice may appear to be unsuccessful with probability no more than 1/p. If the check shows that the solution is incorrect, then another prime element p is chosen. The ring of residues prime modulo p is a field, and the solution in this field may be found for example using the Gauss method. The upper evaluations for numerators and denominators of system solutions is calculated by means of Hadamard inequality. According to them the upper evaluation for p^{k} — the boundary of lifting. Then the solution is lifted modulo p up to p^{k} and a rational solution is constructed. For solving one determined system we apply the algorithm given by Dixon [10], which use a linear p-adic lifting. Its complexity is O(n^{3}(log n + log ∥A∥ + log p)^{2}) bit operations. So the algorithm for getting all the rational solutions of the system Ax = c, when rankA = n, has the complexity O(n^{2}m(m − n + 1)(log n + log ∥A∥ + log p)^{2}) bit operations.

4.2 Solving Linear System in the Principal Ideal Domain

A one iterative step for computing a Diophantine basis is consist in the computing a rational basis and a Diophantine solutions, due to the Corollaries of Theorem 5.

The complexity of the first is O(mn^{β−1}), of the second — O(m(m−r)+C_{G}) operations in the ring R. C_{G} is the amount of operations that are necessary for an expansion of a unit in m − r generators of the unit ideal in the ring R. Such expansion of a unit may be obtained for example with the help of the algorithm of computing the Gröbner basis. The evaluation of the complexity of such algorithm in general is not the subject of this paper. Mention that for Euclidean rings C_{G} it is the amount of operations in Euclidean algorithm, that calculates the GCD for m − r numbers.

The more defined evaluations may be obtained for algorithms, using the linear p-adic lifting. It is known that the complexity of Dixon algorithm [20] in the case of integer numbers Z, for n = m and det A ≠ 0 is bounded by the number D_{Z} = O(n^{3}(log n + log ∥A∥ + log p)^{2} + n(log ∥∥)^{2}) bit operations. In the case of the ring of polynomials F [x] over a field F this complexity is bounded by D_{F [x]} = O(n^{3}(∥A∥ + ∥p∥)^{2} + n(∥c∥)^{2}) operations in the field F. The function ∥ ∥ has the next meaning: ∥α∥ = |α| for α ∈ Z, ∥α∥ = deg α for α ∈ F [x], for matrix A, ∥A∥ = max_{α∈A} ∥α∥.

Computing a basis from m−r solutions needs not more than C_{Z} = (m−r)D_{Z} and C_{F [x]} = (m − r)D_{F [x]} operations for each case correspondingly.

The average number of matrix inversions for computing one rational solution is now m − r times less than this number in the algorithm [20].

The complexity of computing a Diophantine basis for Z and F [x] is the same as the complexity of computing a rational basis. An expected number of iterations is O((m − n + 1)^{−1}(log n + log log ∥(A, c)∥)).

References

1. Abdeljaoued, J.: Berkowitz Algorithm, Maple and computing the characteristic polynomial in an arbitrary commutative ring. Computer Algebra MapleTech 4,

No. 3, Birkhauser Boston (1997)

2. Abdeljaoued, J., Malaschonok, G. I.: Efficient Algorithms for Computing the Characteristic Polynomial in a Domain. Journal of Pure and Applied Algebra (to appear)

3. Akritas, A. G., Akritas, E. K., Malaschonok, G. I.: Various proofs of Sylvester’s (determinant) identity. Mathematics and Computations in Simulation 42 (1996) 585–593

4. Akritas, A. G., Akritas, E. K., Malaschonok, G. I.: Matrix computation of subresultant polynomial remainder sequences in integral domain. Reliable Computing 1, No. 4 (1995) 375–381

5. Bareiss, E. H.: Sylvester’s Identity and Multistep Integer-Preserving Gaussian Elimination. Math. Comp. 22 (103) (1968) 565–578

6. Bareiss, E. H.: Computational solutions of matrix problems over an integral domain. J. Inst. Maths Applics 10 (1972) 68–104

7. Chistov, A. L.: Fast parallel calculation of the rank of matrices over a field of arbitrary characteristic. Proc. FCT’85, Springer Lecture Notes in Computer Science 199 (1985) 147–150

8. Collins, G. E. and Encarnacion, M. J.: Efficient rational number reconstruction. Journal of Symbolic Computation. 20 (1995) 287–297

9. Coppersmith, D. and Winograd, S.: Matrix multiplication via arithmetic progressions. Journal of Symbolic Computation 9 (1990) 251–280

10. Dixon, J.: Exact solution of linear equations using p-adic expansions. Numer. Math. 40 (1982) 137–141

11. Dodgson, C. L.: Condensation of determinants, being a new and brief method for computing their arithmetic values. Proc. Royal Soc. Lond. A.15 (1866) 150–155

12. Faddeev, D. K., Faddeeva, V. M.: Computational methods of linear algebra. W.H. Freeman & Co., San Francisco (1963)

13. Kaltofen, E.: On Computing Determinants of Matrices Without Divisions. In: Wang, P. S. (ed.): Proc. Internat. Symp. Symbolic Algebraic Comput. ISSAC’92. ACM Press (1992) 342–349

14. Malaschonok, G. I.: Solution of a system of linear equations in an integral domain. USSR Journal of computational Mathematics and Mathematical Physics 23 (1983) 1497–1500

15. Malaschonok, G. I.: System of linear equations over a commutative ring. Academy of Sciences of Ukraine, Lvov (1986) (in Russian)

16. Malaschonok, G. I.: On the solution of a linear equation system over commutative rung. Math. Notes of the Acad. Sci. USSR 42, No. 4 (1987) 543–548

17. Malaschonok, G. I.: Algorithms for the solution of systems of linear equations in commutative rings. In: Mora, T. and Traverso, C. (ed.): Effective Methods in Algebraic Geometry. Progress in Mathematics 94, Birkhauser (1991) 289–298

18. Malaschonok, G. I.: Argorithms for the Computing Determinants in a commutative rings. Diskretnaya Matematika 7, No. 4 (1995) 68–76

19. Malaschonok, G. I.: Recursive Method for the Solution of Systems of Linear Equations. In: Sydow, A. (ed.): Computational Mathematics (Proceedings of the 15th IMACS World Congress, Vol. I, Berlin, August 1997). Wissenschaft & Technik Verlag, Berlin (1997) 475–480

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

21. Strassen, V.: Gaussian Elimination is not optimal. Numerische Mathematik 13 (1969) 354–356 22. Sasaki, T. & Murao, H.: Efficient Gaussian elimination method for symbolic de-

terminants and linear systems. A.C.M. Trans. Math. Software 8, No. 4 (1968) 277–289