A New Method for the Bisymmetric Minimum Norm Solution of the Consistent Matrix Equations 𝐴₁𝑋𝐵₁=𝐶₁, 𝐴₂𝑋𝐵₂=𝐶₂

Abstract

We propose a new iterative method to find the bisymmetric minimum norm solution of a pair of consistent matrix equations 𝐴₁𝑋𝐵₁=𝐶₁, 𝐴₂𝑋𝐵₂=𝐶₂. The algorithm can obtain the bisymmetric solution with minimum Frobenius norm in finite iteration steps in the absence of round-off errors. Our algorithm is faster and more stable than Algorithm 2.1 by Cai et al. (2010).

1. Introduction

Let 𝐑^𝑚×𝑛 denote the set of 𝑚×𝑛 real matrices. A matrix 𝑋=(𝑥_𝑖𝑗)∈𝐑^𝑛×𝑛 is said to be bisymmetric if 𝑥_𝑖𝑗=𝑥_𝑗𝑖=𝑥_{𝑛−𝑖+1,𝑛−𝑗+1} for all 1≤𝑖, 𝑗≤𝑛. Let 𝐁𝐒𝐑^𝑛×𝑛 denote 𝑛×𝑛 real bisymmetric matrices. For any 𝑋∈𝐑^𝑚×𝑛, 𝑋^𝑇, tr(𝑋), ‖𝑋‖, and ‖𝑋‖₂ represent the transpose, trace, Frobenius norm, and Euclidean norm of 𝑋, respectively. The symbol vec(⋅) stands for the vec operator; that is, for 𝑋=(𝑥₁,𝑥₂,…,𝑥_𝑛)∈𝐑^𝑚×𝑛, where 𝑥_𝑖(𝑖=1,2,…,𝑛) denotes the 𝑖th column of 𝑋, vec(𝑋)=(𝑥^𝑇₁,𝑥^𝑇₂,…,𝑥^𝑇_𝑛)^𝑇. Let mat(⋅) represent the inverse operation of vec operator. In the vector space 𝐑^𝑚×𝑛, we define the inner product as ⟨𝑋,𝑌⟩=tr(𝑌^𝑇𝑋) for all 𝑋,𝑌∈𝐑^𝑚×𝑛. Two matrices 𝑋 and 𝑌 are said to be orthogonal if ⟨𝑋,𝑌⟩=0. Let 𝑆_𝑛=(𝑒_𝑛,𝑒_𝑛−1,…,𝑒₁) denote the 𝑛×𝑛 reverse unit matrix where 𝑒_𝑖(𝑖=1,2,…,𝑛) is the 𝑖th column of 𝑛×𝑛 unit matrix 𝐼_𝑛; then 𝑆^𝑇_𝑛=𝑆_𝑛,𝑆²_𝑛=𝐼_𝑛.

In this paper, we discuss the following consistent matrix equations: 𝐴₁𝑋𝐵₁=𝐶₁,𝐴₂𝑋𝐵₂=𝐶₂,𝑋∈𝐁𝐒𝐑^𝑛×𝑛, (1) where 𝐴₁∈𝐑^𝑝1×𝑛, 𝐵₁∈𝐑^𝑛×𝑞1, 𝐶₁∈𝐑^𝑝1×𝑞1, 𝐴₂∈𝐑^𝑝2×𝑛, 𝐵₂∈𝐑^𝑛×𝑞2, and 𝐶₂∈𝐑^𝑝2×𝑞2 are given matrices, and 𝑋∈𝐁𝐒𝐑^𝑛×𝑛 is unknown bisymmetric matrix to be found.

Research on solving a pair of matrix equations 𝐴₁𝑋𝐵₁=𝐶₁, 𝐴₂𝑋𝐵₂=𝐶₂ has been actively ongoing for the past 30 or more years (see details in [1–6]). Besides the works on finding the common solutions to the matrix equations 𝐴₁𝑋𝐵₁=𝐶₁, 𝐴₂𝑋𝐵₂=𝐶₂, there are some valuable efforts on solving a pair of the matrix equations with certain linear constraints on solution. For instance, Khatri and Mitra [7] derived the Hermitian solution of the consistent matrix equations 𝐴𝑋=𝐶, 𝑋𝐵=𝐷. Deng et al. [8] studied the consistent conditions and the general expressions about the Hermitian solutions of the matrix equations (𝐴𝑋,𝑋𝐵)=(𝐶,𝐷) and designed an iterative method for its Hermitian minimum norm solutions. Peng et al. [9] presented an iterative method to obtain the least squares reflexive solutions of the matrix equations 𝐴₁𝑋𝐵₁=𝐶₁, 𝐴₂𝑋𝐵₂=𝐶₂. Cai et al. [10, 11] proposed iterative methods to solve the bisymmetric solutions of the matrix equations 𝐴₁𝑋𝐵₁=𝐶₁, 𝐴₂𝑋𝐵₂=𝐶₂.

