Continuously Updating Matrix Factorization Nikolai Chang

• Journal List
• Springer Open Choice
• PMC5684333

J Inequal Appl. 2017; 2017(1): 281.

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

Salman Zeb

Department of Mathematics, University of Malakand, Dir (Lower), Chakdara, Khyber Pakhtunkhwa Pakistan

Muhammad Yousaf

Department of Mathematics, University of Malakand, Dir (Lower), Chakdara, Khyber Pakhtunkhwa Pakistan

Received 2017 Jul 5; Accepted 2017 Oct 16.

Abstract

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.

Keywords: QR factorization, orthogonal transformation, updating, least squares problems, equality constraints

Introduction

We consider the linear least squares problem with equality constraints (LSE)

min x ∥ A x − b ∥ 2 , subject to B x = d ,

1

where A ∈R ^m×n , b ∈R ^m ,B ∈R ^p×n , d ∈R ^p , x ∈R ⁿ 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 in the analysis of large scale structure [1–8]. The assumptions

rank(B) =p and null(A) ∩ null(B) = {0}

2

ensure the existence and uniqueness of solution for problem (1).

The solution of LSE problem (1) 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 −B x) remains of the same size as that of the residual b −A x and the constraints is satisfied effectively. Then the solution of the LSE problem is approximated by solving the weighted LLS problem. In [9], 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 (1), we refer to [2, 3, 7, 10–13].

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 [14], statistics [15], network and structural analysis [16, 17] and discretization of differential equations [18]. Various updating techniques based on matrix factorization for different kinds of problems exist in the literature [3, 11, 13, 19–22]. Hammarling and Lucas [23] discussed updating the QR factorization with applications to LLS problem and presented updating algorithms which exploited LEVEL 3 BLAS. Yousaf [24] studied repeated updating based on QR factorization as a solution tool for LLS problem. Parallel implementation on GPUs of the updating QR factorization algorithms presented in [23] is performed by Andrew and Dingle [25]. Zhdanov and Gogoleva [26] studied augmented regularized normal equations for solving LSE problems. Zeb and Yousaf [27] presented an updating algorithm by repeatedly updating both factors of the QR factorization for the solutions of LSE problems.

In this article, we consider the LSE problem (1) in the following form:

which is an unconstrained weighted LLS problem where γ ≥ ∥A∥₂/∥B∥₂ ϵ _M given in [9] and approximated its solution by updating Householder QR factorization. It is well known that Householder QR factorization is backward stable [11, 12, 28, 29]. The conditions given in (2) ensure that problem (3) is a full rank LLS problem. Hence, there exist a unique solution x(γ) of problem (3) which approximated the solution x _LSE of the LSE problem (1) such that lim _γ→∞ x(γ) =x _LSE . In our proposed technique, we reduced problem (3) 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 (1). 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 2 contains preliminaries related to our main results. In Section 3, we present the QR updating procedure and algorithm for solution of LSE problem (1). The error analysis of the proposed algorithm is provided in Section 4. Numerical experiments and comparison of solutions is given in Section 5, while our conclusion is given in Section 6.

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 (1) 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 (3). 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 (3) fails for large values of γ in general. For details, see [3, 11–13].

QR factorization and householder reflection

Let

be the QR factorization of a matrix X ∈ ℛ ^m×n where Q ∈ ℛ ^m×m and 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 [11, 12].

Here, we briefly discuss the Householder reflection method due to our requirements. For a non-zero Householder vector v ∈ ℛ ⁿ , we define the matrix H ∈ ℛ ^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

such that

where α = ∥u∥₂ and e _k denotes the kth unit vector in ℛ ⁿ in the following form:

e k ( i ) = { 1 if i = k , 0 otherwise .

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 [30]:

v = u 1 − ∥ u ∥ 2 = u 1 2 − ∥ u ∥ 2 2 x 1 + ∥ u ∥ 2 = − ∥ u ∥ 2 2 u 1 + ∥ u ∥ 2 ,

where u ₁ > 0 is the first element of the vector u ∈ ℛ ⁿ . Then we compute the Householder vector v as follows (Algorithm 1).

Algorithm 1

Computation of parameter τ and Householder vector v [12]

Algorithm 2

Calculating the QR factorization and computing matrix U _c

Algorithm 3

Updating R ₂ factor after appending a block of columns G _c

Algorithm 4

Updating R ₁ factor after appending a block of rows G _r

Algorithm 5

Updating QR algorithm for solution of LSE problem

Updating QR factorization procedure for solution of LSE problem

Let

E = ( γ B A ) ∈ R ( m + p ) × n and f = ( γ d b ) ∈ R m + p .

Then we can write problem (3) as

Here, we will reduce problem (7) 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 (7) and from its corresponding right-hand side (RHS) without involving any arithmetics. That is, we consider

E = ( E ( 1 : j − 1 , 1 : n ) E ( j : j + r , 1 : n ) E ( j + r + 1 : m + p , 1 : n ) ) and f = ( f ( 1 : j − 1 ) f ( j : j + r ) f ( j + r + 1 : m + p ) ) ,

8

and removing blocks of rows from both E and f in equation (8) as

G _r =E(j:j +r, 1:n) ∈ ℛ ^r×n and f _r =f(j:j +r),

9

we get

E 1 = ( E ( 1 : j − 1 , 1 : n ) E ( j + r + 1 : m + p , 1 : n ) ) , f 1 = ( f ( 1 : j − 1 ) f ( j + r + 1 : m + p ) ) .

10

Hence, we obtain the following problem:

min x 1 ∥ E 1 x 1 − f 1 ∥ 2 , E 1 ∈ R m 1 × n 1 , f 1 ∈ R m 1 , x 1 ∈ R n 1 ,

11

where m ₁ =m +p −r , n ₁ =n , and f _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 (11) as

E 2 = ( E 1 ( : , 1 : j − 1 ) E 1 ( : , j : j + c ) E 1 ( : , j + c + 1 : n ) ) ,

12

we obtain the following reduced subproblem:

min x 2 ∥ E 2 x 2 − f 2 ∥ 2 , E 2 ∈ R m 2 × n 2 , f 2 ∈ R m 2 , x 2 ∈ R n 2 ,

13

where E ₂ = [E ₁(:, 1:j − 1),E ₁(:,j + 1:n ₁)], n ₂ =n ₁ −c , m ₂ =m ₁ , and f ₁ =f ₂ .

Now, we calculated the QR factorization of the incomplete subproblem (13) is order to reduce it to the upper triangular matrix R ₂ using the following algorithm (Algorithm 2).

Here house denotes the Householder algorithm and the Householder vectors are calculated using Algorithm 1 and V is a self-explanatory intermediary variable.

Hence, we obtain

R ₂ =H _{n 2}  ⋯H ₁ E ₂, g ₂ =H _{n 2}  ⋯H ₁ f ₂

14

and

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 R ₂ [12].

To obtain the solution of problem (7), 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 (14) at the jth position as follows:

R 2 c = ( R 2 ( : , 1 : j − 1 ) G c ( : , j : j + c ) R 2 ( : , j + c + 1 : n ) ) .

15

Here, if the R ₂ c factor in (15) is upper trapezoidal or in upper triangular form then no further calculation is required, and we get R ₁ =R _2c . Otherwise, we will need to reduce equation (15) to the upper triangular factor R ₁ by introducing the Householder matrices H _{n 2+c} , …,H _j :

R ₁ =H _{n 2+c}  ⋯H _j R _2c and g ₁ =H _{n 2+c}  ⋯H _j g ₂,

16

where 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 3).

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.

R r = ( R 1 ( 1 : j − 1 , : ) G r ( 1 : r , : ) R 1 ( j : m 1 , : ) ) , g r = ( g 1 ( 1 : j − 1 ) f r ( 1 : r ) g 1 ( j : m 1 ) ) .

