On the incomplete oblique projections method for solving box constrained least squares problems

The aim of this paper is to extend the applicability of the incomplete oblique projections method (IOP) previously introduced by the authors for solving inconsistent linear systems to the box constrained case. The new algorithm employs incomplete projections onto the set of solutions of the augmented system Ax − r = b, together with the box constraints, based on a scheme similar to the one of IOP, adding the conditions for accepting an approximate solution in the box. The theoretical properties of the new algorithm are analyzed, and numerical experiences are presented comparing its performance with some well-known methods.

computational mechanics and optimization problems. In practice, those systems are often inconsistent, and one usually seeks a point x * ∈ n , l i ≤ x * i ≤ u i , i = 1, . . . , n, that minimizes a certain proximity function. In this paper, for such possibly inconsistent systems we will consider the standard problem where B = {x : x ∈ n , l i ≤ x i ≤ u i , i=1,. . . , n}, A ∈ m×n , b ∈ m , m ≥ n. . D m denotes the norm induced by the positive definite diagonal matrix D m ∈ m×m .
In [19] we have introduced, for inconsistent problems Ax = b, the IOP algorithm that converges to a weighted least squares solution of that system. That algorithm uses an incomplete oblique projections scheme onto the solution set of the augmented system Ax − r = b, in order to solve min x∈ n Ax − b 2 D m . We proved that this problem is equivalent to the problem [19]: min{ p − q 2 D : for all p ∈ P and q ∈ Q}, (2) being P = {p : p = [x; r] ∈ n+m , x ∈ n , r ∈ m , Ax − r = b}, and Q = {q : q = [x; 0] ∈ n+m , x ∈ n , 0 ∈ m }. D is a diagonal matrix of order n + m, whose n first elements are 1's, and the last m coincide with those of D m . That result led us to develop the IOP method for solving (2), applying an alternate projections scheme between the sets P and Q, similar to the one of Csiszár and Tusnády [10], but replacing the computation of the exact projections onto P by suitable incomplete or approximate projections. In [21] the results of [19] were generalized for solving rank deficient problems.
In this paper we add to the IOP alternating algorithm the condition x ∈ B; for that purpose we define the sets: adopting the distance d(p, q) = p − q D , for all p ∈ P b , q ∈ Q b . We consider the following basic scheme:

Algorithm 1 (Basic Alternating Scheme)
Given feasible vectors p k = x k ; r k ∈ P b , q k = x k ; 0 ∈ Q b . Let p = [x; r],p = x;r . Definep Define p k+1 = x k+1 ; r k+1 = x;r , and q k+1 = [x k+1 ; 0] In each step of Algorithm 1 a minimum norm problem with linear constraints (4) needs to be solved. For that purpose we use the HLWB algorithm [9]. This inner iteration needs of course to be terminated at some step. Therefore in practice we must use instead ofp some approximation of it.
In the following sections we present the extended BI OP algorithm based on the basic scheme of Algorithm 1, adding the conditions for accepting an approximate solution of (4). The theoretical properties of the new algorithm are analyzed in Section 3, and numerical experiments are presented in Section 4 comparing its performance with some known methods.

Notations
• . denotes the Euclidean norm, and . D the norm induced by a positive definite matrix D.
• I n denotes the identity matrix in n×n , and e i , the i-th column of I n .
• The upper index T stands for the transpose of a matrix.
• D = diag(d), denotes a diagonal matrix of order n, whose diagonal is d = (d 1 , . . . , d n ). • P q k will denote the exact solutionp = x,r of problem (4), where q k = x k ;0 .

Some basic results
The convex feasibility problem (CFP) is to find a point (any point) in the nonempty intersection C = m i=1 C i = ∅, of a family of closed convex subsets C i ⊂ n , i = 1, 2, . . . , m of the n'-dimensional Euclidean space.
Projection algorithms employ projections onto the individual convex sets in order to reach the required point in the intersection. They employ projections onto the individual sets in various ways [5-7, 12, 18]. In the sequel we denote by P i the orthogonal projection onto C i , and P D i the oblique projection onto C i , where D is a positive definite matrix of n ×n .
The best approximation problem (BAP) is to find the projection of a given point y ∈ n onto the nonempty intersection C. Thus, it seeks a point in the intersection of the convex sets which is closest to the point y.
Aiming at clarifying the applicability of HLWB within the new approach, it is convenient to point out that, it solves the BAP problem [8,9]. It computes the simultaneous projections onto each set, takes a convex combination of the intermediate points, and defines the next iterate.
Next, the simultaneous HLWB algorithm [8] is presented for the sake of completeness.