In this paper, we propose a new iterative algorithm to solve the bisymmetric solution with the minimum Frobenius norm of the consistent matrix equations 𝐴₁𝑋𝐵₁=𝐶₁, 𝐴₂𝑋𝐵₂=𝐶₂, which is faster and more stable than Cai’s algorithm (Algorithm 2.1) in [10].

The rest of the paper is organized as follows. In Section 2, we propose an iterative algorithm to obtain the bisymmetric minimum Frobenius norm solution of (1) and present some basic properties of the algorithm. Some numerical examples are given in Section 3 to show the efficiency of the proposed iterative method.

2. A New Iterative Algorithm

Firstly, we give the following lemmas.

Lemma 1 (see [12]). There is a unique matrix 𝑃(𝑚,𝑛)∈𝐑^{𝑚𝑛×𝑚𝑛} such that vec(𝑋^𝑇)=𝑃(𝑚,𝑛)vec(𝑋) for all 𝑋∈𝐑^𝑚×𝑛. This matrix 𝑃(𝑚,𝑛) depends only on the dimensions 𝑚 and 𝑛. Moreover, 𝑃(𝑚,𝑛) is a permutation matrix and 𝑃(𝑛,𝑚)=𝑃(𝑚,𝑛)^𝑇=𝑃(𝑚,𝑛)⁻¹.

Lemma 2. If 𝑦₀,𝑦₁,𝑦₂,…∈𝐑^𝑚 are orthogonal to each other, then there exists a positive integer 𝑙̂≤𝑚 such that 𝑦̂𝑙=0.

Proof. If there exists a positive integer 𝑙̂≤𝑚−1 such that 𝑦̂𝑙=0, then Lemma 2 is proved.
Otherwise, we have 𝑦_𝑖≠0, 𝑖=0,1,2,…,𝑚−1, and 𝑦₀,𝑦₁,…,𝑦_𝑚−1 are orthogonal to each other in the 𝑚-dimension vector space of 𝐑^𝑚. So 𝑦₀,𝑦₁,…,𝑦_𝑚−1 form a set of orthogonal basis of 𝐑^𝑚.
Hence 𝑦_𝑚 can be expressed by the linear combination of 𝑦₀,𝑦₁,…,𝑦_𝑚−1. Denote 𝑦_𝑚=𝑎₀𝑦₀+𝑎₁𝑦₁+⋯+𝑎_𝑚−1𝑦_𝑚−1 (2) in which 𝑎_𝑖∈𝐑, 𝑖=0,1,2…,𝑚−1. Then 𝑦_𝑖,𝑦_𝑚=𝑎₀𝑦_𝑖,𝑦₀+𝑎₁𝑦_𝑖,𝑦₁+⋯+𝑎_𝑚−1𝑦_𝑖,𝑦_𝑚−1=𝑎_𝑖𝑦_𝑖,𝑦_𝑖+^𝑚−1_{𝑗=1𝑗≠𝑖}𝑎_𝑗𝑦_𝑖,𝑦_𝑗=𝑎_𝑖𝑦_𝑖,𝑦_𝑖,𝑖=0,1,2,…,𝑚−1. (3) From ⟨𝑦_𝑖,𝑦_𝑚⟩=0 and ⟨𝑦_𝑖,𝑦_𝑖⟩≠0, 𝑖=0,1,2,…,𝑚−1, we have 𝑎_𝑖=0,𝑖=0,1,2,…,𝑚−1; that is, 𝑦_𝑚=0. (4) This completes the proof.

Lemma 3. A matrix 𝑋∈𝐁𝐒𝐑^𝑛×𝑛 if and only if 𝑋=𝑋^𝑇=𝑆_𝑛𝑋𝑆_𝑛.

