Updating QR factorization procedure for solution of linear least squares problem with equality constraints

In this article, we present a QR updating procedure as a solution approach for linear least squares problem with equality constraints. We reduce the constrained problem to unconstrained linear least squares and partition it into a small subproblem. The QR factorization of the subproblem is calculated and then we apply updating techniques to its upper triangular factor R to obtain its solution. We carry out the error analysis of the proposed algorithm to show that it is backward stable. We also illustrate the implementation and accuracy of the proposed algorithm by providing some numerical experiments with particular emphasis on dense problems.


Introduction
We consider the linear least squares problem with equality constraints (LSE) where A ∈ R m×n , b ∈ R m , B ∈ R p×n , d ∈ R p , x ∈ R n with m + p ≥ n ≥ p and ·  denotes the Euclidean norm. It arises in important applications of science and engineering such as in beam-forming in signal processing, curve fitting, solutions of inequality constrained least squares problems, penalty function methods in nonlinear optimization, electromagnetic data processing and The solution of LSE problem () can be obtained using direct elimination, the nullspace method and method of weighting. In direct elimination and nullspace methods, the LSE problem is first transformed into unconstrained linear least squares (LLS) problem and then it is solved via normal equations or QR factorization methods. In the method of weighting, a large suitable weighted factor γ is selected such that the weighted residual γ (d -Bx) remains of the same size as that of the residual b -Ax and the constraints is satisfied effectively. Then the solution of the LSE problem is approximated by solving the weighted LLS problem. In [], the author studied the method of weighting for LSE problem and provided a natural criterion of selecting the weighted factor γ such that γ ≥ A  / B  M , where M is the rounding unit. For further details as regards methods of solution for LSE problem (), we refer to [, , , -].
Updating is a process which allow us to approximate the solution of the original problem without solving it afresh. It is useful in applications such as in solving a sequence of modified related problems by adding or removing data from the original problem. Stable and efficient methods of updating are required in various fields of science and engineering such as in optimization and signal processing [ In this article, we consider the LSE problem () in the following form: which is an unconstrained weighted LLS problem where γ ≥ A  / B  M given in [] and approximated its solution by updating Householder QR factorization. It is well known that Householder QR factorization is backward stable [, , , ]. The conditions given in () ensure that problem () is a full rank LLS problem. Hence, there exist a unique solution x(γ ) of problem () which approximated the solution x LSE of the LSE problem () such that lim γ →∞ x(γ ) = x LSE . In our proposed technique, we reduced problem () to a small subproblem using a suitable partitioning strategy by removing blocks of rows and columns. The QR factorization of the subproblem is calculated and its R factor is then updated by appending the removed block of columns and rows, respectively, to approximate the solution of problem (). The presented approach is suitable for solution of dense problems and also applicable for those applications where we are adding new data to the original problem and QR factorization of the problem matrix is already available. We organized this article as follows. Section  contains preliminaries related to our main results. In Section , we present the QR updating procedure and algorithm for solution of LSE problem (). The error analysis of the proposed algorithm is provided in Section . Numerical experiments and comparison of solutions is given in Section , while our conclusion is given in Section .

Preliminaries
This section contains important concepts which will be used in the forthcoming sections.

The method of weighting
This method is based on the observations that while solving LSE problem () we are interested that some equations are to be exactly satisfied. This can be achieved by multiplying large weighted factor γ to those equations. Then we can solve the resulting weighted LLS problem (). The method of weighting is useful as it allows for the use of subroutines for LLS problems to approximate the solution of LSE problem. However, the use of the large weighted factor γ can compromise the conditioning of the constrained matrix. In particular, the method of normal equations when applied to problem () fails for large values of γ in general. For details, see [, -].

QR factorization and householder reflection
Let be the QR factorization of a matrix X ∈ R m×n where Q ∈ R m×m and R ∈ R m×n . This factorization can be obtained using Householder reflections, Givens rotations and the classical/modified Gram-Schmidt orthogonalization. The Householder and Givens QR factorization methods are backward stable. For details as regards QR factorization, we refer to [, ].
Here, we briefly discuss the Householder reflection method due to our requirements. For a non-zero Householder vector v ∈ R n , we define the matrix H ∈ R n×n as which is called the Householder reflection or Householder matrix. Householder matrices are symmetric and orthogonal. For a non-zero vector u, the Householder vector v is simply defined as v = u ± u  e k , such that where α = u  and e k denotes the kth unit vector in R n in the following form: In our present work given in next section, we will use the following algorithm for computation of Householder vector v. This computation is based on Parlett's formula []: where u  >  is the first element of the vector u ∈ R n . Then we compute the Householder vector v as follows (Algorithm ).

Updating QR factorization procedure for solution of LSE problem
Then we can write problem () as Here, we will reduce problem () to a small incomplete subproblem using suitable partition process. In partition process, we remove blocks of rows and columns from the problem matrix E considered in () and from its corresponding right-hand side (RHS) without involving any arithmetics. That is, we consider and removing blocks of rows from both E and f in equation () as we get , .
Algorithm  Calculating the QR factorization and computing matrix U c Input: if k < n then : Hence, we obtain the following problem: where m  = m + pr, n  = n, and f r ∈ R r . Furthermore, by removing block of columns G c = E  (:, j : j + c) from the jth position by considering the partition of E  in the incomplete problem () as we obtain the following reduced subproblem: where E  = [E  (:,  : j -), E  (:, j +  : n  )], n  = n c, m  = m  , and f  = f  . Now, we calculated the QR factorization of the incomplete subproblem () is order to reduce it to the upper triangular matrix R  using the following algorithm (Algorithm ).
Here house denotes the Householder algorithm and the Householder vectors are calculated using Algorithm  and V is a self-explanatory intermediary variable.
Hence, we obtain Here, the QR factorization can also be obtained directly using the MATLAB built-in command qr but it is not preferable due to its excessive storage requirements for orthogonal matrix Q and by not necessarily providing positive sign of diagonal elements in the matrix Algorithm  Updating R  factor after appending a block of columns G c Input: To obtain the solution of problem (), we need to update the upper triangular matrix R  . For this purpose, we append the updated block of columns G c to the R  factor in () at the jth position as follows: Here, if the R  c factor in () is upper trapezoidal or in upper triangular form then no further calculation is required, and we get R  = R c . Otherwise, we will need to reduce equation () to the upper triangular factor R  by introducing the Householder matrices H n  +c , . . . , H j : where R  ∈ R m  ×n  is upper trapezoidal for m  < n  or it is an upper triangular matrix. The procedure for appending the block of columns and updating of the R  factor in algorithmic form is given as follows (Algorithm ).
Here the term triu denotes the upper triangular part of the concerned matrix and V is the intermediary variable. Now, we append a block of rows G r to the R  factor and f r to its corresponding RHS at the jth position in the following manner.
Algorithm  Updating R  factor after appending a block of rows G r Input: if k < n then The equation () is reduced to upper triangular form by constructing the matrices H  , . . . , H n using Algorithm  as given bỹ The procedure is performed in algorithmic form as follows (Algorithm ).
Here qr is the MATLAB command of QR factorization and V is a self-explanatory intermediary variable.
The solution of problem () can then be obtained by applying the MATLAB built-in command backsub for a back-substitution procedure.
The description of QR updating algorithm in compact form is given as follows (Algorithm ).

Algorithm  Updating QR algorithm for solution of LSE problem
Here, Algorithm  for a solution of LSE problem () calls upon the partition process, Algorithms , ,  and MATLAB command backsub for back-substitution procedure, respectively.

Error analysis
In this section, we will study the backward stability of our proposed Algorithm . The mainstay in our presented algorithm for the solution of LSE problem is the updating procedure. Therefore, our main concern is to study the error analysis of the updating steps. For others, such as the effect of using the weighting factor, finding the QR factorization and for the back-substitution procedure, we refer the reader to [, ]. Here, we recall some important results without giving their proofs and refer the reader to [].
We will consider the following standard model of floating point arithmetic in order to analyze the rounding errors effects in the computations: where M is the unit roundoff. The value of M is of order  - in single precision computer arithmetic, while in double precision it is of the order  - . Moreover, for addition and subtraction in the absence of guard integers, we can write model () as follows: Here, we will use the constants γ n andγ cn for convenience as adopted in [] where these are defined by  Here, we represent the Householder transformation as Ivv T , which requires v  = √ . Therefore, by redefining v = √ τ v and τ =  using Lemma ., we havẽ where H = Ivv T and · F denotes the Frobenius norm.

Lemma . ([]) Let H k = Iv k v T k ∈ R m×m be the Householder matrix and we define the sequence of transformations
where X  = X ∈ R m×n . Furthermore, it is assumed that these transformations are performed using computed Householder vectorsṽ k ≈ v k and satisfy Lemma .. Then we havẽ where Q T = H r , . . . , H  , and

Backward error analysis of proposed algorithm
To appreciate the backward stability of our proposed Algorithm , we first need to carry out the error analysis of Algorithms  and . For this purpose, we present the following.

Theorem . The computed factorR in Algorithm
where G r ∈ R r×n is the appended block of rows to the R  factor and Q = H  H  · · · H n .
Proof Let the Householder matrix H j have zeros on the jth column of the matrix r jj G rj . Then using Lemma ., we have Similarly, This implies that where we ignored (γ r+ )  asγ r+ is very small. Continuing in the same fashion till the nth Householder reflection, we havẽ or, by using equation () of Lemma ., we can write which is the required result.
Theorem . The backward error for the computed factor R  in Algorithm  is given by where U c ∈ R (m+p)×c is the appended block of columns to the right-end of R  factor, ẽ c F ≤ γ (m+p)c U c F andQ = H n  + · · · H n  +c .
Proof To prove the required result, we proceed similar to the proof of Theorem ., and obtain This implies that the error in the whole process of appending the block of columns to the R  factor is given by Theorem . Let the LSE problem () satisfy the conditions given in (). Then Algorithm  solves the LSE problem with computedQ = H  · · · H n andR which satisfies Proof To study the backward error analysis of Algorithm , we consider the reduced subproblem matrix E  ∈ R m  ×n  from () and letê  be its backward error in the computed QR factorization. Then from Lemma ., we have As in our proposed Algorithm , we appended a block of columns U c = Q T  G c to the R  factor, then by using Theorem ., we have where e c F ≤γ m  c U c F . Therefore, by simplifying where Q T  = H n  + · · · H n  +c andê  is the error of computing the QR factorization, we obtain Hence, we have which is the total error at this stage of appending the block of columns and its updating. Furthermore, we appended the block of rows G r to the computed factor R  in our algorithm and, by using Theorem ., we obtaiñ where e r F ≤γ (r+)n R  G r F , is the error for appending the block of rows G r . Hence, the total error e for the whole updating procedure in Algorithm  can be written as Now, if the orthogonal factor Q is to be formed explicitly, then the deviation from normality in our updating procedure can be examined. For this, we consider E = I, then from Lemma ., we havẽ where ζ is the error in the computed factorQ  given as where I m  F = √ n  . In a similar manner, the computed factorQ  after appending columns U c is given bỹ Therefore, the total error inQ is given as Similarly, the error in the computed factorQ after appending block of rows is given as So, the total error ζ in Q during the whole updating procedure is given by where m  < m + p, n  < n and c < n. Therefore, the error measure inQ which shows the amount by which it has been deviated from normality is given by From equation (), we havẽ γ (r+)n +γ m  n  +γ m  c ≈γ (m+p+)n =γ (m+p)n , (   ) and using it in equation (), we get the required result (). Also, we have Hence, applying expressions () and () to (), we obtain the required equation () which shows how large is the error in the computed QR factorization.

Numerical experiments
This section is devoted to some numerical experiments which illustrate the applicability and accuracy of Algorithm . The problem matrices A and B and its corresponding right-hand side vectors b and d are generated randomly using the MATLAB built-in commands rand('twister') and rand. These commands generate pseudorandom numbers from a standard uniform distribution in the open interval (, ). The full description of test matrices are given in Table , where we denote the Frobenius norm with · F and the condition number by κ(·). For accuracy of the solution, we consider the actual solution such that x = rand(n, ) and denote the result obtained from Algorithm  by x LSE and that of direct QR Householder factorization with column pivoting by x p . We obtain the relative errors between the solutions given in Table . Moreover, the solution x LSE satisfy the constrained system effectively. The description of the matrix E, the size of the reduced subproblem (SP), value of the weighted factor ω, the relative errors err = xx LSE  / x  and err1 = xx p  / x  are provided in Table . We also carry out the backward error tests of Algorithm  numerically for our considered problems and provide the results in Table , which agrees with our theoretical results.

Conclusion
The solution of linear least squares problems with equality constraints is studied by updated techniques based on QR factorization. We updated only the R factor of the QR  factorization of the small subproblem in order to obtain the solution of our considered problem. Numerical experiments are provided which illustrated the accuracy of the presented algorithm. We also showed that the algorithm is backward stable. The presented approach is suitable for dense problems and also applicable where QR factorization of a problem matrix is available and we are interested in the solution after adding new data to the original problem. In the future, it will be of interest to study the updating techniques for sparse data problems and for those where the linear least squares problem is fixed and the constraint system is changing frequently.