# Norm-Constrained Least-Squares Solutions to the Matrix EquationAXB=C

Authors: An-bao Xu; Zhenyun Peng
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

An iterative method to compute the least-squares solutions of the matrixAXB=C over the norm inequality constraint is proposed. For this method, without the error of calculation, a desired solution can be obtained with finitely iterative step. Numerical experiments are performed to illustrate the efficiency and real application of the algorithm.


---

## Body

## 1. Introduction

Throughout this paper,Rm×n denotes the set of all m×n real matrices. I represents the identity matrix of size implied by context. AT,∥A∥ denote, respectively, the transpose and the Frobenius norm of the matrix A. For the matrices A=(aij)∈Rm×n,B=(bij)∈Rp×q, A⊗B represents the Kronecker production of the matrices A and B defined as A⊗B=(aijB)∈Rmp×nq. The inner product in the matrix set space Rm×n is defined as 〈A,B〉=trace(BTA) for all the matrices A,B∈Rm×n. Obviously, Rm×n is a Hilbert inner product space and the norm of a matrix generated by this inner product space is the Frobenius norm.SolutionsX to the well-known linear matrix equation AXB=C with special structures have been widely studied. See, for example, [1–5] for symmetric solutions, skew-symmetric solutions, centro-symmetric solutions, symmetric R-symmetric solutions, or (R,S)-symmetric solutions. To the best of our knowledge, the solutions to the matrix equation AXB=C subject to the norm inequality constraint, however, have not been studied directly in the literature. In this paper we consider the solutions to the following least-squares problem over the norm inequality constraint:
(1)minX∈Rm×n12∥AXB-C∥2subjectto∥X∥≤Δ,
where A∈Rp×m, B∈Rn×q, C∈Rp×q and Δ is a nonnegative real number.Problem (1) can be regarded as a natural generalization of the unconstrained least-squares problem of the matrix equation AXB=C. In fact, when we let Δ be big enough, the problem will turn out to be the unconstrained least-squares problem of the matrix equation AXB=C. According to [6], moreover, the problem (1) is equivalent to the classical Tikhonov regulation approach of the matrix equation AXB=C(2)minX∈Rm×n12∥AXB-C∥2+ζ∥X∥2,
where ζ>0 is the regularization parameter. While Tikhonov regularization involves the computation of a parameter that does not necessarily have a physical meaning in most problems, the problem (1) has the advantage that, in some applications, the physical properties of the problem either determine or make it easy to estimate an optimal value for the norm constraint Δ. This is the case, for example, in image restoration where Δ represents the energy of the target image [7].In this paper, an iterative method is proposed to compute the solutions of the problem (1). We will use the generalized Lanczos trust region algorithm (GLTR) [8], which is based on Steihaug-Toint algorithm [9, 10], as the frame method for deriving this iterative method. The basic idea is as follows. First, by using the Kronecker production of matrices, we transform the least-squares problem (1) into the trust-region subproblem in vector form which can be solved by the GLTR algorithm. Then, we transform the vector iterative method into matrix form. In the end, numerical experiments are given to illustrate the efficiency and real application of the proposed iteration algorithm.


## 2. Iteration Methods to Solve Problem (1)