Lemma 4. If 𝑌∈𝐑^𝑛×𝑛, then 𝑌+𝑌^𝑇+𝑆_𝑛(𝑌+𝑌^𝑇)𝑆_𝑛∈𝐁𝐒𝐑^𝑛×𝑛.

Next, we review the algorithm proposed by Paige [13] for solving the following consistent problem: 𝑀𝑥=𝑓, (5) with given 𝑀∈𝐑^𝑠×𝑡,𝑓∈𝐑^𝑠.

Algorithm 5 (Paige algorithm). (i) Initialization 𝜏₀=1;𝜉₀=−1;𝜃₀𝑧=0;₀=0;𝑤₀=0;𝛽₁𝑢₁=𝑓;𝛼₁𝑣₁=𝑀^𝑇𝑢₁. (6)
(ii) Iteration. For 𝑖=1,2,…, until {𝑥_𝑖} convergence, do(a)𝜉_𝑖=−𝜉_𝑖−1𝛽_𝑖/𝛼_𝑖; 𝑧_𝑖=𝑧_𝑖−1+𝜉_𝑖𝑣_𝑖;(b)𝜃_𝑖=(𝜏_𝑖−1−𝛽_𝑖𝜃_𝑖−1)/𝛼_𝑖; 𝑤_𝑖=𝑤_𝑖−1+𝜃_𝑖𝑣_𝑖;(c)𝛽_𝑖+1𝑢_𝑖+1=𝑀𝑣_𝑖−𝛼_𝑖𝑢_𝑖; (d)𝜏_𝑖=−𝜏_𝑖−1𝛼_𝑖/𝛽_𝑖+1; (e)𝛼_𝑖+1𝑣_𝑖+1=𝑀^𝑇𝑢_𝑖+1−𝛽_𝑖+1𝑣_𝑖; (f)𝛾_𝑖=𝛽_𝑖+1𝜉_𝑖/(𝛽_𝑖+1𝜃_𝑖−𝜏_𝑖); (g)𝑥_𝑖=𝑧_𝑖−𝛾_𝑖𝑤_𝑖.

It is well known that if the consistent system of linear equations 𝑀𝑥=𝑓 has a solution 𝑥^∗∈𝑅(𝑀^𝑇), then 𝑥^∗ is the unique minimum Euclidean norm solution of 𝑀𝑥=𝑓. It is obvious that 𝑥_𝑖 generated by Algorithm 5 belongs to 𝑅(𝑀^𝑇) and this leads to the following result.

Theorem 6. The solution generated by Algorithm 5 is the minimum Euclidean norm solution of (5).
If 𝑢₁,𝑢₂,… and 𝑣₁,𝑣₂,… are generated by Algorithm 5, then 𝑢^𝑇_𝑖𝑢_𝑗=𝛿_𝑖𝑗, 𝑣^𝑇_𝑖𝑣_𝑗=𝛿_𝑖𝑗 (see details in [13]), in which 𝛿_𝑖𝑗=1,𝑖=𝑗,0,𝑖≠𝑗. (7) If we denote 𝑟_𝑖=𝑓−𝑀𝑥_𝑖, (8) where 𝑥_𝑖 is the approximation solution obtained by Algorithm 5 after the 𝑖th iteration, it follows that 𝑟_𝑖=−𝛽_𝑖+1𝜉_𝑖𝑢_𝑖+1 (see details in [13]). So we have 𝑟^𝑇_𝑖𝑟_𝑗=ℎ_𝑖𝑗𝛿_𝑖𝑗 (9) in which ℎ_𝑖𝑗=𝛽_𝑖+1𝛽_𝑗+1𝜉_𝑖𝜉_𝑗.

Now we derive our new algorithm, which is based on Paige algorithm.

Noting that 𝑋 is the bisymmetric solution of (1) if and only if 𝑋 is the bisymmetric solution of the following linear equations: 𝐴₁𝑋𝐵₁=𝐶₁,𝐵^𝑇₁𝑋𝐴^𝑇₁=𝐶^𝑇₁,𝐴₁𝑆_𝑛𝑋𝑆_𝑛𝐵₁=𝐶₁,𝐵^𝑇₁𝑆_𝑛𝑋𝑆_𝑛𝐴^𝑇₁=𝐶^𝑇₁,𝐴₂𝑋𝐵₂=𝐶₂,𝐵^𝑇₂𝑋𝐴^𝑇₂=𝐶^𝑇₂,𝐴₂𝑆_𝑛𝑋𝑆_𝑛𝐵₂=𝐶₂,𝐵^𝑇₂𝑆_𝑛𝑋𝑆_𝑛𝐴^𝑇₂=𝐶^𝑇₂. (10) Furthermore, suppose (10) is consistent; let 𝑌 be a solution of (10). If 𝑌 is a bisymmetric matrix, then 𝑌 is a bisymmetric solution of (1); otherwise we can obtain a bisymmetric solution of (10) by 𝑋=(𝑌+𝑌^𝑇+𝑆_𝑛(𝑌+𝑌^𝑇)𝑆_𝑛)/4.