We use the permutation matrix P in order to bring the block of rows G _r and f _r at the bottom position if required. Then

The equation (17) is reduced to upper triangular form by constructing the matrices H ₁, …,H _n using Algorithm 1 as given by

and

The procedure is performed in algorithmic form as follows (Algorithm 4).

Here qr is the MATLAB command of QR factorization and V is a self-explanatory intermediary variable.

The solution of problem (7) 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 5).

Here, Algorithm 5 for a solution of LSE problem (1) calls upon the partition process, Algorithms 2, 3, 4 and MATLAB command backsub for back-substitution procedure, respectively.

Error analysis

In this section, we will study the backward stability of our proposed Algorithm 5. 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 [28, 31]. Here, we recall some important results without giving their proofs and refer the reader to [28].

We will consider the following standard model of floating point arithmetic in order to analyze the rounding errors effects in the computations:

f l(x opy) = (x opy)(1 +η),  |η| ≤ϵ _M ,  op =  + , −,  ∗ , /,

18

where ϵ _M is the unit roundoff. The value of ϵ _M is of order 10⁻⁸ in single precision computer arithmetic, while in double precision it is of the order 10⁻¹⁶. Moreover, for addition and subtraction in the absence of guard integers, we can write model (18) as follows:

f l(x ±y) =x(1 +η ₁)±y(1 +η ₂),  |η _i | ≤ϵ _M , i = 1, 2.

Lemma 4.1

([28])

If |η _i | ≤ϵ _M ,δ _i = ±1 fori = 1, 2, …,n andn ϵ _M  < 1, then

where

Here, we will use the constants γ _n and γ ˜ c n for convenience as adopted in [28] where these are defined by

γ n = n ϵ M 1 − n ϵ M , assuming n ϵ M < 1 ,

19

and

γ ˜ c n = n c ϵ M 1 − n c ϵ M , assuming n c ϵ M < 1 ,

20

where c is a small positive integer.

Lemma 4.2

([28])

Let x,y ∈ ℛ ⁿ and consider the inner product x ^T y . Then

f l(x ^T y) =x ^T (y + △y),  |△y| ≤γ _n |y|.

Lemma 4.3

([28])

Let ϕ _k and γ _k be defined for any positive integer k. Then for positive integers j and k the following relations hold:

(1 +ϕ _k )(1 +ϕ _j ) = 1 +ϕ _k+j ,

21

