General solution of full row rank linear systems of equations using a new compression ABS model

This paper describes two new ABS algorithms based on two-step ABS methods for solving general solution of full row rank linear systems of equations. For both of our works, the ith iteration solves the first 2i equations, but for the second algorithm, we compress the space. We investigate the numerical stability of our models and a class of methods having optimal stability is defined. The condition for the residual perturbation to be minimal is also given. Computational complexity and numerical results demonstrate that, for our new methods, we need less work than corresponding two-step ABS algorithms and Huang’s method. The number of multiplications for these new schemes is half of the storage needed by the Gaussian elimination method. Our new version ABS algorithms are computationally better than the classical Gaussian elimination method, having the same arithmetic cost, but using less memory and no pivoting is necessary. The compression ABS model economizes the space more than the first algorithm, via deleting the zero rows of the Abaffian matrix by two in every iteration.


Introduction
ABS methods have been used broadly for solving linear and nonlinear systems of equations comprising large number of constraints and variables. In addition, ABS methods provide a unification of the field of finitely terminating methods for the solution of linear systems of equations. ABS methods were introduced by Abaffy, Broyden, and Spedicato initially for solving a determined or underdetermined linear system and later extended for linear least squares, nonlinear equations, optimization problems, Diophantine equations, and fuzzy linear systems [1,11,13]. These extended ABS algorithms offer some new approaches that are better than classical ones under several respects. In addition, extensive computational experience has shown tNow, since two new equations are hat ABS methods are implementable in a stable way, being often more accurate than the corresponding initial algorithm. ABS methods can be more effective than some of the other traditional methods. See more about ABS algorithms in [15][16][17][18]. A review of ABS algorithms is observed in [19]. The basic ABS algorithm works on a system of the form: where A ¼ ½a 1 ; . . .; a m T ; a i 2 R n ; 1 i m; x 2 R n ; b 2 R m ; m n: There is another notation for (1) as the following form: where A T ¼ ½a 1 ; . . .; a m T ; a i 2 R n ; 1 i m; x 2 R n ; b 2 R m ; m n: Notice that systems (1) and (2) are equivalent. Matrix computations are presented in [14]. The basic ABS methods determine solution of the above systems or signify lack of its existence in at most m iterates. Amini et al. proposed two-step ABS algorithms for the general solution of full row rank linear systems of equations [5,6]. The purpose of our paper is to present two new extended two-step ABS models and study on analysis of error propagation for them. The structure of this paper is organized as follows: Sect. 2 is devoted to the construction of a new two-step ABS model for solving general solution of full row rank linear systems of equations. Rank reducing process is done in two phases, for per iterate. The first phase helps us to have a solution for the ith iteration and the next phase leads to compute general solution of that iteration. In addition, we state and prove related theorems. We present our first twostep ABS algorithm for the general solution of full row rank linear systems of equations in this section. In Sect. 3, to compress required space, we present our second algorithm, such that the number of rows of the Abaffian matrix is reduced by two for per iterate. Thus, considering matrix multiplication, related parameters are selected from proper dimensions. In Sect. 4, we investigate the stability by the backward error analysis techniques due to Broyden's backward error analysis and Galantai's method. The errors in the steps are determined in the terms of projectors constructable from the conjugate directions. A class of methods having optimal stability characteristics is defined. Computational complexity and numerical results are discussed in Sect. 5, in detail. In Sect. 6, we summarize our achievements.