The system of (10) can be transformed into (5) with coefficient matrix 𝑀 and vector 𝑓 as 𝐵𝑀=^𝑇1⊗𝐴₁𝐴₁⊗𝐵^𝑇₁𝐵^𝑇₁𝑆_𝑛⊗𝐴₁𝑆_𝑛𝐴₁𝑆_𝑛⊗𝐵^𝑇₁𝑆_𝑛𝐵^𝑇₂⊗𝐴₂𝐴₂⊗𝐵^𝑇₂𝐵^𝑇₂𝑆_𝑛⊗𝐴₂𝑆_𝑛𝐴₂𝑆_𝑛⊗𝐵^𝑇₂𝑆_𝑛𝐶,𝑓=vec₁vec𝐶^𝑇₁vec𝐶₁vec𝐶^𝑇₁vec𝐶₂vec𝐶^𝑇₂vec𝐶₂vec𝐶^𝑇₂. (11) Therefore, 𝛽₁𝑢₁=𝑓, 𝛼₁𝑣₁=𝑀^𝑇𝑢₁, 𝛽_𝑖+1𝑢_𝑖+1=𝑀𝑣_𝑖−𝛼_𝑖𝑢_𝑖, and 𝛼_𝑖+1𝑣_𝑖+1=𝑀^𝑇𝑢_𝑖+1−𝛽_𝑖+1𝑣_𝑖 can be written as 𝛽₁𝑢₁=vec𝐶₁vec𝐶^𝑇₁vec𝐶₁vec𝐶^𝑇₁vec𝐶₂vec𝐶^𝑇₂vec𝐶₂vec𝐶^𝑇₂,𝛼₁𝑣₁=𝐵^𝑇₁⊗𝐴₁𝐴₁⊗𝐵^𝑇₁𝐵^𝑇₁𝑆_𝑛⊗𝐴₁𝑆_𝑛𝐴₁𝑆_𝑛⊗𝐵^𝑇₁𝑆_𝑛𝐵^𝑇₂⊗𝐴₂𝐴₂⊗𝐵^𝑇₂𝐵^𝑇₂𝑆_𝑛⊗𝐴₂𝑆_𝑛𝐴₂𝑆_𝑛⊗𝐵^𝑇₂𝑆_𝑛^𝑇𝑢₁,𝛽_𝑖+1𝑢_𝑖+1=𝐵^𝑇₁⊗𝐴₁𝐴₁⊗𝐵^𝑇₁𝐵^𝑇₁𝑆_𝑛⊗𝐴₁𝑆_𝑛𝐴₁𝑆_𝑛⊗𝐵^𝑇₁𝑆_𝑛𝐵^𝑇₂⊗𝐴₂𝐴₂⊗𝐵^𝑇₂𝐵^𝑇₂𝑆_𝑛⊗𝐴₂𝑆_𝑛𝐴₂𝑆_𝑛⊗𝐵^𝑇₂𝑆_𝑛𝑣_𝑖−𝛼_𝑖𝑢_𝑖𝛼,𝑖=1,2,…,_𝑖+1𝑣_𝑖+1=𝐵^𝑇₁⊗𝐴₁𝐴₁⊗𝐵^𝑇₁𝐵^𝑇₁𝑆_𝑛⊗𝐴₁𝑆_𝑛𝐴₁𝑆_𝑛⊗𝐵^𝑇₁𝑆_𝑛𝐵^𝑇₂⊗𝐴₂𝐴₂⊗𝐵^𝑇₂𝐵^𝑇₂𝑆_𝑛⊗𝐴₂𝑆_𝑛𝐴₂𝑆_𝑛⊗𝐵^𝑇₂𝑆_𝑛^𝑇𝑢_𝑖+1−𝛽_𝑖+1𝑣_𝑖,𝑖=1,2,…. (12) From (12), we have 𝑢_𝑖=vec𝑈_𝑖1vec𝑈^𝑇_𝑖1vec𝑈_𝑖1vec𝑈^𝑇_𝑖1vec𝑈_𝑖2vec𝑈^𝑇_𝑖2vec𝑈_𝑖2vec𝑈^𝑇_𝑖2,𝑣_𝑖=vec𝑉_𝑖, (13) where 𝑈_𝑖1∈𝐑^𝑝1×𝑞1, 𝑈_𝑖2∈𝐑^𝑝2×𝑞2, 𝑉_𝑖∈𝐑^𝑛×𝑛, and 𝑉_𝑖 is a bisymmetric matrix.