( 1 + ϕ k ) ( 1 + ϕ j ) = { 1 + ϕ k + j , j ≤ k , 1 + ϕ k + 2 j , j > k ,

22

Lemma 4.4

([28])

Considering the construction of τ ∈ ℛ and v ∈ ℛ ⁿ given in Section 2, then the computed $\tilde{\tau }$ and ṽ in floating point arithmetic satisfy

and

τ ˜ = τ ( 1 + ϕ ˜ n ) and v 1 ˜ = v 1 ( 1 + ϕ ˜ n ) ,

where ϕ ˜ n ≤ γ ˜ n .

Here, we represent the Householder transformation as I −v v ^T , which requires ∥ v ∥ 2 = 2 . Therefore, by redefining v = τ v and τ = 1 using Lemma 4.4, we have

v ˜ = v + △ v , | △ v | ≤ γ ˜ m v for v ∈ R m , ∥ v ∥ 2 = 2 .

26

Lemma 4.5

([29])

Considering the computation of y = H ˜ b = ( I − v ˜ v ˜ T ) b = b − v ˜ ( v ˜ T b ) for b ∈ ℛ ^m and v ˜ ∈ R m satisfies (26). Then the computed ỹ satisfies

where H =I −v v ^T and ∥⋅∥ _F denotes the Frobenius norm.

Lemma 4.6

([28])

Let H k = I − v k v k T ∈ R m × m be the Householder matrix and we define the sequence of transformations

where X ₁ =X ∈ ℛ ^m×n . Furthermore, it is assumed that these transformations are performed using computed Householder vectors v ˜ k ≈ v k and satisfy Lemma 4.4. Then we have

where Q ^T =H _r , …,H ₁ , and

Theorem 4.7

([28])

Let R ˜ ∈ R m × n be the computed upper trapezoidal QR factor of X ∈ ℛ ^m×n (m ≥n) obtained via Householder QR algorithm by Lemma 4.6. Then there exists an orthogonal matrix Q ∈ ℛ ^m×m such that

where

∥ △ x j ∥ 2 ≤ γ ˜ m n ∥ x j ∥ 2 , j = 1 : n .

28

Lemma 4.8

([28])

Let X ∈ ℛ ^m×n and H =I −τ v v ^T  ∈R ^m×m be the Householder matrix. Also, assuming that the computation of HX is performed using computed $\tilde{\tau }$ and ṽ such that it satisfies Lemma 4.4. Then, from Theorem 4.7, we have

f l(H X) =H(X + △X),  ∥△X∥ _F  ≤γ _{c m} ∥X∥ _F .

29

Backward error analysis of proposed algorithm

To appreciate the backward stability of our proposed Algorithm 5, we first need to carry out the error analysis of Algorithms 3 and 4. For this purpose, we present the following.

Theorem 4.9

The computed factor $\tilde{R}$ in Algorithm 4 satisfies

R ˜ = Q T [ ( R 1 G r ) + e r ] and ∥ e r ∥ F ≤ γ ˜ n ( r + 1 ) ∥ ( R 1 G r ) ∥ F ,

where G _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 ˜ j j G r j ) . Then using Lemma 4.6, we have

H 1 [ ( R 1 G r ) + e 1 ] = ( R ˜ 1 r G ˜ 1 r ) ,

30

where

Similarly,

H 2 [ ( R ˜ 1 r G ˜ 1 r ) + e 2 ] = ( R ˜ 2 r G ˜ 2 r ) = H 2 H 1 [ ( R 1 G r ) + e r ( 2 ) ] ,

and

∥ e 2 ∥ F ≤ γ ˜ r + 1 ∥ ( R ˜ 1 r G ˜ 1 r ) ∥ F ,

where

∥ e r ( 2 ) ∥ F = ∥ e 1 + e 2 ∥ F , ∥ e r ( 2 ) ∥ F ≤ ∥ e 1 ∥ F + ∥ e 2 ∥ F , ≤ γ ˜ r + 1 ∥ ( R 1 G r ) ∥ F + γ ˜ r + 1 ∥ ( R ˆ 1 r G ˆ 1 r ) ∥ F .

From equation (30), we have

∥ ( R ˆ 1 r G ˆ 1 r ) ∥ F ≤ ∥ ( R 1 G r ) ∥ F + e 1 ≤ ∥ ( R 1 G r ) ∥ F + γ ˜ r + 1 ∥ ( R 1 G r ) ∥ F .

This implies that

∥ e r ( 2 ) ∥ F ≤ γ ˜ r + 1 ∥ ( R 1 G r ) ∥ F + γ ˜ r + 1 ( 1 + γ ˜ r + 1 ) ∥ ( R 1 G r ) ∥ F ≤ ( γ ˜ r + 1 + γ ˜ r + 1 ( 1 + γ ˜ r + 1 ) ) ∥ ( R 1 G r ) ∥ F ≤ 2 γ ˜ r + 1 ∥ ( R 1 G r ) ∥ F ,

where we ignored ( γ ˜ r + 1 ) 2 as γ ˜ r + 1 is very small.

Continuing in the same fashion till the nth Householder reflection, we have

where

∥ e r ∥ F = ∥ e r ( n ) ∥ F ≤ n γ ˜ r + 1 ∥ ( R 1 G r ) ∥ F ,