A new class of two-step ABS model
The basic ABS algorithm starts with an initial vector x 0 2 R n and a nonsingular matrix H 0 2 R nÂn (Spedicato's parameter). Given that x i is a solution of the first i equations, the ABS algorithm computes x i as the solution of the first i þ 1 equations as the following steps [1]: where w i 2 R n (Abaffy's parameter) is chosen, so that w T i H i a i 6 ¼ 0. Here, we are motivated to study on a method that satisfies two new equations at a time. We consider the system (1) under the assumption that A is full rank in row, i.e., rank(A) = m and m n. Suppose that m ¼ 2l (if m is odd, we can add a trivial equation to the system). Take A 2i ¼ ½a 1 ; . . .; a 2i T ; b 2i ¼ ½b 1 ; . . .; b 2i T and r j ðxÞ ¼ a T j x À b j ðj ¼ 1; . . .; mÞ: Assume that we are at the ith iteration and x i satisfies ( Suppose that r 2iÀ1 ðx iÀ1 Þ 6 ¼ 0 and r 2i ðx iÀ1 Þ 6 ¼ 0. Then, k i must be nonzero and the above systems are compatible if and only if we take where r 2iÀ1 ðx iÀ1 Þ ¼ r 2i ðx iÀ1 Þ ¼ r 2iÀ1 ðx iÀ1 Þr 2i ðx iÀ1 Þ. There are various ways to satisfy (3). We consider the following model: and Now, assume As relations (4)-(6), we will construct (7) and (8), such that and We compute H iþ1 from H l i , such that the relations (7) and (8) hold and proceed inductively. We define & Therefore, we have ( Thus, we must define g 2iþ1 ; d 2iþ1 2 R n in such a way that and H l i c j ¼ 0; j ¼ 1; . . .; 2i. The condition (9) is satisfied for j 2i À 1 by the induction hypothesis. By taking j ¼ 2i þ 1 in (9), we get We consider the choice g 2iþ1 ¼ ÀH l i c 2iþ1 with d T 2iþ1 c 2iþ1 ¼ 1; which clearly holds in (10). Now, we define Later, as Theorem 2.3, we will conclude that the above system has solution and H l i is well defined. Therefore, the updating formula for H i is given by the following: where w 2iþ1 can be any vector satisfying (11). Now, to satisfy (12) and complete the induction, H 1 should be chosen, so that H 1 a 1 ¼ H 1 a 2 6 ¼ 0 or Let H 0 be an arbitrary nonsingular matrix. We obtain H 1 from H 0 using a rank one update. Take and we choose w 1 2 R n satisfying the next condition: Clearly, (14) can be held with a proper choice of w 1 2 R n , whenever a 1 and a 2 are linearly independent. Thus, we have a rank one update as follows: where w 1 is an arbitrary vector satisfying (14). To compute the general solution and update the second phase for the ith iteration, we introduce a matrix H l i with properties H l i c j ¼ H l i a j ¼ 0; j ¼ 1; . . .; 2i: Therefore, we define the matrix H l i by a rank one update as the next formula: Notice that w 2i 2 R n is an arbitrary vector satisfying the following condition: Hence, the general solution of the ith iteration is given by where s 2 R n is arbitrary. Clearly, the general solution for the last iteration is presented by the following: where s 2 R n is arbitrary. Therefore, we proved the following theorem.
Theorem 2.2 Assume that we have m ¼ 2l arbitrary linearly independent vectors a 1 ; . . .; a m 2 R n and an arbitrary nonsingular matrix H 0 2 R nÂn . Let H 1 be generated by (15) and w 1 is satisfying (14), and the sequence of matrices H i , i ¼ 2; . . .; l, be generated by the following: with w 2iÀ1 2 R n satisfying the following condition: In addition, let the sequence of matrices H l 1 ; . . .; H l l be generated by (16) with w 2i 2 R n satisfying (17). Then, when we are at the ith iteration, the following properties (1)-(4) hold for i ¼ 1; . . .; l.
we are at the ith iteration. matrices H i be generated by (19) with w 2iÀ1 2 R n satisfying (20). Then, for all i, 1 i l, and j, 2i j m, the vectors H i a j are nonzero and linearly independent.
Proof We proceed by induction.
. . .; a m are nonzero and linearly independent and H 0 is nonsingular, H 0 a j for 1 j m are nonzero and linearly independent.
Therefore, the vectors H 1 a j , for 2 j m, are nonzero and linearly independent. Now, we assume that the theorem is true up to 1 i l À 1, and then, we prove it for i þ 1. From (12), We need to show that the relation implies that a j ¼ 0, for 2i þ 2 j m. Using (21), we can write (22) as follows: By the induction hypothesis, the vectors H i a j , for 2i j m, are nonzero and linearly independent. Now, as the assumption of the induction, we are going to prove that the vectors H l i a j are nonzero and linearly independent for 2i þ 1 j m. Using (16), we have We must prove the relation implies that a 0 j ¼ 0, for j ! 2i þ 1. Using (24), we can write (25) as follows: As the linearly independence of H i a j , for 2i j m, we conclude that a 0 (i) When we are at the ith iteration, we have (14) and (20) has solution.
(It is noted again that when we choose j ¼ 2i, we are at the ith step. In fact for the ith iteration, we have two choices j ¼ 2i À 1 and j ¼ 2i.) Theorem 2.5 For the matrices H i generated by (15) and (19) and the matrices H l i given by (16), we have The notations R and N stand for the range and nullspace, respectively.
Then, the matrix L i defined by L i ¼ ðA i Þ T P i is lower triangular and nonsingular.
Proof As ðL i Þ j;k ¼ a T jpk , Theorem 2.2 implies that ðL i Þ j;k is zero for j\k, and thus, L i is lower triangular. As Theorem 2.3 and Corollary 2.4, we have a T Remark 2.8 For the matrices H i generated by (15) and (19), we have . . .; lÞ; where x i is a solution of the first 2i equations of the system and k i will be discussed in Algorithms 1 and 2. The matrix H 0 can be any arbitrary matrix, but to reduce computational complexity, and for stability reasons, we choose H 0 unitary matrix I n . Numerical stability and computational complexity will be discussed in Sects. 4 and 5, respectively.
. . .; a m Ã T is a full row rank linear matrix, where m ¼ 2l and Ax ¼ b is a compatible full row rank linear system, x 2 R n and b 2 R m .
1. Let x 0 2 R n be an arbitrary vector and choose H 0 ¼ & (c) If r 1 ðx 0 Þr 2 ðx 0 Þ ¼ 0 and one of the residual values is nonzero, without loss of generality, we assume that r 2 ðx 0 Þ 6 ¼ 0 and we let & (d) If r 1 ðx 0 Þr 2 ðx 0 Þ ¼ 0, and both of the residual values are zero, x 1 will be x 0 and go to (5).
9. Stop (x l is a solution of the system). From (18), we can compute general solution of the system after the final iterate by x l l ¼ x l À H T l l s where s 2 R n is arbitrary.

Remark 2.9
To reduce computational complexity of Algorithm 1, we propose the following tactics: 1. Taking x 0 is proper in step (1) as Algorithm 1.
2. To compute w 1 2 R n , we take t j ¼ H 0 c 1 , and then, we define w 1 as follows: Then, we take t j M ¼ maxfjt j j : j 2 f1; . . .; ng; suchthatw T 1 t j ¼ 1g: To determine the other w 2iÀ1 2 R n , i ¼ 2; . . .; l, we let t 0 j ¼ H l iÀ1 c 2iÀ1 , and then, we continue by the similar way. 3. We can choose w 2i ¼ z i 2 R n , i ¼ 1; . . .; l, by defining the next parameters d j ¼ H i a 2i , and Now, we take d j M ¼ maxfjd j j : j 2 f1; . . .; ng; suchthat w T 2i d j ¼ z T i d j ¼ 1g: Now, assume that the vectors a 1 ; . . .; a m are linearly independent. According to Theorem 2.5, we have NðH i Þ ¼ 2i À 1 and NðH l i Þ ¼ 2i. Therefore, 2i À 1 rows of matrix H i are depend on the other rows of H i and 2i rows of the matrix H l i are depend on the other rows of H l i . As we have defined the matrix H i and H l i in such a way that exactly 2i À 1 rows of H i and 2i rows of matrix H l i are the zero vectors, we can delete these zero vectors using appropriate operator. Next, we will show how to reduced the number of rows of the Abaffian matrix by two in every iterate and economize the space needed for our new model.

New compression two-step ABS algorithm
Algorithm 2 Assume that A mÂn ¼ ½a 1 ; a 2 ; . . .; a m T is a full row rank linear matrix, where m ¼ 2l and Ax ¼ b is a compatible full row rank linear system, x 2 R n and b 2 R m .
1. Let x 0 2 R n be an arbitrary vector and take H 0 ¼ I n . Set & (c) If r 1 ðx 0 Þr 2 ðx 0 Þ ¼ 0 and one of the residual values is nonzero, without loss of generality, we assume that r 2 ðx 0 Þ 6 ¼ 0 and we let & (d) If r 1 ðx 0 Þr 2 ðx 0 Þ ¼ 0, and both the residual values are zero, x 1 will be x 0 and go to (5).
9. Stop (x l is a solution of the system). General solution of the system is obtained after the final iterate by x l l ¼ x l À H T l l s, where s 2 R n is arbitrary.

Remark 3.1
We can present the stepsize for Algorithms 1 and 2 as the following form: Notice that a 2iÀ1 ; a 2i ; b 2iÀ1 , b 2i are updated in steps 2 and 6 andr i is defined the updated residual value of the ith step. Thus, we have the exact values presented in 3(d) and 7(d) for the stepsize. Notice that dimension of H i and z i is compressed in Algorithm 2.
Remark 3.2 To reduce computational complexity of Algorithm 2, we propose the following tactics.
1. Taking x 0 is proper in step (1) as our new algorithm. 2. To compute w 1 2 R n , we take t j ¼ H 0 c 1 , and then, we define w 1 as follows: Then, we take t j M ¼ maxfjt j j : j 2 f1; . . .; ng; suchthatw T 1 t j ¼ 1g: To compute the other w 2iÀ1 2 R nÀð2iÀ2Þ , i ¼ 2; . . .; l, we let t 0 j ¼ H l iÀ1 c 2iÀ1 , and then, we continue the similar way. 3. We can choose w 2i ¼ z i 2 R nÀð2iÀ1Þ , i ¼ 1; . . .; l, by defining the next parameters d j ¼ H i a 2i , and Now, we take d j M ¼ maxfjd j j : j 2 f1; . . .; n À ð2iÀ 1Þg; such that w T 2i d j ¼ z T i d j ¼ 1g: Notice that by reducing the numbers of rows of the Abaffian matrix by two, for per step, we choose the other parameters as proper dimensions.

Analysis of error propagation in the new versions of two-step ABS models
We investigate the stability of Algorithm 1, for nonsingular system (2), by the backward error analysis technique due to Broyden and Galantai as [8,9,12]. Some basic results given by Broyden [9] and Galantai [12] are extended here. Galantai's method has important role in our work. Systems (1) and (2) are equivalent, but for proceeding with Galantai's method, we choose system (2) in all of this section. We prove the analysis of error propagation for nonsingular systems, and then, we will conclude that it is valid for every full row rank system. Galantai studied the numerical stability of the ABS class for linear systems of the form: n where the matrix A T is nonsingular. P ¼ ½p 1 ; . . .; p n and V ¼ ½v 1 ; . . .; v n are n Â n-type nonsingular matrices with column vectors p j and v j ðj ¼ 1; . . .; nÞ. Denote by I the unit matrix of n Â n type. The ABS class has the following form:

Algorithm 3
Let x 0 2 R n be arbitrary For k ¼ 1; . . .; n; Compute: where the ABS update algorithm is given by Set H 1 ¼ I; For k ¼ 1; . . .; n; Compute: Algorithm 3 is finitely terminated in n steps. The matrix V is scaling the system A T x ¼ b. As [20], the pair (P, V) is said to be the A T -conjugate if V T A T P ¼ L is the lower triangular. The ABS algorithm generates all A T -conjugate directions for suitable choices of parameters [2]. Therefore, the ABS class of methods coincides with those studied by Stewart [20]. Results on the numerical stability of conjugate direction methods are given in [9,21]. A stability analysis for decent methods is given in [7]. See also [3,4,10]. As Corollary 2.7, we conclude that the pair (P, I) is A T -conjugate for Algorithm 1, considering nonsingular matrix of system (2). Using the basic idea of Broyden's backward error analysis method, we present Algorithm 1 as the following form: where X k is the solution of the problem. Assume that an error n j occurs at the jth step and that is error propagates further. It is also assumed that no other source of error occurs. The exact solution X n 2 is given by the following: while the perturbed solution X 0 n 2 is given by the following: If the quantity jjX n 2 À X 0 n 2 jj is large, then the algorithm (26) must be very unstable. As Broyden ideas, it must be small for stable algorithms. Consequently, jjX n 2 À X 0 n 2 jj is a measure of stability for our new version of algorithm. Now, we recall some of the important parameters of Algorithm 1 for the system A T x ¼ b, A 2 R nÂn , where the matrix A T is nonsingular considering system (2). Let Math Sci (2017) 11:333-343 339 P ¼ ½p 1 ; . . .;pn 2 , with column vectors defined from Theorem 2.6 as the formp i ¼ ðp 0 i ; p i Þ. We know that, p i ¼ H T i z i is obtained the first-phase rank reducing process of the ith iteration. Thus, P is n Â n-type nonsingular matrix. Here, we take V ¼ I n ¼ ½v 1 ; . . .; vk 2 with v k ¼ e k ¼ ðe 2kÀ1 ; e 2k Þ, ðk ¼ 1; . . .; n 2 Þ. Therefore, V ¼ I n ¼ ½e 1 ; . . .; e n is an n Â n identify matrix. As Corollary 2.4, we have p T k Ae 2kÀ1 ¼ p T k Ae 2k 6 ¼ 0. Taking v k for the stepsize means that we can arbitrary choose one of the values of e 2kÀ1 or e 2k . Also from Remark 3.1, we have We choose H 0 ¼ I n and we present k k ; x k ; p k , and H k as the following form: x k ¼ x kÀ1 À k k p k ; ð30Þ Updating vectors and residual values can be computed as Algorithm 1, for k ¼ 1; . . .; n 2 . Next, we investigate the stability of the conjugate direction methods of the form (29) and (30) with a special emphasize on the ABS update (31) and (32). It is noted again that for (29)-(32) and the others related parameters, we choose V the unitary matrix and the pair (P, I) is A T conjugate for a nonsingular matrix as system (2). Considering these points, Algorithm 1 coincides with those studies in [12] for a unitary V and n 2 steps. Some basic results of [9,12,20] are extended here. Let us introduce the following notation: It is easily verified that P 2 k ¼ P k ; i:e:; P k is a projector. The matrix P k is a rank one projector onto the space spanned by p k along the orthogonal complement of the space spanned by Av k . Therefore, we have RðP k Þ ¼ Rðp k Þ and Nðp k Þ ¼ R ? ðAv k Þ. Note that the projectors P 1 ; . . .; P n are called conjugate by Stewart [20], if we have P t P j ¼ 0ðt\jÞ. The following lemma is immediate consequent of Stewart's definition [20]. Lemma 4.1 Let P 1 ; P 2 ; . . . be conjugate projectors. for i k we have: P i ðI À P k ÞðI À P kÀ1 Þ Á Á Á ðI À P 1 Þ ¼ 0Á ð 34Þ Proof By conjugacy, for i\j, we have and since P i is a projector, Together (35) and (36), we have (34). h With the notation P k as (33), the method of the class (29) and (30) has the following form (motivating techniques in [12,20] and extending them for our model): Denote by x Ã the solution of the linear system. By taking E k ¼ x Ã À x k , we have the recursion E k ¼ ðI À P k ÞE kÀ1 with the following solution: Now, we introduce the notation ðI À P k Þ Á Á Á ðI À P j Þ ¼ Q k;j ; for k ! j. Now, let us suppose that an error occurs at the ðk À 1Þth step and only at it. The perturbed results of the ðk À 1Þth step are denoted by x 0 kÀ1 . The perturbed results of further steps are denoted by x 0 j ðj ¼ k; . . .; n 2 Þ. Then, we have from which it follows that the error occurring in the final iteration: The matrix Q n 2 ;k can be considered as the error matrix. Hence, we have the error bound jjx n 2 À x 0 n 2 jj jjQ n 2 ;k jjjjx kÀ1 À x 0 kÀ1 jj: ð37Þ A method of the classes (29) and (30) is considered to be optimal in the sense of Broyden [9], if jjQ n 2 ;k jj is minimal for all k. First, we characterize Q n 2 ;k . Using the A T -conjugate property, we have P t P j ¼ 0ðt\jÞ. Therefore, P t Q n 2 ;k ¼ 0ðk t n 2 Þ and ðI À P t ÞQ n 2 ;k ¼ Q n 2 ;k ðk t n 2 Þ, which means that Q n 2 ;k Q n 2 ;k ¼ Q n 2 ;k , i.e., Q n 2 ;k is a projector. By observing that RðI À P k Þ ¼ NðP k Þ and NðI À P k Þ ¼ RðP k Þ, we conclude that  (1) and (2) are the equivalent. Moreover, the points 2 and 3, as Remarks 2.9 and 3.2, are useful to increase stability of our algorithms.

Computational complexity and numerical results
At first, considering Remark 2.9, we compute the number of multiplications, for Algorithm 1, but for arbitrary nonsingular matrix H 0 , as follows (Table 1): Hence, the total number of multiplications for the l iterates is þ ð2n þ 1Þðn À 2i þ 1Þ þ 2n; N ¼ 2mn 2 À m 2 n þ OðmnÞ À Oðm 2 Þ þ OðmÞ: The total number for Huang's method is 3 2 mn 2 þ OðmnÞ multiplications and 3mn 2 À 7 4 m 2 n þ 1 6 m 3 þ OðmnÞ þ Oðm 2 Þ þ Oðn 2 Þ for the two-step methods presented by Amini et al. [5,6]. Our computations indicate that our numerical results are better than Amini et al.'s two-step algorithms. In addition, when n\2m, we need cheaper number of multiplications up to the corresponding Huang's algorithm. In addition, when m tends to n, the total number for our two-step method is n 3 multiplications, while it is 3 2 n 3 for Huang's algorithm and 17 12 n 3 for those two-step methods given by Amini et al. [5,6]. Obviously, when m and n are not too large, the lower order terms have influence on our results.
Remark 5.1 Now, considering Algorithm 1 with Remark 2.9, setting of H 0 the identity matrix I n , we have the better results. In this case, for the ith iteration, the number of multiplications of the matrix H i is ðn À 2i þ 2Þð4i À 1Þ and it is ðn À 2i þ 1Þð4i þ 1Þ multiplications for the matrix H l i . Therefore, our algorithm needs nm 2 À 2 3 m 3 þ OðnmÞ þ Oðm 2 Þ multiplications. Hence, for m ¼ n, 1 3 n 3 , multiplications plus lower order terms are needed. Notice that the main storage requirement is the storage of H l i , which has at most 1 4 n 2 positions. This is half of the storage needed by the Gaussian elimination method. Our new version ABS algorithm is computationally better than the classical Gaussian elimination method, having the same arithmetic cost but using less memory and no pivoting is necessary.

Remark 5.2
The operator D is used to compress and economize our required space, via deleting the zero rows of the Abaffian matrix by two for per iteration. Therefore, it has not any effect on the number of multiplications, and Algorithms 1 and 2 have the same numbers of multiplications. Thus, both of our works are better than Amini et al.'s two-step ABS algorithms and Huang's method. Furthermore, we have good competition with the classical Gaussian elimination method.

Conclusion
This paper presents two new version of ABS models for the general solution of compatible full row rank linear systems of equations in at most mþ1 2 Â Ã iterates. For the second algorithm, the number of rows of the Abaffian matrix is reduced by two for every iteration and the space is compressed as much as possible. A theory of conjugate projectors is developed and used for the analysis of error propagation. We evaluate the numerical stability of our schemes using Broyden's backward error analysis method and Galantai's technique. In addition, both our models need less computational complexity up to those corresponding two-step ABS algorithms and Huang's method. Our new two-step ABS algorithms is computationally better than the classical Gaussian elimination method, having the same arithmetic cost but using less memory and no pivoting is necessary.