Algorithms for the solution of systems of linear equations in commutative ring
aa r X i v : . [ c s . S C ] N ov Algorithms for the solution of systems of linearequations in commutative ring ∗ Gennadi.I.Malashonok
Abstract
Solution methods for linear equation systems in a commutative ringare discussed. Four methods are compared, in the setting of several dif-ferent rings: Dodgson’s method [1], Bareiss’s method [2] and two meth-ods of the author - method by forward and back-up procedures [3] anda one-pass method [4].We show that for the number of coefficient operations, or for thenumber of operations in the finite rings, or for modular computation inthe polynomial rings the one-pass method [4] is the best. The method offorward and back-up procedures [3] is the best for the polynomial ringswhen we make use of classical algorithms for polynomial operations.
Among the set of known algorithms for the solution of systems of linearequations, there is a subset, which allows us to carry out computationswithin the commutative ring generated by the coefficients of the system.Recently, interest in these algorithms grew due to computer algebra com-putations. These algorithms may be used (a) to find exact solutions ofsystems with numerical coefficients, (b) to solve systems over the ringsof polynomials with one or many variables over the integers or over thereals, (c) to solve systems over finite fields.Let Ax = a (1) ∗ This paper was published in the book
Effective Methods in Algebraic Geometry , ed. byT. Mora and C. Traverso, Progress in Mathematics 94, Birkhauser, 1991, 289–298. No partof this materials may be reproduced, stored in retrieval system, or transmitted, in any formwithout prior permission of the copyright owner. e the given system of linear equations over a commutative ring R ,A ∈ R n × ( m − , a ∈ R n , x ∈ R m − , n < m , with extended coefficientmatrix A ∗ = ( a ij ), i = 1 , . . . , n, j = 1 , . . . , m .The solution of such a system may be written according to Crammer’srule x i = ( δ n im − P m − j = n +1 x j δ nij ) /δ n , i = 1 . . . n, where x j , j = n + 1 . . . m ,are free variables and δ n = 0 . δ k =— a ij | , i = 1 . . . k, j = 1 . . . k, k =1 . . . n, denote corner minors of the matrix A of order k, δ kij denotesminors obtained after a substitution in the minors δ k of the column i bythe column j of the matrix A ∗ , k = 1 . . . n, i = 1 . . . k, j = 1 . . . m. We examine four algorithms [1]-[4] for solving system (1) over thefraction field of the ring R , assuming that R does not have zero divisorsand all corner minors δ k , k = 1 . . . n , of the matrix A are different fromzero. Each of the algorithms is in fact a method for omputing the minors δ n and δ nij in the ring R . For each of the algorithms we evaluate:1.The general time for the solution taking into consideration onlyarithmetic operations and assuming moreover, the execution time formultiplication, division and addition/subtraction of two operands, thefirst of which is a minor of order i and the second one, a minor of order j , will be M ij , D ij , A ij correspondingly.2. The exact number of operations of multiplications, division andaddition/subtraction over the coefficient of the system.3. The number of operations of multiplication/division ( M R ), when R = R [ x . . . x r ] is a ring of polynomials with r variables with realcoefficients and only one computer word is required for storing any oneof the coefficients.4. The number of operations of multiplication/division ( M Z ), when R = Z [ x . . . x r ] is a ring of polynomials with r variables with integercoefficients and these coefficient are stored in as many computer wordsas are needed.5. The number of operations of multiplication/division ( M ZM ), when R = Z [ x . . . x r ] is a ring of polynomials with r variables with integercoefficients, but for the solution the modular method is applied, whichis based on the remainder theorem. Dodgson’s algorithm
Dodgson’s algorithm [1], the first algorithm to be examined, appearedmore than 100 years before the others.The first part of the algorithm consists of n − a ij = a i − ,j − a ij − a i − ,j a i,j − , i = 2 . . . n, j = 2 . . . m, formed from four neighboring elements, located at the intersection oflines i − i and of columns j − j .In the k -th step, k = 2 . . . n − , according to formulaˆ a k +1 ij = (ˆ a ki − ,j − ˆ a kij − ˆ a ki − ,j ˆ a ki,j − ) / ˆ a k − i − ,j − , (2) i = k + 1 . . . n, j = k + 1 . . . m, all minors ˆ a k +1 ij of order k + 1 are computed, these minors are formedby the elements located at the intersection of the lines i − k . . . i andcolumns j − k . . . j of the matrix A ∗ .In this way the minors ˆ a nnj = δ nnj , j = n . . . m , will be computed.Note that, it is essential that all minors ˆ a kij , k = 1 . . . n − , i = k . . . n − j = k . . . m −
1, appearing in the denominator of expression (2),be different from 0. This, of course, narrows down the set of solvableproblems.The second part of Dodgson’s algorithm is feasible only under theassumption that m = n + 1 and all the searched for unknown variables x i belong to R . Then the value x n = δ nnm /δ nnn may be computed andsubstituted in the initial system. After eliminating the last equationand recalculating the column of free members, a system of order n − x n − . During this process, it suffices to recalculate onlythose minors in which the column of free members appears. This processcontinues until all solutions are obtained.Let us evaluate the computing time of the algorithm. The first partof the algorithm is executed in time T = ( n − m − M , + A , )+ n − X i =2 ( n − i )( m − i )(2 M i,i + A i, i + D i,i − ) . (3)The second part of the algorithm, when m = n + 1 , is executed in time T = n − X i =1 i ( M , + A , ) + + n − X j =1 j (2 M , + A , )+ − X i =2 i X j =1 (2 M i,i +1 + A i +1 , i +1 + D i +1 ,i − ) . In this we suppose that the time needed for the recalculation of one freemember is M , + A , . The forward procedure of Bareiss’s algorithm [2] differs from Dodgson’salgorithm only in the selection of the leading (pivoting) minors.In the first step all minors of second order are computed a ij = a a ij − a j a i , i = 2 . . . n, j = 2 . . . m, (4)which surround the corner element a . At the k -th step, k = 2 . . . n − , and according to the formula a k +1 ij = ( a kkk a kij − a kik a kkj ) /a k − k − ,k − , (5) i = k + 1 . . . n, j = k + 1 . . . m, all the minors a k +1 ij of order k + 1 are computed, which are formed bysurrounding the corner minor δ k by row i and column j , that is, minorswhich are formed by the elements, located at the intersection of row1 . . . k, i and of columns 1 . . . k, j . Obviously, a nnj = δ nnj , j = n . . . m ,holds.Comparing the above procedure with Dodgson’s algorithm, it is seenthat here, it suffices that only the corner minors δ k = a kkk , k = 1 . . . n − , be different from zero and zero divisors . In order to do so, they can becontrolled by choice of the pivot row or column.The back-up procedure consists of n − k -th step, k = 1 . . . n − , all the minors δ k +1 ij = ( a k +1 k +1 ,k +1 δ ki,j − a k +1 k +1 ,j δ ki,k +1 ) /a kk,k , (6) i = k + 1 . . . n, j = k + 1 . . . m, are computed.The computing time of the forward procedure is the same as that ofthe first part of Dodgson’s algorithm (3), T = T .The computing time of the back-up procedure is T = n − X i =1 i ( m − i − M i,i +1 + A i +1 , i +1 + D i +1 ,i ) . Algorithm of Forward and Back-upProcedures
The forward procedure in this algorithm [3] is the same as the forwardprocedure in Bareiss’s algorithm, and is based on formulae (4) and (5).The back-up procedure is more economical, and consists of im medi-ate computation of the values δ nij according to the formulae δ nij = a nnn a iij − P nk = i +1 a iik δ nkj a iii (7) i = n − , . . . , j = n + 1 . . . m. The computing time of the forward procedure is the same as that of thefirst part of Dodgson’s algorithm (3), T = T .The computing time of the back-up procedure is T = ( m − n ) n − X i =1 (( n + 1 − i ) M n,i + ( n − i ) A i + n,i + n + D i + n,i ) . Let us note that when m = n + 1, T — is a quantity of order n , and T — a quantity of order n . Algorithms [2] and [3] consist of two parts (two passes). They make zeroelements under the main diagonal of the coefficient matrix during thefirst pass. And during the second pass, they make zero elements up tothe main diagonal.Algorithm [4] requires one pass, consisting of n − δ j = a a j − a a j , j = 2 . . . m,δ j = a j a − a j a , j = 3 . . . m. In the k -th step, k = 2 . . . n − , the minors of order k + 1 are computedaccording to the formulae (and see the Appendix) δ k +1 k +1 ,j = a k +1 ,k +1 δ kkk − k X p =1 a k +1 ,p δ kpj , j = k + 1 . . . m, (8) k +1 ij = ( δ k +1 k +1 ,k +1 δ ki,j − δ k +1 k +1 ,j δ ki,k +1 ) /δ kk,k , (9) i = 1 . . . k, j = k + 2 . . . m. In this way, at the k -th step the coefficients of the first k + 1 equationsof the system take part.The general computing time of the solution is T = (2 m − M , + A , ) + n − X k =2 ( m − k )(( k + 1) M k, + kA k +1 ,k +1 )++ n − X k =1 k ( m − k − M k,k +1 + A k +1 , k +1 + D k +1 ,k ) . We begin the comparison of the algorithms considering the general num-ber of multiplication
N M m , divisions N M d and additions/btractions N M a , which are necessary for the solution of the system of linear equa-tions (1) of order n × m . Moreover, we will not make any assumptionsregarding the computational complexity of these operations; that is wewill consider that during the execution of the whole computational pro-cess, all multiplications of the coefficients are the same, as are the sameall divisions and all additions/subtractions.The quantity of operations, necessary for Bareiss’s algorithm will be N M m = 2 n m − n − nm + n,N M d = (2 n m − n − nm + 2 m + 3 n − / ,N M a = (2 n m − n − nm + n ) / . The quantity of operations, necessary for the algorithm of forward andback-up procedure, is
N M m = (9 n m − n − nm − n − m + 8 n ) / ,N M d = (3 n m − n − nm − n + 13 n − / N M a = (6 n m − n − nm + 3 n + n ) / . The quantity of operations, necessary for the one-pass algorithm, is
N M m = (9 n m − n − nm − m + 6 n ) / , M d = (3 n m − n − nm − m + 2 n + 12) / N M a = (6 n m − n − nm + 3 n + n ) / . In the case when the number of equations and unknowns in the systemis the same and equal to n , we can compare all four algorithms.quantity of operations (2 n − n − n )2 ( n − n − n − n − n )2 n − n ( n − n + n )2 ( n − n )2 (4 n +3 n − n − n − n +10 n − n +3 n − n )6 ( n +2 n − n − n − n +6)6 (2 n +3 n − n )6 In this way, according to this evaluation, the fourth algorithm (one-pass) is to be preferred. Bareiss’s and Dodgson’s algorithms are approxi-mately equal regarding the quantity of operations, and each one of themrequires three times more divisions and two times more multiplications,as does the one-pass algorithm. The third algorithm lies somewhere inbetween.If we evaluate according to the general quantity of multiplication anddivision operations, considering only the third power, then we obtain theevaluation 3 n / n / n : 2 n / . [ x , x . . . x r ] Let R be the ring of polynomials of r variables over an integral do-main and let us suppose that every coefficient a ij of the system (1) is apolynomial of degree p in each variable a ij = p X u =0 p X v =0 . . . p X w =0 a iju,v...w x u x v . . . x wr . (10)Then it is possible to define, how much time is required for the executionof the arithmetic operations over polynomials which are minors of order i and j of the matrix A ∗ A ij = ( jp + 1) r a ij ,M ij = ( ip + 1) r ( jp + 1) r ( m ij + a i + j,i + j ) ,D ij = ( ip − jp + 1) r ( d ij + ( jp + 1) r ( m i − j,j + a ii )) . ere we assume, that the classical algorithms for polynomial multiplica-tion and division are used. And also, we consider that the time necessaryfor execution of the arithmetic operations of the coefficients of the poly-nomials, - when the first operand is coefficient of the polynomial, whichis a minor of order i , and the second - of order j , - is m ij , d ij , a ij , for theoperations of multiplication, division and addition/subtraction, respec-tively.Let us evaluate the computing time for each of the four algorithms,considering that the coefficients of the polynomials are real numbersand each one is stored in one computer word. We will assume that a ij = 0 , m ij = d ij = 1 , A ij = 0 , M ij = i r j r p r , D ij = ( i − j ) r j r p r , andwe will consider only the leading terms in m and n : M R = ρn r (cid:18) m r + 1 − n r + 2 (cid:19) ,M R = ρn r (cid:18) m r + 1 − n + 3 m r + 2 + 3 n r + 3 + m − nr + 1 − m − nr + 2 (cid:19) ,M R = ρn r (cid:18) m r + 2 − n r + 3 (cid:19) + ρ (cid:18) mr + 2 − nr + 3 (cid:19) , where ρ = n r +2 p r .For m = n + 1 it is possible to compare all four algorithms: N R =3 σ, N R = (2 r + 3) σ, N R = 2 σ, N R = (2 r + 1) σ + ρn/ ( r + 2)( r + 3),where σ = 3 n r +3 p r / (2 r + 1)(2 r + 2)(2 r + 3). For r = 0 we obtainthe same evaluation as in the previous section. For r = 0 we obtain3 : (2 r + 3) : 2 : (2 r + 1). [ x , x . . . x r ] , Classical Case As before we suppose that every coefficient of the system is a polynomialof the form (10). However, the coefficients of these polynomials are nowintegers and each one of these coefficients a ijuv...w is stored in l computerwords. Then, the coefficients of the polynomial, which is a minor oforder i , are integers of length il of computing words.Under the assumption that classical algorithms are used for the arith-metic operations on these long integers, we obtain: a ij = 2 jla , m ij = ijl ( m + 2 a ), d ij = ( il − jl + 1)( d + jl ( m + 2 a )), where a, m, d are the ex-ecution time of the single-precision operations of addition/subtraction,multiplication, and division. ssuming that a = 0 , m = d = 1 , we obtain the following evaluationof the execution times of polynomial operations: M ij = ijl ( ijp ) r , D ij =( i − j ) r +1 j r +1 l p r , A ij = 0 . In this way, the evaluation of the time for solution will be the sameas that for the ring R = R [ x , x , . . . , x r ] (section 6), if we replace ev-erywhere r by r + 1 and p r by lp r .Therefore, for m = n + 1 we obtain N Z = 3 ψ, N Z = (2 r + 5) ψ, N Z =2 ψ, N Z = (2 r + 3) ψ, where ψ = 3 n r +5 p r l / (2 r + 3)(2 r + 4)(2 r + 5) ,l ≥ , r ≥ . [ x , x . . . x r ] , Modular Case Let us evaluate the time for the solution of the same problem, for the ringof polynomials with r variables with integer coefficients R = Z [ x . . . x r ],when the modular method is applied, based on the remainder theorem.In this case we will not take into consideration the operations for trans-forming the problem in the modular form and back again.It suffices to define the number of moduli, since the exact quantityof operations on the system coefficients for the case of a finite field hasalready been obtained in section 4.We will consider that every prime modulus m i is stored in exactly onecomputer word, so that, in order to be able to recapture the polynomialcoefficients, which are minors of order n , n ( l + log( np ) / m i ) moduliare needed. It is easy to see due to Hadamar’s inequality.Further, we need up moduli for each unknown x j , which appears withmaximal degree np . There are r such unknowns, and therefore, in all, µ = prn ( l + log( np ) / m i ) moduli are needed.If we now make use of the table in section 5, denote the time formodular multiplication by m and the time for modular division by d ,then not considering addition/subtraction and considering only leadingterms in n , we obtain for m = n + 1: N ZM = (6 m + 3 d ) ν, N ZM =(6 m + 3 d ) ν, N ZM = (4 m + 2 d ) ν, N ZM = (3 m + d ) ν , where ν = µn / . Conclusion
We see, that modular methods are better then non-modular ones, asusual. And the one-pass method [4] is the best for modular computation.The method of forward and back-up procedures [3] are better fornon-modular computation in polynomial rings. And Dodgson’s method
1] stands nearly it for such polynomial computations.
Appendix: Foundation of the One-passAlgorithm