or, by using equation (23) of Lemma 4.3, we can write

which is the required result. □

Theorem 4.10

The backward error for the computed factor R ₁ in Algorithm 3 is given by

where U _c  ∈ ℛ^(m+p)×c is the appended block of columns to the right-end of R ₂ factor, ∥ e ˜ c ∥ F ≤ γ ˜ ( m + p ) c ∥ U c ∥ F and Q ˜ = H n 2 + 1 ⋯ H n 2 + c .

Proof

To prove the required result, we proceed similar to the proof of Theorem 4.9, and obtain

∥ e c 2 ∥ F ≤ 2 γ ˜ c 1 ∥ [ ( R 2 , U c ) + e c ˜ ] ∥ F .

This implies that the error in the whole process of appending the block of columns to the R ₂ factor is given by

R 1 = H n 2 + c ⋯ H n 2 + 1 [ ( R 2 , U c ) + e c ˜ ] ,

where

 □

Theorem 4.11

Let the LSE problem (1) satisfy the conditions given in (2). Then Algorithm5 solves the LSE problem with computed Q ˜ = H 1 ⋯ H n and $\tilde{R}$ which satisfies

and

∥ E − Q ˜ R ˜ ∥ F ≤ n γ ˜ ( m + p ) n ∥ E ∥ F .

32

Proof

To study the backward error analysis of Algorithm 5, we consider the reduced subproblem matrix E ₂ ∈ ℛ ^{m ₂×n ₂} from (13) and let e ˆ 2 be its backward error in the computed QR factorization. Then from Lemma 4.8, we have

Q 2 T ( E 2 + e ˆ 2 ) = R 2 and ∥ e ˆ 2 ∥ F ≤ γ ˜ m 2 n 2 ∥ E 2 ∥ F .

As in our proposed Algorithm 5, we appended a block of columns U c = Q 2 T G c to the R ₂ factor, then by using Theorem 4.10, we have

R ₁ =H _{n 2+1} ⋯H _{n 2+c} [(R ₂,U _c ) +e _c ],

where ∥ e c ∥ F ≤ γ ˜ m 2 c ∥ U c ∥ F .

Therefore, by simplifying

∥ e ˆ 1 ∥ F = ∥ Q 1 T Q 2 T e ˆ 2 + Q 1 T e c ∥ F ,

where Q 1 T = H n 2 + 1 ⋯ H n 2 + c and e ˆ 2 is the error of computing the QR factorization, we obtain

∥ e ˆ 1 ∥ F = ∥ Q 1 T Q 2 T e ˆ 2 + Q 1 T e c ∥ F ≤ ∥ Q 1 T ∥ F ∥ Q 2 T ∥ F ∥ e ˆ 2 ∥ F + ∥ Q 1 T ∥ F ∥ e c ∥ F ≤ ∥ e ˆ 2 ∥ F + ∥ e c ∥ F ≤ γ ˜ m 2 n 2 ∥ E 2 ∥ F + γ ˜ m 2 c ∥ U c ∥ F ≤ [ γ ˜ m 2 n 2 + γ ˜ m 2 c ] max ( ∥ E 2 ∥ F , ∥ U c ∥ F ) .

Hence, we have

∥ e ˆ 1 ∥ F ≤ [ γ ˜ m 2 n 2 + γ ˜ m 2 c ] max ( ∥ E 2 ∥ F , ∥ U c ∥ F ) ,

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 4.9, we obtain

where

is the error for appending the block of rows G _r .

Hence, the total error e for the whole updating procedure in Algorithm 5 can be written as

∥ e ∥ F ≤ ∥ e r ∥ F + ∥ e ˆ 1 ∥ F ≤ γ ˜ ( r + 1 ) n ∥ ( R 1 G r ) ∥ F + ( γ ˜ m 2 n 2 + γ ˜ m 2 c ) max ( ∥ E 2 ∥ F , ∥ U c ∥ F ) ≤ ( γ ˜ ( r + 1 ) n + γ ˜ m 2 n 2 + γ ˜ m 2 c ) max ( ∥ E 2 ∥ F , ∥ U c ∥ F , ∥ ( R 1 G r ) ∥ F ) .

