A New Iterative Algorithm for Solving a Class of Matrix Nearness Problem

Abstract

Based on the alternating projection algorithm, which was proposed by Von Neumann to treat the problem of finding the projection of a given point onto the intersection of two closed subspaces, we propose a new iterative algorithm to solve the matrix nearness problem associated with the matrix equations 𝐴𝑋𝐵=𝐸, 𝐶𝑋𝐷=𝐹, which arises frequently in experimental design. If we choose the initial iterative matrix 𝑋₀=0, the least Frobenius norm solution of these matrix equations is obtained. Numerical examples show that the new algorithm is feasible and effective.

1. Introduction

Denoted by 𝑅^𝑚×𝑛 be the set of 𝑚×𝑛 real matrices, 𝐴^𝑇 and 𝐴⁺ be the transpose and Moore-Penrose generalized inverse of the matrix 𝐴, respectively. For 𝐴,𝐵∈𝑅^𝑚×𝑛,  ⟨𝐴,𝐵⟩=trace(𝐵^𝑇𝐴) denotes the inner product of the matrix 𝐴 and 𝐵. The induced norm is the so-called Frobenius norm, that is, ‖𝐴‖=⟨𝐴,𝐴⟩^1/2; then 𝑅^𝑚×𝑛 is a real Hilbert space. Let 𝑀 be a closed convex subset in a real Hilbert space 𝐻 and 𝑢 be a point in 𝐻; then the point in 𝑀 nearest to 𝑢 is called the projection of 𝑢 onto 𝑀 and denoted by 𝑃_𝑀(𝑢); that is to say, 𝑃_𝑀(𝑢) is the solution of the following minimization problem (see [1, 2]) min_𝑥∈𝑀‖𝑥−𝑢‖,(1) that is, ‖‖𝑃_𝑀(𝑢)−𝑢‖‖=min_𝑥∈𝑀‖𝑥−𝑢‖.(2)

The problem of finding a nearness matrix 𝑋 in a constraint matrix set to a given matrix 𝑋 is called the matrix nearness problem. Because the preliminary estimation 𝑋 is frequently obtained from experiments, it may not satisfy the given restrictions. Hence it is necessary to find a nearness matrix 𝑋 in this constraint matrix set to replace the estimation 𝑋 [3]. In the area of structure design, finite element model updating and control theory, and so forth, the matrix set is always the (constraint) solution set or the least square (constraint) solution set of some matrix equations [4–6]. Thus, the problem mentioned above is also called the matrix nearness problem associated with the matrix equation. Recently, there are many discussions on the matrix nearness problem associated with some matrix equations. For instance, see [4, 6–14].

In this paper, we consider the following problem.

Problem 1. Given matrices  𝐴∈𝑅^𝑝×𝑛,𝐵∈𝑅^𝑛×𝑞,𝐶∈𝑅^𝑠×𝑛,𝐷∈𝑅^𝑛×𝑡,𝐸∈𝑅^𝑝×𝑞,𝐹∈𝑅^𝑠×𝑡, and 𝑋∈𝑅^𝑛×𝑛, find 𝑋∈Ω such that ‖‖𝑋−𝑋‖‖=min_𝑋∈Ω‖‖𝑋−𝑋‖‖,(3) where Ω=𝑋∈𝑅^𝑛×𝑛∣𝐴𝑋𝐵=𝐸,𝐶𝑋𝐷=𝐹.(4) Obviously, Ω is the solution set of the matrix equations 𝐴𝑋𝐵=𝐸,𝐶𝑋𝐷=𝐹,(5) and 𝑋 is the optimal approximate solution of (5) to the given matrix 𝑋. In particular, if 𝑋=0, then the solution 𝑋 of Problem 1 is just the least Frobenius norm solution of the matrix equations (5). It is easy to verify that Ω is convex set; then the solution of Problem 1 is unique.

The matrix equations (5) and its matrix nearness problem have been extensively studied for the past 40 or more years. Wang [15] and Navarra et al. [16] gave some conditions for the existence of a solution and some representations of the general common solution to (5). By the projection theorem and matrix decompositions, Liao et al. [6] gave an analytical expression of the optimal approximate least square symmetric solution of (5). However, these direct methods may be less efficient for the large sparse coefficient matrices due to the limited storages and the speeds of the computers. Therefore, iterative methods for solving the matrix equations (5) have attracted much interests recently. Peng et al. [11] and Chen et al. [7] proposed some iterative methods to compute the symmetric solutions and optimal approximate symmetric solution of (5). An efficient iterative method was presented to solve the matrix nearness Problem 1 associated with the matrix equations (5) in [13]. Ding et al. [17] considered the unique solution of the matrix equations (5) and used gradient-based iterative algorithm to compute the unique solution. The (least square) solution and the optimal approximate (least square) solution of (5), which is constrained as bisymmetric, reflexive, generalized reflexive, generalized centrosymmetric, were studied in [7–10, 12].