Algorithm 2
Initialization: Let y 0 = y be the given point y whose projection The sequence {σ j } must satisfy the following conditions: In the new algorithm, HLWB will be used for obtaining an approximation to the solution of problem (4), where the former convex set C will be the intersection of . . , n, and the initial point at the k-th iteration of Algorithm 1, will be y 0 = x k ; 0 .
For each constraint of Az − r = b, we will denote by Given an arbitrary y = [z; r] ∈ n+m , we denote for each i = 1, . . . , m, , e i ∈ m , and the oblique projection of y onto C i , by We define Given an arbitrary y = [z; r] ∈ n+m , for each i = 1, . . . n we denote r m+i (y) = max(0, max(z i − u i , l i − z i )), and P D m+i (y) is calculated by In particular, in this paper we will use a diagonal weighting matrix D = I n 0 0 D m , therefore the last previous projections will be computed by Analogously, the projections (5) onto C i , i = 1, . . . , m, will be computed by , where e i ∈ m .

Incomplete oblique projections algorithm
As said in the Introduction, for possibly inconsistent systems Ax = b, x ∈ B, A ∈ m×n , b ∈ m , m ≥ n, we will consider the standard problem (1). We will follow the alternating scheme of Algorithm 1 for solving the equivalent problem that minimizes the distance between the sets described in (3). For that purpose we will use a diagonal weighting matrix D ∈ (n+m)×(n+m) , as defined in the previous Section.
Given p k = x k ; r k ∈ P b , and its projection q k = x k ; 0 onto Q b , we will denote by P q k the projection of q k onto P b , which is the solution of problem (4).
As mentioned before, instead of defining p k+1 as the exact solution of (4), we define it by means of the incomplete resolution of that problem.
In order to define an inexact projection, we consider the following: Aiming at obtaining properties of the sequence p k generated by the new algorithm that guarantees convergence to the solution of (1) we establish an "acceptance condition" that an approximation to P q k must satisfy, using the iterative HLWB algorithm for solving (4).
Lemma 1 Applying the HLWB algorithm for solving (4) it is possible to find an iterate p j = [z j , v j ] satisfying conditions (8) and (9).
Proof Since the sequence obtained by HLWB converges to the solution of (4), the sequence {p j } = {[z j ; v j ]} tends to P q k = ẑ;v , then P r (p j ) − P q k 2 D goes to zero when σ j tends to zero.
Then, given 0 < β k < 1, as σ j tends to zero, a σ j exists such that σ j < β k . Hence, it is possible to satisfy condition (8).
If P q k = p k , P q k − q k 2 D < p k − q k 2 D , and considering that p j − P q k 2 D tends to zero, there exists an j such that p j − q k 2 D < p k − q k 2 D . Furthermore, considering that the sequences P r (p j ) − P q k 2 D and P r (p j ) − p j 2 D go to zero, it is possible to satisfy condition (9).
Our proposal is to replace in the new algorithm the exact solution P q k of (4) by P r (p j ) if it satisfies the above conditions. Thus, we will define (8) and (9).
Notation In the following, the index j corresponding to the accepted approximation p j = [z j , v j ] will be denoted by j k to indicate that it corresponds to the k-th iteration of the new algorithm.

Algorithm 3 Bounded Incomplete Oblique Projections (BIOP)
Initialization: Given γ , 0 < γ < 1, a positive definite diagonal matrix D m of order m and β k ∞ 0 , such that β k > 0 and β k tends to zero, if k tends to ∞.
• Calculate an approximation to the solution of (4) satisfying the conditions of definition 2, applying HLWB from the initial point y 0 = [z 0 ; v 0 ] = q k , and using the sequence {σ j }: , until finding an index j k such that y j k = z j k , v j k and P r y j k satisfy conditions (8) and (9).