33

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 4.8, we have

where ζ is the error in the computed factor Q ˜ 2 given as

∥ ζ ∥ F ≤ γ ˜ m 2 n 2 ∥ I m 2 ∥ F = n 2 γ ˜ m 2 n 2 ,

where ∥ I m 2 ∥ F = n 2 . In a similar manner, the computed factor Q ˜ 1 after appending columns U _c is given by

where

Therefore, the total error in $\tilde{Q}$ is given as

where

∥ ζ 1 ∥ F ≤ ∥ ζ 2 ∥ F + ∥ ζ c ∥ F ≤ n 2 γ ˜ m 2 n 2 + c γ ˜ m 2 c .

Similarly, the error in the computed factor $\tilde{Q}$ after appending block of rows is given as

where

So, the total error ζ in Q during the whole updating procedure is given by

∥ ζ ∥ F = ∥ ζ 1 ∥ F + ∥ ζ r ∥ F ≤ n γ ˜ ( r + 1 ) n + n 2 γ ˜ m 2 n 2 + c γ ˜ m 2 c ,

34

where m ₂ <m +p , n ₂ <n and c <n .

Therefore, the error measure in $\tilde{Q}$ which shows the amount by which it has been deviated from normality is given by

∥ I − Q ˜ T Q ˜ ∥ F ≤ n ( γ ˜ ( r + 1 ) n + γ ˜ m 2 n 2 + γ ˜ m 2 c ) .

35

From equation (20), we have

γ ˜ ( r + 1 ) n + γ ˜ m 2 n 2 + γ ˜ m 2 c ≈ γ ˜ ( m + p + 1 ) n = γ ˜ ( m + p ) n ,

36

and using it in equation (35), we get the required result (31).

Also, we have

∥ E − Q ˜ R ˜ ∥ F = ∥ ( E − Q R ˜ ) + ( ( Q − Q ˜ ) R ˜ ) ∥ F ≤ n ( γ ˜ ( r + 1 ) n + γ ˜ m 2 n 2 + γ ˜ m 2 c ) max ( ∥ E 2 ∥ F , ∥ E c ∥ F , ∥ ( R 1 G r ) ∥ F ) .

37

As ∥X∥ _F = ∥QR∥ _F = ∥R∥ _F , therefore, we can write

max ( ∥ E 2 ∥ F , ∥ U c ∥ F , ∥ ( R 1 G r ) ∥ F ) = ∥ ( R 1 G r ) ∥ F = ∥ E ∥ F .

38

Hence, applying expressions (36) and (38) to (37), we obtain the required equation (32) 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 5. 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 (0, 1). The full description of test matrices are given in Table1, 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, 1) and denote the result obtained from Algorithm 5 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 Table2. 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 = ∥x−x _LSE∥₂/∥x∥₂ and err1 = ∥x−x _p ∥₂/∥x∥₂ are provided in Table2. We also carry out the backward error tests of Algorithm 5 numerically for our considered problems and provide the results in Table3, which agrees with our theoretical results.

Table 1

Description of test problems

Problem	Size of (A)	κ ( A )	∥ A ∥ F	Size of (B)	κ ( B )	∥ B ∥ F
1.	10×8	1.3667e+02	2.0006e+02	6×8	7.4200e+01	1.0216e+02
2.	100×90	2.9303e+03	2.1395e+03	90×90	3.3687e+03	1.3735e+03
3.	800×700	6.2106e+03	1.6872e+04	600×700	1.6164e+03	9.9000e+03
4.	1,000×500	1.1602e+03	1.5943e+04	500×500	1.2883e+05	7.6360e+03
5.	2,000×1,000	1.6727e+03	3.1884e+04	1,000×1,000	1.7430e+06	1.5272e+04