The alternating projection algorithm dates back to von Neumann [18], who treated the problem of finding the projection of a given point onto the intersection of two closed subspaces. Later, Bauschke and Borwein [1] extended the analysis of Von Neumann’s alternating projection scheme to the case of two closed affine subspaces. There are many variations and extensions of the alternating projection algorithm, and we can use them to find the projection of a given point onto the intersection of 𝑘(𝑘≥2) closed subspaces [19] and 𝑘(𝑘≥2) closed convex sets [20, 21]. For a complete discussion on the alternation projection algorithm see Deutsch [2].

In this paper, we propose a new algorithm to solve Problem 1. We state Problem 1 as the minimization of a convex quadratic function over the intersection of two closed affine subspaces in the vector space 𝑅^𝑛×𝑛. From this point of view, Problem 1 can be solved by the alternating projection algorithm. If we choose the initial iterative matrix 𝑋₀=0, the least Frobenius norm solution of the matrix equations 𝐴𝑋𝐵=𝐸,𝐶𝑋𝐷=𝐹 is obtained. In the end, we use some numerical examples to show that the new algorithm is faster and lower computational cost for each step than the algorithm proposed by Sheng and Chen [13] to solve Problem 1. Especially, the CPU time and iteration steps of our algorithm increase slowly as the dimension of the matrix is increasing; so our algorithm is suitable for large-scale problems.

2. Alternating Projection Algorithm for Solving Problem 1

In this section, we apply the alternating projection algorithm to solve Problem 1. We begin with two lemmas.

Lemma 1 (see [1, Theorem  4.1]). Let 𝑀₁=𝑎+𝑀₁,𝑀₂=𝑏+𝑀₂ be closed affine subspaces in a Hilbert space 𝐻 and u be a point in 𝐻. Here, 𝑀₁ and 𝑀₂ are closed subspaces and 𝑎,𝑏∈𝐻. If 𝑀₁∩𝑀₂≠∅, then the sequences {𝑥_𝑘} and {𝑦_𝑘} generated by the alternating projection algorithm 𝑥_𝑘+1𝑥₀=𝑢,𝑦_𝑘+1=𝑃_𝑀1𝑥_𝑘,=𝑃_𝑀2𝑦_𝑘+1,𝑘=0,1,2,…(6) both converge to the projection of the point 𝑢 onto the intersection of 𝑀₁ and 𝑀₂, that is, 𝑥_𝑘⟶𝑃_{𝑀1∩𝑀2}(𝑢),𝑦_𝑘⟶𝑃_{𝑀1∩𝑀2}(𝑢),𝑘⟶+∞.(7)

Lemma 2 (see [22, Theorem  9.3.1]). Let 𝐴∈𝑅^𝑝×𝑛,𝐵∈𝑅^𝑛×𝑞 and 𝐸∈𝑅^𝑝×𝑞, be known matrices. Then the matrix equation 𝐴𝑋𝐵=𝐸 has a solution if and only if 𝐴𝐴⁺𝐸𝐵⁺𝐵=𝐸,(8) and the representation of the solution is 𝑋=𝐴⁺𝐸𝐵⁺+𝑌−𝐴⁺𝐴𝑌𝐵𝐵⁺,(9) where 𝑌∈𝑅^𝑛×𝑛 is arbitrary.

Lemma 3 (see [22, Theorem  9.3.2]). Given 𝑍∈𝑅^𝑛×𝑛, set ℜ=𝑋∈𝑅^𝑛×𝑛∣𝐴𝑋𝐵=𝐸,𝐴∈𝑅^𝑝×𝑛,𝐵∈𝑅^𝑛×𝑞,𝐸∈𝑅^𝑝×𝑞,(10) then the solution 𝑋 of the following problem min_𝑋∈ℜ‖𝑋−𝑍‖(11) is 𝑋=𝑍+𝐴⁺(𝐸−𝐴𝑍𝐵)𝐵⁺,(12) that is, ‖‖𝑋−𝑍‖‖=min_𝑋∈ℜ‖𝑋−𝑍‖.(13)