Convergence of the BIOP algorithm
We will consider the set of solutions to the problem (1): The corresponding set in m+n : In the following we analyze the convergence properties of the sequence given by Algorithm 3 studying the relationship between two consecutive iterates p k and p k+1 .

Lemma 2
Let p k = x k ; r k be the sequence generated by the Algorithm 3, then (ii) The sequence r k D m is decreasing and bounded, therefore it converges.
(iii) The following three sequences tend to zero: where P q k is the exact solution of (4). (iv) The sequence A T D m r k −μ k goes to zero, μ k being the vector of the Karush-Kuhn Tucker (KKT) multipliers associated to the bound constraints at the solution of the problem (4).
Proof We get, as a consequence of the definition of p k+1 = x k+1 ; r k+1 , that p k+1 = P r z j k ; v j k satisfies the condition (9), is decreasing and bounded, and therefore it converges. Hence, x k+1 − x k 2 tends to zero. As a consequence of that, we also obtain that p k+1 − p k 2 D goes to zero.
By the conditions required in Definition 2 to accept p j k and P r ( Also, from the conditions in Definition 2 to accept p j k = [z j k ; v j k ], we obtain that Thus, the sequence of parameters {σ j k } converges to zero because {β k } tends to zero. Therefore, from the properties of the HLWB algorithm [9], p j k − P q k D tends to zero, when k tends to ∞. Hence, we also obtain that p k+1 − P q k 2 D tends to zero. Then, also p k − P q k 2 D tends to zero. In order to prove (iv), the components of the exact solution P q k will be denoted by x k ;r k , to exhibit their dependence of the k-th iterative step.
The solution of a convex problem with linear equalities and inequalities constraints satisfies KKT conditions [5]. Therefore, the exact solution P q k = x k ;r k of problem (4), satisfiesx where μ k is the vector of multipliers associated to the bound constraints at P q k . It is known that, if l i <x k i < u i then μ k we get that A T D m r k − μ k goes to zero. Thus, (iv) follows. Furthermore, as a consequence of the previous result it turns out that μ k is bounded, and considering that x k −x k tends to zero, we get that μ k T x k −x k also tends to zero.
In order to prove (v), it is known that for each i ∈ {1, 2, . . . , n} such that l i <x k i < u i , the product μ k i min Hence, since x k i −x k i converges to zero, an index k u exists such that for k > k u the product is equal to μ k i x k i − x k i . Hence, we get (v).

Lemma 3
Let {p k } = x k ; r k be the sequence generated by the Algorithm 3, then Proof Considering that x k ∈ B, we know a subsequence {x k s } exists convergent to x in B. Hence, the subsequence {p k s } = {[x k s ; r k s ]}, satisfying r k s = Ax k s − b, is convergent to [x; r] ∈ P b . By (iv) of the Lemma 2, we know that A T D m r k s − μ k s tends to zero. Furthermore, by (iii) of the Lemma 2 x k s −x k s and r k s −r k s D m tend to zero, then x k s − x tends to zero, as well as r k s − r D m .
Therefore, μ k s = x k s −x k s + A T D m r k s converges to μ = A T D m (r), satisfying r = Ax − b. Moreover, by the KKT conditions of problem (4), Hence, for k s > k i , μ k s i = 0, so we get that μ i = 0. If x i , satisfies x i = l i , since |x k s i − x i | tends to zero, the same happens with |x k s i −l i |, and therefore there exists k L i , such that for k s > k L i , |x k s i −l i | < 0.5 * |u i −l i |, then μ k s i ≥ 0, and it follows that μ i ≥ 0. Likewise we deduce that if x i = u i , then μ i ≤ 0.
Hence, μ i min (x i − l i , u i − x i ) = 0, with μ i ≥ 0 if x i − l i = 0, and μ i = 0 if l i < x i < u i , and μ i ≤ 0 when x i = u i , taking into account thatx k s tends to x. As a consequence of the last result, [x; r] ∈ P b , satisfies [x; r] ∈ S D Therefore, [x; r] satisfies the optimality conditions of the problem (1). Then, (i) follows.
From the proof of (i), we get (ii). Every convergent subsequence has its limit point in S D .
In order to prove (iii), under the hypothesis that the set S D has only one element, we consider that the domain is bounded, and since the sequence x k+1 − x k goes to zero, applying Theorem 14.1.3 in J.M. Ortega and W. C. Rheinboldt [15], we obtain that the sequence x k ; r k converges to p * ∈ S D .

Numerical experiments
The objectives of the following experiments are two-fold. First we compare our algorithm with other methods, in relation to the rate of decrease of the norm of the residual, reporting the number of iterations needed for satisfying the convergence conditions and the corresponding CPU time. In the second part we analyze the BI OP behavior when solving image reconstruction problems comparing it with IOP (Scolnik et al. [19]). All visualizations have been obtained using MATLAB 7.0 For the first purpose, we have implemented our algorithm in Fortran 77 and MATLAB 7. We made some comparisons of the BI OP algorithm with BVLS, implemented in Fortran 77, which solves linear least-squares problems with upper and lower bounds on the variables, using an active set strategy. It is documented by Stark and Parker [23]. We have also compared our algorithm with one method of the SLS system, a software in MATLAB for sparse least squares problems. That method is SBLS2, for solving sparse box constrained least squares problems, in the version of Mikael Lundquist, Linkoping University, Sweden, http://www.mai.liu.se/ ∼ milun/sls, which computes the solution by a version of the block principal pivoting algorithm described by Portugal et al. in [17].

Test problems
The problems used for these numerical experiments are some from the RRA (real, rectangular, assembled) part of Harwell-Boeing Collection. Those least squares problems are from the set LSQ. In particular we used four matrices, such that the second and fourth matrix have the same pattern as the first and third respectively but are much more ill-conditioned. The matrices in this set are: WELL1033 and ILLC1033(real unsymmetric, 1033 by 320, 4732 entries), WELL1850 and ILLC1850 (real unsymmetric, 1850 by 712, 8758 entries). Also, aiming at testing problems with larger dimensions, we tested some randomly generated dense matrices with DQRT15(LAPACK) which, using the parameter RKSEL = 2, generates a rank deficient matrix A(mxn) such that m is the number of rows and n is the number of columns of A, and rank = 3 4 min(m, n). Other systems arise from the two image reconstruction problems used by Popa and Zdunek in [16] and other problems from the SNARK System [2]. The first simulate real objects in electro-magnetic geotomography, leading to problems whose data comes from projections made with a limited angular range. Those problems, A1 and A2, lead to inconsistent systems, and the corresponding matrix has deficient rank due to the angle limitations of the projections. Also aiming at testing the algorithm with larger dimensions, we considered other problems from the SNARK System, with m ≥ n. They are B1, B4, B6 and B7. The dimensions of these matrices are showed in Table 1. For A1 and A2 we analyze the results with a system Ax = b + δb, arising from simulating noisy perturbations of the right hand side b = Ax exact . The problems here analyzed exactly correspond to those mentioned as case 2 in [16]. Starting from the knowledge of b = Ax exact , a perturbation δb is defined satisfying δb / b ≈ 5.5 %. Since δb = δb A +δb A ⊥ , where δb A ∈ R(A) and δb A ⊥ ∈ R(A) ⊥ , where R(A) denotes the subspace spanned by the columns of A, and R(A) ⊥ stands for the orthogonal subspace to R(A). Those perturbations are applied to each problem according to:  Perturbations for both problems were computed using a standard procedure in MATLAB.

Numerical results
In these numerical experiments we have considered in the implementation of BI OP the matrix D m = I m , and the sequence {σ j } = 1 j +1 , for j > 0 and parameters w i , equal to 1 n+m , for i = 1, 2, . . . , n + m. BIOP was implemented sequentially, with parameters γ = 10 −1 and β k = 1 k+1 . The stopping condition for the three algorithms was: if | r k+1 − r k | < max( r 0 , 1), with = 10 −6 . The convergence criterion allows to terminate the iterations when decreasing stagnates.
In Table 2 we compare the results obtained with BIOP and BVLS, implemented in Fortran, in regard to the number of iterations(Iter) and the CPU time required by them for each problem p. R E (p) denotes the final value of the norm of the residual obtained by the algorithms for each problem p, when using = 10 −6 for the stopping criterion.
We report for BIOP the number of iterations(Iter) and CPU time required for each problem p, for reaching the stopping condition using = 10 −6 . The number of iterations (Iter) of BIOP is the total of inner iterations required by the process since in each one of the main iterations another inner iterative process is needed. We report for BVLS the number of iterations(Iter) and CPU time required for each problem p, for reaching the stopping condition or a norm of the residual less or equal than the one obtained by BIOP.
We also give, the corresponding CPU time in seconds. The dimensions of the test problems are shown using m for the number of rows, n for the number of columns.
It is necessary to point out that the solver BVLS does not deal with sparsity. Thus, as observed in the previous Table 2, its CPU times are quite significant. Taking this into account it is fair to compare our algorithm with another able to solve sparse problems under the same computational framework.
In Table 3 we compare the number of iterations and the CPU time for BIOP and for the SBLS2, both implemented in MATLAB. We report the results using the same notation of the previous table. It follows from those results that the SBLS2 is unable to solve several problems. The reason is that matrices A1, A2 are rank deficient. On the other hand, BIOP solved all test problems, a convenient property of the projection methods. Moreover, when both algorithms solve a given problem BIOP is faster than SBLS2.
In Fig. 1, we analyze the performance data using the performance profiles of Dolan and Moré, as described in [11]. We compare the performance of BIOP and SBLS2 methods in relation to the CPU time required for each problem when the convergence criterion is met. The percentage of the test problems for which a given method is best, is given on the vertical axis of the plot.
It follows from those results that the BIOP method is also able to reduce the residual more rapidly.

Image reconstruction problems
In our previous papers we have compared IOP with other algorithms for solving image reconstruction problems. For instance, in [19] we compared IOP with ART (underrelaxed) [5], which is one of the most important class of methods in this field, and CAV [6] using test problems from SNARK [2], showing that is very efficient. Those results showed that IOP requires more CPU time than ART in the sequential implementations, but it obtains better images in some of the problems. In regard to CAV our algorithms are superior both from the points of view of the CPU time and the images' qualities as shown in [19]. In [21] we have also compared a version of IOP called EIOP for dealing with rank deficient problems with KERP (Popa and Zdunek [16]) and LANDW [14]. The numerical results given there demonstrate that EIOP outperforms KERP and LANDW, both in the rate of convergence and the required CPU time. In the case of image reconstruction problems, EIOP obtains the closest image to the original one in much less CPU time than KERP. In that paper we used problems A1, A2 with the perturbations proposed in [16], and it turned out that EIOP gave better results than KERP both in terms of CPU time and distance to the true image.
In all those papers we have considered unconstrained linear least squares problems. Now, in this paper we use BIOP with bounds on the variables, for testing if it is possible to improve in image reconstruction problems the final distance to the true solution when using x ≥ 0, and suitable upper bounds. Hence, in this part we analyze the BI OP behavior when solving image reconstruction problems compared with IOP. From the theoretical results we observe that IOP requires γ ≤ 0.5 for convergence in the case of rank deficient problems [21]. We have used γ = 0.01 in the initial iteration, then γ = 10 −1 . In Fig. 2 we compare the performances of BIOP and IOP by means of the distance of the reconstructed images for problem A1(case2) and A2(case2), from [16]. In Fig. 3 we compare the performances of BIOP and IOP by means of the distance of the reconstructed images for problem B1, B4, B6 and B7, from SNARK.
From the point of view of the number of iterations BIOP reaches faster than IOP the image closest to the original for all test problems. BIOP also obtains the minimum distance, and keeps it below the one corresponding to IOP. It reaches distances similar to those obtained with regularization techniques [20,22].

Conclusions
We have presented in this paper the results obtained when including bound constraints into the IOP algorithm. The need of including these sort of constraints appear in, for instance, image reconstruction applications where the solution is positive and the density values are bounded. Also in nonlinear optimization problems, like in the SQP method, it is necessary to solve least squares problems in a box.
As it can be seen in the graphs A1 − A2 and B1 − B7 when x ≥ 0 is used, the convergence to the true image is faster, from the point of view of the number of iterations, than in the case of using IOP. In previous papers we had to use regularization techniques in order to obtain results similar to the ones here presented, but at a higher computational cost [20,22].
As the numerical results show, the BIOP effectiveness is remarkable in several problems as compared to SBLS2 because in particular it does not fail to converge.