Table 2

Results comparison

Problem	Size of (E)	γ	Size of SP	err1	err
1.	16×8	8.9564e+15	3×3	1.3222e−14	1.4585e−15
2.	190×90	7.1258e+15	3×3	1.2628e−13	5.5294e−14
3.	1,400×700	7.7993e+15	3×3	1.2821e−12	4.2522e−13
4.	1,500×500	9.5551e+15	3×3	1.9377e−12	1.3559e−12
5.	3,000×1,000	9.5549e+15	3×3	1.0828e−10	8.5181e−12

Table 3

Backward error analysis results

Problem	∥ E − Q ˜ R ˜ ∥ F ∥ E ∥ F	∥ I − Q ˜ T Q ˜ ∥ F
1.	4.4202e−16	1.3174e−15
2.	4.7858e−16	9.0854e−15
3.	1.0450e−15	4.9428e−14
4.	9.0230e−16	3.8711e−14
5.	9.9304e−16	6.4026e−14

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.

Authors' contributions

All authors contributed equally to this work. All authors read and approved the final manuscript.

Notes

Competing interests

The authors declare that they have no competing interests.

Footnotes

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Contributor Information

Salman Zeb, Email: moc.liamg@namlasbez.

Muhammad Yousaf, Email: moc.liamg@moufasuoym.

References

1. Auton, JR, Van Blaricum, ML: Investigation of procedures for automatic resonance extraction from noisy transient electromagnetics data. Final Report, Contract N00014-80-C-029, Office of Naval Research Arlington, VA 22217 (1981)

2. Barlow JL, Handy SL. The direct solution of weighted and equality constrained least-squares problems. SIAM J. Sci. Stat. Comput. 1988;9:704–716. doi: 10.1137/0909046. [CrossRef] [Google Scholar]

3. Lawson CL, Hanson RJ. Solving Least Squares Problems. Philadelphia: SIAM; 1995. [Google Scholar]

4. Eldén L. Perturbation theory for the least squares problem with linear equality constraints. SIAM J. Numer. Anal. 1980;17:338–350. doi: 10.1137/0717028. [CrossRef] [Google Scholar]

5. Leringe, O, Wedin, PÅ: A comparison between different methods to compute a vector x which minimizes ∥A x−b∥₂∥A x−b∥₂ when G x =h G x =h. Tech. Report, Department of computer science, Lund University (1970)

6. Van Loan CF. A generalized SVD analysis of some weighting methods for equality constrained least squares. In: Kågström B, Ruhe A, editors. Matrix Pencils. Heidelberg: Springer; 1983. pp. 245–262. [Google Scholar]

7. Van Loan CF. On the method of weighting for equality-constrained least-squares problems. SIAM J. Numer. Anal. 1985;22:851–864. doi: 10.1137/0722051. [CrossRef] [Google Scholar]

8. Wei M. Algebraic properties of the rank-deficient equality-constrained and weighted least squares problems. Linear Algebra Appl. 1992;161:27–43. doi: 10.1016/0024-3795(92)90003-S. [CrossRef] [Google Scholar]

9. Stewart GW. On the weighting method for least squares problems with linear equality constraints. BIT Numer. Math. 1997;37:961–967. doi: 10.1007/BF02510363. [CrossRef] [Google Scholar]

10. Anderson E, Bai Z, Dongarra J. Generalized QR factorization and its applications. Linear Algebra Appl. 1992;162:243–271. doi: 10.1016/0024-3795(92)90379-O. [CrossRef] [Google Scholar]

11. Björck Å. Numerical Methods for Least Squares Problems. Philadelphia: SIAM; 1996. [Google Scholar]

12. Golub GH, Van Loan CF. Matrix Computations. Baltimore: Johns Hopkins University Press; 1996. [Google Scholar]