Now we begin to use the alternating projection algorithm (6) to solve Problem 1. Firstly, we define two sets ΩΩ₁=𝑋∈𝑅^𝑛×𝑛∣𝐴𝑋𝐵=𝐸₂,=𝑋∈𝑅^𝑛×𝑛∣𝐶𝑋𝐷=𝐹.(14) It is easy to know that Ω=Ω₁∩Ω₂, and if the set Ω is nonempty, then Ω=Ω₁∩Ω₂≠∅.(15) And by Lemma 2, the sets Ω₁ and Ω₂ can be equivalently written as Ω₁=={𝑋𝑋∣=𝐴𝐴𝑋₊𝐵=𝐸}=𝐴₊𝐸𝐵⁺+𝑌−𝐴⁺𝐴𝑌𝐵𝐵⁺∣𝑌∈𝑅^𝑛×𝑛Ω₂𝐸𝐵⁺+𝑌−𝐴⁺𝐴𝑌𝐵𝐵⁺∣𝑌∈𝑅^𝑛×𝑛,=={𝑋𝑋∣=𝐶𝐶𝑋₊𝐷=𝐹}=𝐶₊𝐹𝐷⁺+𝑌−𝐶⁺𝐶𝑌𝐷𝐷⁺∣𝑌∈𝑅^𝑛×𝑛𝐹𝐷⁺+𝑌−𝐶⁺𝐶𝑌𝐷𝐷⁺∣𝑌∈𝑅^𝑛×𝑛.(16) Hence, Ω₁ and Ω₂ are closed affine subspaces.

After defining the sets Ω₁ and Ω₂, Problem 1 can be rewritten as finding 𝑋∈Ω=Ω₁∩Ω₂, such that ‖‖𝑋−𝑋‖‖=min_{𝑋∈Ω1∩Ω2}‖‖𝑋−𝑋‖‖.(17) Noting the equalities (17) and (2), it is easy to find that 𝑋=𝑃_Ω1∩Ω2𝑋.(18) Therefore, Problem 1 can be converted equivalently into finding the projection 𝑃_Ω1∩Ω2(𝑋) of a given matrix 𝑋 onto the intersection set Ω₁∩Ω₂. Now we will use alternating projection algorithm (6) to compute the projection 𝑃_Ω1∩Ω2(𝑋). Consequently, we can get the solution 𝑋 of Problem 1.

By (6) we can see that the key problems to realize the alternating projection algorithm (6) are how to compute the projections 𝑃_Ω1(𝑍),𝑃_Ω2(𝑍) of a matrix 𝑍 onto Ω₁ and Ω₂, respectively. Such problems are perfectly solvable in the following theorems.

Theorem 1. Suppose that the set Ω₁ is nonempty. For a given 𝑛×𝑛 matrix 𝑍, we have 𝑃_Ω1(𝑍)=𝑍+𝐴⁺(𝐸−𝐴𝑍𝐵)𝐵⁺.(19)

Proof. By (1) and (2), we know that the projection 𝑃_Ω1(𝑍) is the solution of the following minimization problem min_𝑋∈Ω1‖𝑋−𝑍‖,(20) and according to Lemma 3 we know that the solution of the minimization problem (20) is 𝑍+𝐴⁺(𝐸−𝐴𝑍𝐵)𝐵⁺. Hence, 𝑃_Ω1(𝑍)=𝑍+𝐴⁺(𝐸−𝐴𝑍𝐵)𝐵⁺.(21)

Theorem 2. Suppose that the set  Ω₂ is nonempty. For a given 𝑛×𝑛 matrix 𝑍, we have 𝑃_Ω2(𝑍)=𝑍+𝐶⁺(𝐹−𝐶𝑍𝐷)𝐷⁺.(22)

Proof. The proof is similar to that of Theorem 1 and is omitted here.

By the alternation projection algorithm (6) and Theorems 1 and 2, we get a new algorithm to solve Problem 1 which can be stated as follows.

Algorithm 1. One has(1)set 𝐴=𝐴⁺,𝐵=𝐵⁺,𝐶=𝐶⁺,𝐷=𝐷⁺;(2)set 𝑋₀=𝑋;(3)for 𝑘=0,1,2,3,…𝑋𝑌_𝑘+1=𝑃_Ω1𝑋_𝑘=𝑋_𝑘+𝐴𝐸−𝐴𝑋_𝑘𝐵𝐵,_𝑘+1=𝑃_Ω2𝑌_𝑘+1=𝑌_𝑘+1+𝐶𝐹−𝐶𝑌_𝑘+1𝐷𝐷,(23) end.
By Lemma 1 and (15) and (16), we get the convergence theorem for Algorithm 1.