And so, the vector form of 𝛽₁𝑢₁=𝑓, 𝛼₁𝑣₁=𝑀^𝑇𝑢₁, 𝛽_𝑖+1𝑢_𝑖+1=𝑀𝑣_𝑖−𝛼_𝑖𝑢_𝑖, and 𝛼_𝑖+1𝑣_𝑖+1=𝑀^𝑇𝑢_𝑖+1−𝛽_𝑖+1𝑣_𝑖 in Algorithm 5 can be rewritten as matrix form. Then we now propose the following matrix-form algorithm.

Algorithm 7. (i) Initialization  𝜏₀=1; 𝜉₀=−1;𝜃₀=0; 𝑍₀=0(∈𝐑^𝑛×𝑛); 𝑊₀=𝑍₀;  𝛽₁=2‖𝐶₁‖²+‖𝐶₂‖²; 𝑈_1𝑗=𝐶_𝑗/𝛽₁,𝑗=1,2;  𝑇₁=𝐴^𝑇₁𝑈₁₁𝐵^𝑇₁+𝐴^𝑇₂𝑈₁₂𝐵^𝑇₂; 𝑉₁=𝑇₁+𝑇^𝑇₁+𝑆_𝑛(𝑇₁+𝑇^𝑇₁)𝑆_𝑛; 𝛼₁=‖𝑉₁‖; 𝑉₁=𝑉₁/𝛼₁.
(ii) Iteration. For 𝑖=1,2,…, until {𝑋_𝑖} convergence, do(a)𝜉_𝑖=−𝜉_𝑖−1𝛽_𝑖/𝛼_𝑖;𝑍_𝑖=𝑍_𝑖−1+𝜉_𝑖𝑉_𝑖; (b)𝜃_𝑖=(𝜏_𝑖−1−𝛽_𝑖𝜃_𝑖−1)/𝛼_𝑖;𝑊_𝑖=𝑊_𝑖−1+𝜃_𝑖𝑉_𝑖; (c)𝑈_𝑖+1,𝑗=𝐴_𝑗𝑉_𝑖𝐵_𝑗−𝛼_𝑖𝑈_𝑖𝑗,𝑗=1,2;   𝛽_𝑖+1=2‖𝑈_𝑖+1,1‖²+‖𝑈_𝑖+1,2‖²;   𝑈_𝑖+1,𝑗=𝑈_𝑖+1,𝑗/𝛽_𝑖+1,𝑗=1,2;(d)𝜏_𝑖=−𝜏_𝑖−1𝛼_𝑖/𝛽_𝑖+1;(e)𝑇_𝑖+1=𝐴^𝑇₁𝑈_𝑖+1,1𝐵^𝑇₁+𝐴^𝑇₂𝑈_𝑖+1,2𝐵^𝑇₂;   𝑉_𝑖+1=𝑇_𝑖+1+𝑇^𝑇_𝑖+1+𝑆_𝑛(𝑇_𝑖+1+𝑇^𝑇_𝑖+1)𝑆_𝑛−𝛽_𝑖+1𝑉_𝑖;  𝛼_𝑖+1=‖𝑉_𝑖+1‖; 𝑉_𝑖+1=𝑉_𝑖+1/𝛼_𝑖+1;(f)𝛾_𝑖=𝛽_𝑖+1𝜉_𝑖/(𝛽_𝑖+1𝜃_𝑖−𝜏_𝑖);(g)𝑋_𝑖=𝑍_𝑖−𝛾_𝑖𝑊_𝑖.