13. Stewart GW. Matrix Algorithms. Philadelphia: SIAM; 1998. [Google Scholar]

14. Alexander ST, Pan CT, Plemmons RJ. Analysis of a recursive least squares hyperbolic rotation algorithm for signal processing. Linear Algebra Appl. 1988;98:3–40. doi: 10.1016/0024-3795(88)90158-9. [CrossRef] [Google Scholar]

15. Chambers JM. Regression updating. J. Am. Stat. Assoc. 1971;66:744–748. doi: 10.1080/01621459.1971.10482338. [CrossRef] [Google Scholar]

16. Haley SB, Current KW. Response change in linearized circuits and systems: computational algorithms and applications. Proc. IEEE. 1985;73:5–24. doi: 10.1109/PROC.1985.13106. [CrossRef] [Google Scholar]

17. Haley SB. Solution of modified matrix equations. SIAM J. Numer. Anal. 1987;24:946–951. doi: 10.1137/0724061. [CrossRef] [Google Scholar]

18. Lai SH, Vemuri BC. Sherman-Morrison-Woodbury-formula-based algorithms for the surface smoothing problem. Linear Algebra Appl. 1997;265:203–229. doi: 10.1016/S0024-3795(97)80366-7. [CrossRef] [Google Scholar]

19. Björck Å, Eldén L, Park H. Accurate downdating of least squares solutions. SIAM J. Matrix Anal. Appl. 1994;15:549–568. doi: 10.1137/S089547989222895X. [CrossRef] [Google Scholar]

20. Daniel JW, Gragg WB, Kaufman L, Stewart GW. Reorthogonalization and stable algorithms for updating the Gram-Schmidt QR factorization. Math. Comput. 1976;30:772–795. [Google Scholar]

21. Gill PE, Golub GH, Murray W, Saunders M. Methods for modifying matrix factorizations. Math. Comput. 1974;28:505–535. doi: 10.1090/S0025-5718-1974-0343558-6. [CrossRef] [Google Scholar]

22. Reichel L, Gragg WB. Algorithm 686: FORTRAN subroutines for updating the QR decomposition. ACM Trans. Math. Softw. 1990;16:369–377. doi: 10.1145/98267.98291. [CrossRef] [Google Scholar]

23. Hammarling, S, Lucas, C: Updating the QR factorization and the least squares problem. Tech. Report, The University of Manchester (2008). http://www.manchester.ac.uk/mims/eprints

24. Yousaf, M: Repeated updating as a solution tool for linear least squares problems. Dissertation, University of Essex (2010)

25. Andrew R, Dingle N. Implementing QR factorization updating algorithms on GPUs. Parallel Comput. 2014;40:161–172. doi: 10.1016/j.parco.2014.03.003. [CrossRef] [Google Scholar]

26. Zhdanov AI, Gogoleva SY. Solving least squares problems with equality constraints based on augmented regularized normal equations. Appl. Math. E-Notes. 2015;15:218–224. [Google Scholar]

27. Zeb S, Yousaf M. Repeated QR updating algorithm for solution of equality constrained linear least squares problems. Punjab Univ. J. Math. 2017;49:51–61. [Google Scholar]

28. Higham NJ. Accuracy and Stability of Numerical Algorithms. Philadelphia: SIAM; 2002. [Google Scholar]

29. Wilkinson JH. Error analysis of transformations based on the use of matrices of the form I − 2w w ^H I − 2w w ^H . In: Rall LB, editor. Error in Digital Computation. New York: Wiley; 1965. pp. 77–101. [Google Scholar]

30. Parlett BN. Analysis of algorithms for reflections in bisectors. SIAM Rev. 1971;13:197–208. doi: 10.1137/1013037. [CrossRef] [Google Scholar]

31. Cox, AJ: Stability of algorithms for solving weighted and constrained least squares problems. Dissertation, University of Manchester (1997)

vernonwoemorre.blogspot.com

Source: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5684333/

Share this post