Theorem 3. If the set Ω is nonempty, then the matrix sequences {𝑋_𝑘} and {𝑌_𝑘} generated by Algorithm 1 both converge to the projection 𝑃_Ω1∩Ω2(𝑋) of  𝑋 onto the intersection of Ω₁ and Ω₂, that is, 𝑋_𝑘⟶𝑃_Ω1∩Ω2𝑋,𝑌_𝑘⟶𝑃_Ω1∩Ω2𝑋,𝑘⟶+∞.(24)

Proof. If the set Ω is nonempty, by (15) we have Ω₁∩Ω₂≠∅.(25) And noting (16), we know that the sets Ω₁ and Ω₂ are closed affine subspaces in Hilbert space 𝑅^𝑚×𝑛. Hence, by Lemma 1 we derive that the matrix sequences {𝑋_𝑘} and {𝑌_k} generated by Algorithm 1 both converge to the projection 𝑃_Ω1∩Ω2(𝑋) of 𝑋 onto the intersection of Ω₁ and Ω₂, that is, 𝑋_𝑘⟶𝑃_Ω1∩Ω2𝑋,𝑌_𝑘⟶𝑃_Ω1∩Ω2𝑋,𝑘⟶+∞.(26)

Combining Theorem 3 and the equalities (18) and (17), we have the following.

Theorem 4. If the set Ω is nonempty, then the matrix sequence {𝑋_𝑘} and {𝑌_𝑘} generated by Algorithm 1 both converge to the unique solution of Problem 1. Moreover, if the initial matrix 𝑋₀=𝑋=0, then the matrix sequence {𝑋_𝑘} and {𝑌_𝑘} both converge to the least Frobenius norm solution of the matrix equations 𝐴𝑋𝐵=𝐸,𝐶𝑋𝐷=𝐹.

3. Numerical Experiments

In this section, we give some numerical examples to illustrate that the new algorithm is feasible and effective to solve Problem 1. All programs are written in MATLAB 7.8. We denote Error=‖𝐸−𝐴𝑋𝐵‖+‖𝐹−𝐶𝑋𝐷‖,(27) and use the practical stopping criterion Error≤1.0×10⁻¹⁰.

Example 1. Consider Problem 1 with 𝐹𝐸𝐷==𝐴𝐵=𝐶⎛⎜⎜⎜⎜⎝⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝1==33−−⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝1541=100−1851066⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝5331⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝7557⎛⎜⎜⎜⎜⎝255−−−1−5371391217733115−7−−1455345197159611−−−−8008−−495343600321871−−33553754−−99154111−−−31−53718511391661331157557−55222953217792−−−−−−−−94653431971921871−11558008445137600155−−−1557753713911115−−⎞⎟⎟⎟⎟⎠1541−12⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠34518517766,967,,75575555⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠5⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠⎞⎟⎟⎟⎟⎠,,.(28) Here we use ones(𝑛) and zeros(𝑛) to stand for 𝑛×𝑛 matrix of ones and zeros. It is easy to verify that 𝑋=ones(5) is a solution of the matrix equations (5); that is to say, the set Ω is nonempty. Therefore we can use Algorithm 1 to solve Problem 1.

Let 𝑋₀=𝑋=zeros(5). After 5 iterations of Algorithm 1, we get the optimal approximate solution 𝑋≈𝑋₅=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝−00000....0069.611048519912317700000.....88997112988780158073−−−00000..49...52283589940613236700000.....9999966699399887980500000.....78899811290877081387⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠,(29) which is also the least Frobenius norm solution of the matrix equations (5), and its residual errorError≈‖‖𝐸−𝐴𝑋₅𝐵‖‖+‖‖𝐹−𝐶𝑋₅𝐷‖‖=2.78×10⁻¹¹.(30) By concrete computations, we know that the distance from 𝑋 to the solution set Ω is min_𝑋∈Ω‖‖𝑋−𝑋‖‖=‖‖𝑋−𝑋‖‖≈‖‖𝑋₅−𝑋‖‖=3.9057.(31)