In this section we first give the necessary and sufficient conditions for the problem (1) to have a solution. Then we propose an iteration method to compute the solution to the problem. And some properties of this algorithm are also given.Obviously, problem (1) is equivalent to the following problem
(3)minX∈Rm×n12〈AXB,AXB〉-〈AXB,C〉subjectto∥X∥≤Δ.
This equivalent form of the problem (1) makes us more convenient to prove the following theorem.Theorem 1.
MatrixX* is a solution of the problem (3) if and only if there is a scalar λ*≥0 such that the following conditions are satisfied:
(4)ATAX*BBT+λ*X*=ATCBT,λ*(∥X*∥-Δ)=0.Proof.
Assume that there is a scalarλ*≥0 such that the conditions (4) are satisfied. Let
(5)φ(X)=12〈AXB,AXB〉-〈AXB,C〉,φ∧(X)=12〈AXB,AXB〉+12λ*〈X,X〉-〈AXB,C〉=φ(X)+12λ*〈X,X〉.
For any matrix W∈Rm×n, we have
(6)φ∧(X*+W)=12〈A(X*+W)B,A(X*+W)B〉+12λ*〈(X*+W),(X*+W)〉-〈A(X*+W)B,C〉=(12〈AX*B,AX*B〉+12λ*〈X*,X*〉-〈AX*B,C〉)+〈AWB,AX*B〉+λ*〈W,X*〉-〈AWB,C〉+12〈AWB,AWB〉+12λ*〈W,W〉=φ∧(X*)+〈W,(ATAX*BBT+λ*X*-ATCBT)〉+12〈AWB,AWB〉+12λ*〈W,W〉=φ∧(X*)+12〈AWB,AWB〉+12λ*〈W,W〉≥φ∧(X*).
This implies that X* is a global minimizer of the function φ∧(X). Since φ∧(X)≥φ∧(X*) for all X∈Rm×n, we have
(7)φ(X)≥φ(X*)+12λ*(〈X*,X*〉-〈X,X〉).
The equality λ*(∥X*∥-Δ)=0 implies that λ*(〈X*,X*〉-Δ2)=0. Consequently, the following inequality always holds:
(8)φ(X)≥φ(X*)+12λ*(Δ2-〈X,X〉).
Hence, from λ*≥0, we have φ(X)≥φ(X*) for all X∈Rm×n with ∥X∥≤Δ. And so X* is a global minimizer of (3).
Conversely, assuming thatX* is a global solution of the problem (3), we show that there is a nonnegative λ* such that satisfies conditions (4). For this purpose we consider two cases: ∥X*∥<Δ and ∥X*∥=Δ.
In case∥X*∥<Δ, X* is certainly an unconstrained minimizer of φ(X). So X* satisfies the stationary point condition ∇φ(X*)=0; that is, ATAX*BBT-ATCBT=0. This implies that the properties (4) hold for λ*=0. In the case ∥X*∥=Δ, the second equality is immediately satisfied, and X* also solves the constrained problem
(9)minX∈Rm×nφ(X)subjectto∥X∥=Δ.
By applying optimality conditions for constrained optimization to this problem, we know that there exists a scalar λ* such that the Lagrangian function defined by
(10)ζ(X,λ)=φ(X)+12λ(〈X,X〉-Δ2)
has a stationary point at X*. By setting ∇Xζ(X*,λ*) to zero, we obtain
(11)ATAX*BBT-ATCBT+λ*X*=0.
Now the proof is concluded by showing that λ*≥0. Since the equality (11) holds, then X* minimizes φ∧(X), and so we have
(12)φ(X)≥φ(X*)+12λ*(〈X*,X*〉-〈X,X〉)
for all X∈Rm×n. Suppose that there are only negative values of λ* that satisfy (11). Then we have from (12) that
(13)φ(X)≥φ(X*)whenever∥X∥≥∥X*∥=Δ.
Since we already know that X* minimizes φ(X) for ∥X∥≤Δ, it follows that X* is in fact a global, unconstrained minimize of φ(X). Therefore conditions (11) hold when λ*=0, which contradicts our assumption that only negative values of λ* can satisfy condition (11). The proof is completed.We give an iteration method to solve problem (1) as in Algorithm 2.Algorithm 2.
(i) Given matricesX0=0, Q-1=0 and a small tolerance ε>0. ComputeR0=-ATCBT,t0=-ATCBT,γ0=∥R0∥,P0=-R0,T-1=[].Setk←0.
(ii) ComputingQk=tk/γk, δk=∥AQkB∥2, tk+1=ATAQkBBT-δkQk-γkQk-1, γk+1=∥tk+1∥Tk=[Tk-1ΓkΓkTδk], where Γk=(0,…,0,γk)T∈Rk.
(iii) IfAPkB≠0, compute αk=∥Rk∥2/∥APkB∥2.If∥Xk+αkPk∥≤Δ, computing Rk+1=Rk+αkATAPkBBT, βk=∥Rk+1∥2/∥Rk∥2, Xk+1=Xk+αkPk,Pk+1=-Rk+1+βkPk, else, go to Step 4.If∥Rk+1∥<ε, stop, else set k←k+1 and go to Step 2.
(iv) Find the solutionhk to the following optimization problem:
(14)minh∈Rk+112hTTkh+γ0hTe1subjectto∥h∥≤Δ.
(v) Ifγk+1|〈ek+1,hk〉|<ε (here ek+1 represents the last column of identity matrix I), set
X
~
k
=
(
Q
0
,
Q
1
,
…
,
Q
k
)
(
h
k
⊗
I
) and then stop, else set k←k+1 and go to Step 2.The basic iteration route of Algorithm2 to solve problem (1) includes two cases: First, using CG method (Step 3) to compute the solution of problem (1) in feasible region. When the first case is failure, the solution of problem (1) in feasible region cannot be obtained by using CG method, and then the solution of problem (1) on the boundary can be obtained by solving the optimization problem (14). The properties about Algorithm 2 are given as follows.Theorem 3.
Assume that the sequences{Ri}, {Pi}, and {APiB} are generated by Algorithm 2; then the following equalities hold for all i≠j, 0≤i,j≤k:
(15)〈Ri,Rj〉=0,〈Pi,Rj〉=0,〈APiB,APjB〉=0.Proof.
Since〈A,B〉=〈B,A〉 holds for all matrices A and B, we only need to prove that the conclusion holds for all 0≤i<j≤k. Using induction and two steps are required.
Step1. Show that 〈Ri,Ri+1〉=0, 〈Pi,Ri+1〉=0 and 〈APiB,APi+1B〉=0 hold for all i=0,1,2,…k. We also use the principle of mathematical induction to prove these conclusions. When i=0, we have
(16)〈R0,R1〉=〈R0,R0+α0ATAP0BBT〉=〈R0,R0〉+〈R0,R0〉〈AP0B,AP0B〉〈R0,ATAP0BBT〉=〈R0,R0〉+〈R0,R0〉〈AP0B,AP0B〉〈AR0B,AP0B〉=0,〈P0,R1〉=〈P0,R0+α0ATAP0BBT〉=〈P0,R0〉+〈R0,R0〉〈AP0B,AP0B〉〈P0,ATAP0BBT〉=〈P0,R0〉+〈R0,R0〉〈AP0B,AP0B〉〈AP0B,AP0B〉=0,〈AP0B,AP1B〉=〈AP0B,A(-R1+β0P0)B〉=-〈AP0B,AR1B〉+〈R1,R1〉〈R0,R0〉〈AP0B,AP0B〉=-〈ATAP0BBT,R1〉+〈R1,R1〉〈R0,R0〉〈AP0B,AP0B〉=-〈AP0B,AP0B〉〈R0,R0〉〈R1-R0,R1〉+〈R1,R1〉〈R0,R0〉〈AP0B,AP0B〉=0.
Assume that conclusion holds for all i≤s(0<s<k); then
(17)〈Rs,Rs+1〉=〈Rs,Rs+αsATAPsBBT〉=〈Rs,Rs〉+〈Rs,Rs〉〈APsB,APsB〉〈Rs,ATAPsBBT〉=〈Rs,Rs〉+〈Rs,Rs〉〈APsB,APsB〉×〈(-Ps+βs-1Ps-1),ATAPsBBT〉=〈Rs,Rs〉+〈Rs,Rs〉〈APsB,APsB〉×(〈-Ps,ATAPsBBT〉+〈βs-1Ps-1,ATAPsBBT〉)=〈Rs,Rs〉+〈Rs,Rs〉〈APsB,APsB〉〈-Ps,ATAPsBBT〉=0,〈Ps,Rs+1〉=〈Ps,Rs+αsATAPsBBT〉=〈Ps,Rs〉+〈Rs,Rs〉〈APsB,APsB〉〈Ps,ATAPsBBT〉=〈-Rs+βs-1Ps-1,Rs〉+〈Rs,Rs〉〈APsB,APsB〉〈APsB,APsB〉=0,〈APsB,APs+1B〉=〈APsB,A(-Rs+1+βsPs)B〉=-〈APsB,ARs+1B〉+〈Rs+1,Rs+1〉〈Rs,Rs〉〈APsB,APsB〉=-〈ATAPsBBT,Rs+1〉+〈Rs+1,Rs+1〉〈Rs,Rs〉〈APsB,APsB〉=-〈APsB,APsB〉〈Rs,Rs〉〈Rs+1-Rs,Rs+1〉+〈Rs+1,Rs+1〉〈Rs,Rs〉〈APsB,APsB〉=0.
By the principle of induction, 〈Pi,Ri+1〉=0, 〈Ri,Ri+1〉=0, and 〈APiB,APi+1B〉=0 hold for all i=0,1,2,…k.
Step2. Assume that 〈Pi,Ri+l〉=0, 〈APiB,APi+lB〉=0, and 〈Ri,Ri+l〉=0 for all 0≤i≤k and 1<l<k, show that 〈Pi,Ri+l+1〉=0, 〈APiB,APi+l+1B〉=0, and 〈Ri,Ri+l+1〉=0. The proof is as follows:
(18)〈Pi,Ri+l+1〉=〈Pi,Ri+l+αi+lATAPi+lBBT〉=〈Pi,Ri+l〉+αi+l〈Pi,ATAPi+lBBT〉=0+αi+l〈APiB,APi+lB〉=0,〈APiB,APi+l+1B〉=〈Pi,ATAPi+l+1BBT〉=1αi+l〈Pi,Ri+l+1-Ri+l〉=0,〈Ri,Ri+l+1〉=〈Ri,Ri+l+αi+lATAPi+lBBT〉=αi+l〈Ri,ATAPi+lBBT〉=αi+l〈(-Pi+βi-1Pi-1),ATAPi+lBBT〉=αi+l(〈-APiB,APi+lB〉+βi-1〈APi-1B,APi+lB〉)=0.
From Steps 1 and 2, we have by principle induction that 〈Ri,Rj〉=0, 〈Pi,Rj〉=0, 〈APiB,APjB〉=0 hold for 0≤i<j≤k.Theorem 4.
Assume that the sequence{Qi} is generated by Algorithm 2; then the following equalities hold:
(19)〈Qi,Qj〉={1,i=j=0,1,2,…,0,i≠j,i,j=0,1,2,….Proof.
By the definition ofQi, we immediately know that 〈Qi,Qi〉=1(i=0,1,2,…). Similar to the proof of Theorem 3, we also use the principle of mathematical induction to prove this conclusion with the two following cases.
Step1. Show that 〈Qi,Qi+1〉=0 for all i=0,1,2,…k.
Wheni=0, we have
(20)〈Q0,Q1〉=1γ1〈Q0,ATAQ0BBT-δ0Q0〉=1γ1(〈Q0,ATAQ0BBT〉-〈Q0,δ0Q0〉)=1γ1(〈AQ0B,AQ0B〉-〈AQ0B,AQ0B〉〈Q0,Q0〉)=0.
Assume that conclusion holds for all i≤s(0<s<k); then
(21)〈Qs,Qs+1〉=1γs+1〈Qs,ATAQsBBT-δsQs-γsQs-1〉=1γs+1(〈Qs,ATAQsBBT〉-〈Qs,δsQs〉-〈Qs,γsQs-1〉)=1γs+1(〈AQsB,AQsB〉-〈AQsB,AQsB〉×〈Qs,Qs〉-γs〈Qs,Qs-1〉)=0.
By the principle of induction, 〈Qi,Qi+1〉=0 holds for all i=0,1,2,…k.
Step2. Assume that 〈Qi,Qi+l〉=0 for all 0≤i≤k and 1<l<k show that 〈Qi,Qi+l+1〉=0. The proof is as follows:
(22)〈Qi,Qi+l+1〉=1γi+l+1〈Qi,ATAQi+lBBT-δi+lQi+l-γi+lQi+l-1〉=1γi+l+1(〈AQiB,AQi+lB〉-〈Qi,δi+lQi+l〉-〈Qi,γi+lQi+l-1〉)=1γi+l+1(〈ATAQiBBT,Qi+l〉-0-0)=1γi+l+1〈-γi+1Qi+1-δiQi-γiQi-1,Qi+l〉=0.
From steps 1 and 2, we have by the principle of mathematical induction that 〈Qi,Qj〉=0 hold for all i,j=0,1,2,…k,i≠j.Theorem 5.
Assume that the sequences{γk}, {Tk}, and {Qi} are generated by Algorithm 2. Let
(23)X~=Q0h0+Q1h1+⋯+Qkhk=(Q0,Q1,…,Qk)(h⊗I),h=(h0,h1,…,hk)T∈Rk+1.
Then the following equality holds:
(24)12〈AX~B,AX~B〉-〈AX~B,C〉=12hTTkh+γ0hTe1,
where e1 represents the first column of identity matrix I.Proof.
By the definition ofTk and Qk(k=0,1,2,…), we have
(25)(ATAQ0BBT,ATAQ1BBT,…,ATAQkBBT)=(Q0,Q1,…,Qk)(Tk⊗I)+(0,…,0,γk+1Qk+1).
Hence, we have
(26)12〈AX~B,AX~B〉-〈AX~B,C〉=12〈X~,ATAX~BBT〉-〈X~,ATCBT〉=12〈(Q0h0+Q1h1+…+Qkhk)BBT(Q0,Q1,…,Qk)(h⊗I),ATA(Q0h0+Q1h1+⋯+Qkhk)BBT〉-〈(Q0,Q1,…,Qk)(h⊗I),γ0Q0〉=12〈(Q0,Q1,…,Qk)(h⊗I),ATAQ0BBTh0+ATAQ1BBTh1+⋯+ATAQkBBThk〉+γ0trace(h0Q0TQ0)=12〈(Q0,Q1,…,Qk)(h⊗I),(ATAQ0BBT,ATAQ1BBT,…,ATAQkBBT)×(h⊗I)(Q0,Q1,…,Qk)(h⊗I)(Q0,Q1⋯Qk)(h⊗I)〉+γ0h0=12〈(Q0,Q1,…,Qk)(h⊗I),(Q0,Q1,…,Qk)×[(Tk⊗I)+(0,…,0,γk+1Qk+1)]×(h⊗I)(Q0,Q1,…,Qk)(h⊗I)(Q0,Q1,…,Qk)(h⊗I)〉-γ0hTe1=12〈(Q0,Q1,…,Qk)(h⊗I),(Q0,Q1,…,Qk)(Tk⊗I)(h⊗I)(Q0,Q1,…,Qk)(h⊗I)〉-γ0hTe1=12trace[(hTTk⊗I)(Q0,Q1,…,Qk)T×(Q0,Q1,…,Qk)(h⊗I)(hTTk⊗I)(Q0,Q1,…,Qk)T]+γ0hTe1=12hTTkh+γ0hTe1.
So the equality (24) holds. In addition, from above equality, we havehTTkh=〈AX~B,AX~B〉 for allh∈Rk+1. SoTkis positive semi-definite. The proof is completed.Theorem 6.
Assume that the sequences{Qk}, {Rk}, {γk}, {δk}, {αk}, and {βi} are generated by Algorithm 2, then the following equalities hold for all k=0,1,2,…:
(27)Qk=(-1)kRk∥Rk∥,δk={1αk,k=0,1αk+βk-1αk-1,k>0,γk=βk-1αk-1.Proof.
(the proof of the first equality in (27)). By the definition of Qk and Rk, we have
(28)Qk=ak(ATA)k(ATCBT)(BBT)k+ak-1(ATA)k-1(ATCBT)(BBT)k-1+⋯+a0(ATCBT),Rk=(-1)kbk(ATA)k(ATCBT)(BBT)k+bk-1(ATA)k-1(ATCBT)(BBT)k-1+⋯+b0(ATCBT),
where ai, bi(i=0,1,2,…,k) are positive numbers. These equalities imply that Qk and Rk belong to the same space
(29)Kk=span{(ATA)k(ATCBT)(BBT)k,(ATA)k-1(ATCBT)(BBT)k-1,…,(ATCBT)}.
And furthermore we can have
(30)span{Qk-1,Qk-2,…,Q0}=Kk-1=span{Rk-1,Rk-2,…,R0}.
By Theorems 3 and 4, we have Qk⊥Kk-1 and Rk⊥Kk-1. Hence Qk and Rk must be linear correlation, so there exists a real number ck such that Qk=ckRk. Noting that ∥Qk∥=1, we have by (28) that Qk=(-1)kRk/∥Rk∥.
(The proof of the second equality in (27)). Noting that the first equality in (27) holds, then, when k=0, we have
(31)δ0=〈AQ0B,AQ0B〉=〈AR0B,AR0B〉〈R0,R0〉=〈AP0B,AP0B〉〈R0,R0〉=1α0.
When k>0, we have
(32)δk=〈AQkB,AQkB〉=〈ARkB,ARkB〉〈Rk,Rk〉=〈A(-Pk+βk-1Pk-1)B,A(-Pk+βk-1Pk-1)B〉〈Rk,Rk〉=1αk+βk-12〈APk-1B,APk-1B〉〈Rk,Rk〉=1αk+βk-1〈APk-1B,APk-1B〉〈Rk-1,Rk-1〉=1αk+βk-1αk-1.
(The proof of the third equality in (27)). By the definition of γk, we have
(33)γk2=〈tk,tk〉=〈ATAQk-1BBT-δk-1Qk-1-γk-1Qk-2,γk-1Qk〉=γk〈ATAQk-1BBT,Qk〉=-γk〈ATARk-1BBT,Rk〉(∥Rk-1∥∥Rk∥)=-γk〈ATA(-Pk-1+βk-2Pk-2)BBT,Rk〉(∥Rk-1∥∥Rk∥)=-γk〈(Rk-1-Rk)αk-1+(βk-2αk-2)(Rk-1-Rk-2),Rk〉×(∥Rk-1∥∥Rk∥)-1=γkβk-1αk-1.
Hence the third equality in (27) holds. The proof is completed.Remark 7.
This theorem shows the relationship between the sequences{Qk}, {Rk}, {γk}, {δk}, {αk}, and {βi} to lower down the cost of calculation.


## 3. The Main Results and Improvement of the Iteration Method

We will show that the solution of the problem (1) can be obtained within finite iteration steps in the absence of round-off errors. And we give the detail to solve the problem (14) in order to complete Algorithm 2. By discussing the characterization of the proposed iteration method, the further optimization method for the proposed iteration method is given at the end of this section.Theorem 8.
Assume that the sequences{Xk}, {Rk} are generated by Algorithm 2. Then the following equalities hold for all k=0,1,2,…:
(34)ATAXkBBT-ATCBT=Rk.Proof.
We use the principle of mathematical induction to prove this conclusion. Whenk=0, obviously, the conclusion holds. Assume that the conclusion holds for k-1; then
(35)ATAXkBBT-ATCBT=ATA(Xk-1+αk-1Pk-1)BBT-ATCBT=ATAXk-1BBT-ATCBT+αk-1ATAPk-1BBT=Rk-1+αk-1ATAPk-1BBT=Rk.
This implies that the conclusion holds for k. By the principle of mathematical induction, we know that the conclusion holds for all k=0,1,2,….Remark 9.
For Theorem3, the sequences R0,R1,R2,… are orthogonal to each other in the finite dimension matrix space Rn×n; it is certain that there exists a positive number k+1≤n2 such that Rk+1=0. So without the error of calculation, the first stopping criterion in the algorithm will perform with finite steps. From Theorem 8, we get ATAXk+1BBT-ATCBT=0. According to Theorem 1, when we set λ*=0, Xk+1 is a solution of the problem (3).Theorem 10.
Assume that the sequences{Qk}, {γk}, and {hk} are generated by Algorithm 2. Let
(36)X~k=Q0hk0+Q1hk1+⋯+Qkhkk,hk=(hk0,hk1,…,hkk).
Then, for all k=0,1,2,…, there exists a nonnegative number λk such that
(37)ATAX~kBBT+λkX~k-ATCBT=Qk+1γk+1hkk,λk(∥X~k∥-Δ)=0.Proof.
Assume thathk is the solution of optimization problem (14); then there exists a nonnegative number λk such that the following optimality Karush-Kuhn-Tucker (KKT) conditions are satisfied:
(38)(Tk+λkI)hk=-γ0e1,λk(∥hk∥-Δ)=0.
Noting that ∥X~k∥=∥hk∥ and the second equality in (38) hold, we know that the second equality in (37) holds.
Since the first equality in (38) can be rewritten as
(39)(Tk⊗I)(hk0I,hk1I,…,hkkI)T+λk(hk0I,hk1I,…,hkkI)T+(γ0I,0,…,0)T=0,
so we have
(40)(Q0,Q1,…,Qk)×[(Tk⊗I)(hk0I,hk1I,…,hkkI)T+λk(hk0I,hk1I,…,hkkI)T+(γ0I,0,…,0)T]=0.
Hence, we have
(41)ATAX~kBBT+λkX~k-ATCBT-γk+1hkkQk+1=0.
The proof is completed.Theorem 11.
Assume thatγ0,γ1,…,γk≠0, and γk+1=0. Then X~k=Q0hk0+Q1hk1+⋯+Qkhkk is the solution of the problem (1).Proof.
Sinceγk+1=0 and X~k=Q0hk0+Q1hk1+⋯+Qkhkk, we have by Theorem 10 that
(42)ATAX~kBBT+λkX~k=ATCBT,λk(∥X~k∥-Δ)=0,
which implies that X~k is the solution of the problem (1).Remark 12.
According to Theorem4, the sequences Q0,Q1,Q2,… are orthogonal each other in the finite dimension matrix space Rn×n; it is certain that there exists a positive number k≤n2 such that Qk=0. Since tk=γkQk=0, then γk=〈tk,tk〉=0. So without the error of calculation, the second stopping criterion in the algorithm also performs with finite steps.Remark 13.
According to Remarks9 and 12, we have that, without the error of calculation, a desired solution can be obtained with finitely iterative step by Algorithm 2.Theorem 14.
The solutionhk of the problem (14) obtained by Algorithm 2 is on the boundary. In other words, hk is the solution of the following optimization problem:
(43)minh∈Rk+112hTTkh+hT(γ0e1)subjectto∥h∥=Δ.Proof.
Assuming that the solutionhk of the problem (14) obtained by Algorithm 2 is inside the boundary, we have by (38) that Tkhk=-γ0e1. By Theorem 5, we know Tk is a positive semidefinite matrix. If Tk is positive definite, then hk=-Tk-1(γ0e1) with ∥hk∥<Δ is a unique solution of the problem (14). Hence, we have by Theorem 5 that X=(Q0,Q1,…,Qk)(hk⊗I) with ∥X∥=∥hk∥<Δ is a unique solution of the problem (1). In this case, the step of solving the problem (14) in Algorithm 2 cannot be implemented. If Tk is positive semidefinite and not positive definite, then there exists a matrix Z such that Tk(hk+Z)=-γ0e1 and ∥hk+Z∥=Δ which implies that hk+Z is a solution to the problem (1) on the boundary. This contradicts our assumption.Now we use the following Algorithm15, which was proposed by More and Sorensen in paper [11], to solve the problem (43).Algorithm 15.
(I) Let a suitable starting valueλk0 and Δ>0 be given.
(II) Fori=0,1,… until convergence.(a)
FactorizeTk+λkiI=QΛQT, where Q and Λ are unit bidiagonal and diagonal matrix, respectively.(b)
SolveQΛQTh=-γ0e1.(c)
SolveQw=h.(d)
Setλki+1=λki+1+((∥h∥-Δ)/Δ)(∥h∥2/〈w,Λ-1w〉).
In the implementation of Algorithm15, the initial secular λk0 can be chosen by the following principles: If ∥hk(λk-1)∥≥Δ, let λk0=λk-1; else let λk0=0, where λk-1 is obtained by the (k-1)th iterative steps of Algorithm 2. The stopping criteria can be used as |λki+1-λki|≤ε, where ε is a small tolerance.By fully using the result of Theorem6, Algorithm 2 can be optimized as in Algorithm 16.Algorithm 16.
(i) Given matricesX0=0, Q-1=0 and a small tolerance ε>0.
ComputingR0=-ATCBT,t0=-ATCBT.
Setγ0=∥R0∥,P0=-R0,T-1=[], and k←0.
(ii) IfAPkB≠0, compute
(44)Qk=(-1)kRk∥Rk∥,αk=∥Rk∥2∥APkB∥2,Rk+1=Rk+αkATAPkBBT,βk=∥Rk+1∥2∥Rk∥2,δk={1αk,k=0,1αk+βk-1αk-1,k>0,γk+1=βkαk,Tk=[Tk-1ΓkΓkTδk],
where Γk=(0,…,0,γk)T∈Rk.
Else, computingQk=tk/γk (the first one Qk=(-1)kRk/∥Rk∥),
(45)δk=∥AQkB∥2,tk+1=ATAQkBBT-δkQk-γkQk-1,γk+1=∥tk+1∥,Tk=[Tk-1ΓkΓkTδk].
(iii) If∥Xk+1+αk+1Pk+1∥≤Δ, computing Xk+1=Xk+αkPk, Pk+1=-Rk+βkPk.
If∥Rk+1∥≤ε, stop. Else, setting k←k+1 and go to Step 2.
Else, go to Step 4.
(iv) Using Algorithm15 to compute the solution hk of the problem (43).
(v) Ifγk+1|〈ek+1,hk〉|<ε, setting X~k=(Q0,Q1,…,Qk)(hk⊗I), then stop.
Else, settingk←k+1 and go to step 2.


## 4. Numerical Experiments

In this section, we present numerical examples to illustrate the availability and the real application of the proposed iteration method. All tests are performed using MATLAB 7.1 with a 32-bit Windows XP operating system. Our experiments are performed on an FOUNFER computer of mode E520 with 2.8 GHz CPU and 3.25 G RAM. Because of the error of calculation, the iteration will not stop with finite steps. Hence, we regard the approximation solutionXk as the solution of problem (1) if the t(k)≤10-10, where
(46)t(k)={∥Rk+1∥F,∥Xk+αkPk∥F<Δγk+1|〈ek+1,hk〉|,∥Xk∥F=Δ.Example 17.
Given the matricesA,B,C as follows:(47)
A
=
(
2.1133
-
2.3711
-
1.5741
-
1.0450
-
2.3090
2.0660
-
1.4522
3.6855
-
3.5974
0.6188
-
1.4764
2.8853
3.1969
0.0771
3.9125
-
2.0298
3.0998
-
1.2121
-
0.4192
-
2.5485
0.2237
1.5758
-
1.8922
0.8742
0.7832
0.6938
-
4.2669
0.7283
-
2.3402
2.1912
-
1.0682
0.1450
3.9923
0.7929
-
2.2801
)
,
B
=
(
1.8289
-
0.0255
-
0.1382
2.8240
0.8011
1.2860
0.4890
-
0.1238
-
0.0766
-
0.1251
-
0.8981
1.5505
0.2457
0.4984
2.4399
2.7879
-
1.5069
-
2.2386
-
1.6603
-
0.3753
0.6156
-
1.2884
-
1.4745
-
1.9062
2.2008
0.1921
-
3.3952
-
4.9563
-
2.0966
-
0.4614
)
,
C
=
(
-
141.9226
-
7.5851
25.3523
34.6519
236.6062
80.9125
167.3048
8.9417
-
29.8864
-
40.8493
-
278.9222
-
95.3833
31.8605
1.7028
-
5.6914
-
7.7791
-
53.1162
-
18.1642
-
46.6795
-
2.4948
8.3386
11.3973
77.8216
26.6127
44.4786
2.3772
-
7.9454
-
10.8599
-
74.1526
-
25.3580
)
.WhenΔ=40, using Algorithm 16 and iterate 43 steps, we obtain the approximation solution(48)X43=(0.77190.56780.36930.82180.46081.52392.74873.60391.28703.09186.956610.166912.35336.347811.0186-1.1825-1.4474-1.5816-1.1406-1.47638.24919.738810.49058.01809.90194.91487.04608.55584.50007.66355.98929.669012.37305.255110.8217).WhenΔ=10, using the Algorithm 16 and iterative 23 steps, we obtain the approximation solution(49)X23=(0.28090.41830.94990.20190.8336-0.2275-0.1577-0.17260.08370.24340.14170.75322.26830.84223.1792-0.1560-0.2950-0.7358-0.1975-0.78441.01181.87275.14841.36855.51550.27030.74961.96220.63752.4602-0.17690.29151.04170.57702.0664).Given a nonnegative real numberΔ=1000, with iterate 45 steps, we obtain the approximation solution(50)X45=(0.43490.43490.43490.43490.43493.17803.17803.17803.17803.178011.237311.237311.237311.237311.2373-1.4841-1.4841-1.4841-1.4841-1.48419.95479.95479.95479.95479.95477.82127.82127.82127.82127.821211.114311.114311.114311.114311.1143).Example 18.
We work with a 2D first-kind Fredholm integral equation of the generic form(51)∬01κ(x-x′)ω(y-y′)f(x′,y′)dx′dy′=g(x,y),
where κ and ω are function. Based on [12], we have that the discretization of the problem (51) leads to the linear relation AFA-T=G between the discrete solution F and the discrete data G, where
(52)Aik=m-1κ(xi-xk′),A-jl=n-1ω(yj-yl′)Fkl=f(xk′,yk′),Gij=g(xi,yj),i,k=1,2,…,m,j,l=1,2,…,n.
An example of such problem is image denoising with a Gaussian point spread function:(53)κ(t)=ω(t)=12πσexp(-12(tσ)2),
which is used as a model for out-of-focus as well as atmospheric turbulence blur [13]. In Figure 1(a), the original image is the standard test image of Lena with size 256×256, which is also a 256×256 matrix F′. After the image was blurred by Gaussian kernel (53) with σ=0.01, we get Figure 1(b); that is the matrix G=AF′A-T. In image denoising, our target is to get the solution of AFA-T=G. Tikhonov regularization is needed to treat this problem in order to control the effect of the noise on the solution. As we have said in Section 1, Tikhonov regularization is equivalent to over the norm inequality constraint matrix equation
(54)minF∈Rm×n12∥AFA-T-G∥2subjectto∥F∥≤Δ.(a) The original image. (b) The noisy image. (c) The recovered image.
(a)
(b)
(c)Based on [7], Δ represents the energy of the target image, so we get Δ=∥F′∥. Solving the above problem by Algorithm 16, we get the recovered image F* in Figure 1(c). It means our algorithm is suitable for image denoising.



---

*Source: 781276-2013-07-22.xml*