Remark 8. The stopping criteria on Algorithm 7 can be used as 𝐶₁−𝐴₁𝑋_𝑖𝐵₁+𝐶₂−𝐴₂𝑋_𝑖𝐵₂𝜉≤𝜖,_𝑖≤𝜖or𝑋_𝑖−𝑋_𝑖−1≤𝜖, (14) where 𝜖 is a small tolerance.

Remark 9. As 𝑉_𝑖, 𝑍_𝑖, and 𝑊_𝑖 in Algorithm 7 are bisymmetric matrices, we can see that 𝑋_𝑖 obtained by Algorithm 7 are also bisymmetric matrices.

Some basic properties of Algorithm 7 are listed in the following theorems.

Theorem 10. The solution generated by Algorithm 7 is the bisymmetric minimum Frobenius norm solution of (1).

Theorem 11. The iteration of Algorithm 7 will be terminated in at most 𝑝₁𝑞₁+𝑝₂𝑞₂ steps in the absence of round-off errors.

Proof. By (8) and (11), we have by simple calculation that 𝑟_𝑖=vec𝑅_𝑖1vec𝑅^𝑇_𝑖1vec𝑅_𝑖1vec𝑅^𝑇_𝑖1vec𝑅_𝑖2vec𝑅^𝑇_𝑖2vec𝑅_𝑖2vec𝑅^𝑇_𝑖2, (15) in which 𝑅_𝑖1=𝐶₁−𝐴₁𝑋_𝑖𝐵₁, and 𝑅_𝑖2=𝐶₂−𝐴₂𝑋_𝑖𝐵₂, where 𝑋_𝑖 is the approximation solution obtained by Algorithm 7 after the 𝑖th iteration.
By Lemma 1, we have that vec𝑅^𝑇_𝑖1=𝑃𝑝₁,𝑞₁vec𝑅_𝑖1,vec𝑅^𝑇_𝑖2=𝑃𝑝₂,𝑞₂vec𝑅_𝑖2, (16) where 𝑃(𝑝₁,𝑞₁)∈𝐑^{𝑝1𝑞1×𝑝1𝑞1} and 𝑃(𝑝₂,𝑞₂)∈𝐑^{𝑝2𝑞2×𝑝2𝑞2} are permutation matrices. For simplicity, we denote 𝑃₁=𝑃(𝑝₁,𝑞₁), 𝑃₂=𝑃(𝑝₂,𝑞₂). Then 𝑃^𝑇₁𝑃₁=𝐼_𝑝1𝑞1,𝑃^𝑇₂𝑃₂=𝐼_𝑝2𝑞2.
Hence 𝑟^𝑇_𝑖𝑟_𝑗=vec𝑅_𝑖1vec𝑅^𝑇_𝑖1vec𝑅_𝑖1vec𝑅^𝑇_𝑖1vec𝑅_𝑖2vec𝑅^𝑇_𝑖2vec𝑅_𝑖2vec𝑅^𝑇_𝑖2^𝑇vec𝑅_𝑗1vec𝑅^𝑇_𝑗1vec𝑅_𝑗1vec𝑅^𝑇_𝑗1vec𝑅_𝑗2vec𝑅^𝑇_𝑗2vec𝑅_𝑗2vec𝑅^𝑇_𝑗2=vec𝑅_𝑖1𝑃₁vec𝑅_𝑖1vec𝑅_𝑖1𝑃₁vec𝑅_𝑖1vec𝑅_𝑖2𝑃₂vec𝑅_𝑖2vec𝑅_𝑖2𝑃₂vec𝑅_𝑖2^𝑇vec𝑅_𝑗1𝑃₁vec𝑅_𝑗1vec𝑅_𝑗1𝑃₁vec𝑅_𝑗1vec𝑅_𝑗2𝑃₂vec𝑅_𝑗2vec𝑅_𝑗2𝑃₂vec𝑅_𝑗2=2vec𝑅_𝑖1^𝑇vec𝑅_𝑗1+2vec𝑅_𝑖1^𝑇𝑃^𝑇₁𝑃₁vec𝑅_𝑗1+2vec𝑅_𝑖2^𝑇vec𝑅_𝑗2+2vec𝑅_𝑖2^𝑇𝑃^𝑇₂𝑃₂vec𝑅_𝑗2=4vec𝑅_𝑖1^𝑇vec𝑅_𝑗1+vec𝑅_𝑖2^𝑇vec𝑅_𝑗2𝑅=4vec_𝑖1vec𝑅_𝑖2^𝑇vec𝑅_𝑗1vec𝑅_𝑗2. (17)
If we let 𝑡_𝑖=^{vec(𝑅𝑖1)}_{vec(𝑅𝑖2)}∈𝐑^{𝑝1𝑞1+𝑝2𝑞2}, then we have by (9) that 𝑡₀,𝑡₁,𝑡₂,… are orthogonal to each other in 𝐑^{𝑝1𝑞1+𝑝2𝑞2}. By Lemma 2, there exists a positive integer 𝑙̂≤(𝑝₁𝑞₁+𝑝₂𝑞₂) such that 𝑡̂𝑙=0. Hence 𝑅̂𝑙¹𝑙̂=𝑅₂=0, (18) that is, the iteration of Algorithm 7 will be terminated in at most 𝑝₁𝑞₁+𝑝₂𝑞₂ steps in the absence of round-off errors.