Let 𝑋₀=𝑋=ones(5). After 6 iterations of Algorithm 1, we get the optimal approximate solution 𝑋≈𝑋₆=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝54421.....7995243560694317286311111.....8772072283458877981176631.....2110204646166233387011111....111.044220750068259211111....8772.372387770820879⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠,(32) and its residual error Error≈‖‖𝐸−𝐴𝑋₆𝐵‖‖+‖‖𝐹−𝐶𝑋₆𝐷‖‖=3.89×10⁻¹¹.(33) By concrete computations, we know that the distance from 𝑋 to the solution set Ω is min_𝑋∈Ω‖‖𝑋−𝑋‖‖=‖‖𝑋−𝑋‖‖≈‖‖𝑋₆−𝑋‖‖=16.0902.(34)

Example 1 shows that Algorithm 1 is feasible to solve Problem 1.

Example 2. Consider Problem 1 with 𝐷𝐸𝐴===r𝐴arna∗dndo(𝑛n(,1e1s020(0𝑛,))𝑛,∗),𝐵,𝐹𝐵=𝐶=𝐶=r∗anroadnn(ed𝑛s,((71𝑛05),0∗𝑛)𝐷),,,(35) where rand(𝑠,𝑡) stand for 𝑠×𝑡 random matrix. Let the given matrix 𝑋=zeros(𝑛). It is easy to verify that 𝑋=ones(𝑛) is the solution of the matrix equations (5); that is, the set Ω is nonempty; therefore, we can use Algorithm 1 and the following algorithm proposed by Sheng and Chen [13] to solve Problem 1.

Algorithm 2. One has(1)input 𝐴,𝐵,𝐶,𝐷,𝐸,𝐹 and 𝑋₀;(2)calculate 𝑅𝑟₀=𝐸−𝐴𝑋₀𝐵,₀𝑃=𝐹−𝐶𝑋₀𝐷,₀𝑄₀=𝐴^𝑇𝑅₀𝐵^𝑇,=𝐶^𝑇𝑟₀𝐷^𝑇;(36)(3)if 𝑅_𝑘=0,𝑟₀=0, then stop; else, 𝑘∶=𝑘+1; (4)calculate 𝑋_𝑘=𝑋_𝑘−1+‖𝑅^𝑘−1‖²+‖𝑟_𝑘−1‖²‖𝑃_𝑘−1+𝑄_𝑘−1‖²𝑅_𝑘𝑃_𝑘−1+𝑄_𝑘−1,𝑟=𝐸−𝐴𝑋_𝑘𝐵,_𝑘𝑃=𝐹−𝐶𝑋_𝑘𝐷,_𝑘=𝐴^𝑇𝑅_𝑘𝐵^𝑇+‖𝑅^𝑘‖²+‖𝑟_𝑘‖²‖𝑅_𝑘−1‖²+‖𝑟_𝑘−1‖²𝑄_𝑘,=𝐶^𝑇𝑟_𝑘𝐷^𝑇+‖𝑅^𝑘‖²+‖𝑟_𝑘‖²‖𝑅_𝑘−1‖²+‖𝑟_𝑘−1‖²𝑄_𝑘−1;(37)(5)go to step 3.

It is easy to see that Algorithm 1 has lower computational cost for each step in the comparison with Algorithm 2. Experiments show that Algorithm 1 and Algorithm 2 are feasible to solve Problem 1. We list the iteration steps (denoted by IT), CPU time (denoted by CPU), residual error (denoted by ERR), and the distance ‖𝑋_IT−𝑋‖ (denoted by DIS) in Table 1.

Table 1 

From Table 1, we can see that Algorithm 1 outperforms Algorithm 2 in iteration step and CPU time. Therefore our algorithm is faster than the algorithm proposed by Sheng and Chen [13]. Especially, the CPU time and iteration steps of our algorithm increase slowly as the dimension 𝑛 is increasing; so our algorithm is suitable for large-scale problems.

4. Conclusion

The alternating projection algorithm dates back to von Neumann [18], who treated the problem of finding the projection of a given point onto the intersection of two closed subspace. In this paper, we first apply the alternating projection algorithm to solve Problem 1, which occurs frequently in experimental design [23]. If we choose the initial matrix 𝑋₀=0, the least Frobenius norm solution of the matrix equations 𝐴𝑋𝐵=𝐸,𝐶𝑋𝐷=𝐹 can be obtained. Numerical examples show that the new algorithm is faster and lower computational cost for each step than the algorithm proposed by Sheng and Chen [13] to solve Problem 1. Especially, the CPU time and iteration steps of the new algorithm increase slowly as the matrix’s dimension is increasing; so the alternating projection algorithm is suitable for large-scale matrix nearness problems.

Acknowledgments

The work was supported in part by National Science Foundation of China (11101100; 10861005) and Natural Science Foundation of Guangxi Province (0991238; 2011GXNSFA018138).