3. Numerical Examples

In this section, we use some numerical examples to illustrate the efficiency of our algorithm. The computations are carried out at PC computer, with software MATLAB 7.0. The machine precision is around 10⁻¹⁶.

We stop the iteration when ‖𝑅_𝑖1‖+‖𝑅_𝑖2‖≤10⁻¹².

Example 12. Given matrices 𝐴₁, 𝐵₁, 𝐶₁, 𝐴₂, 𝐵₂, and 𝐶₂ as follows: 𝐴₁=,𝐵1−4−2−101−331−13−1−214−3−32−1−1−22514−1−34−14210−13−3−11−312−1₁=,𝐶−32−13−212−3−1−23−4−1101−110110−12123−1−253−30−33−30−1−101−2₁=,𝐴−19301119−3041−5547−855−4739−7477374−7780−3617−1936−17−219−30−11−1930−4155−478−5547−39₂=,𝐵3−2−11−40−10−31−3231−2−41−303103−13−2−3−11−60−2−430₂=,𝐶213−2−3−1−43123−10440−20−221−5−4−1−1−2−31₂=33107140−3317−34−17−1727−29−2−27−173417176078138−60 (19) then (1) is consistent, for one can easily verify that it has a bisymmetric solution: 𝑋=.1−1121−11−131111−1110−2−11121−21−21211−1−2011−111113−11−1121−11 (20) We choose the initial matrix 𝑋₀=0, then using Algorithm 7 and iterating 13 steps, we have the unique bisymmetric minimum Frobenius norm solution of (1) as follows:𝑋₁₃=,0.4755−0.68220.62741.45860.2774−1.2112−0.1053−0.68222.66280.40460.07161.01330.4001−1.21120.62740.4046−1.0215−2.2128−1.61761.01330.27741.45860.0716−2.2128−1.1548−2.21280.07161.45860.27741.0133−1.6176−2.2128−1.02150.40460.6274−1.21120.40011.01330.07160.40462.6628−0.6822−0.1053−1.21120.27741.45860.6274−0.68220.4755 (21)with ‖𝑅_13,1‖+‖𝑅_13,2‖=6.2303𝑒−013.

Figure 1 illustrates the performance of our algorithm and Cai’s algorithm [10]. From Figure 1, we see that our algorithm is faster than Cai’s algorithm.

Figure 1 

Convergence curves of log10(‖𝑅𝑖1‖+‖𝑅𝑖2‖).

Example 13. Let 𝐴₁=hilb(7),𝐵₁=pascal(7),𝐴₂=rand(7,7),𝐵₂=rand(7,7), (22) with hilb, pascal, and rand being functions in Matlab. And we let 𝐶₁=𝐴₁𝑋𝐵₁, 𝐶₂=𝐴₂𝑋𝐵₂, in which 𝑋 is defined in Example 12. Hence (1) is consistent.

We choose the initial matrix 𝑋₀=0; Figure 2 illustrates the performance of our algorithm and Cai’s algorithm [10]. From Figure 2, we see that our algorithm is faster and more stable than Cai’s algorithm.

Figure 2 

Convergence curves of log10(‖𝑅𝑖1‖+‖𝑅𝑖2‖).

Acknowledgments

This research was supported by the Natural Science Foundation of China (nos. 10901056, 11071079, and 11001167), the Shanghai Science and Technology Commission “Venus” Project (no. 11QA1402200), and the Natural Science Foundation of Zhejiang Province (no. Y6110043).