A Note on the Inversion of Sylvester Matrices in Control Systems

Abstract

We give a sufficient condition (the solvability of two standard equations) of Sylvester matrix by using the displacement structure of the Sylvester matrix, and, according to the sufficient condition, we derive a new fast algorithm for the inversion of a Sylvester matrix, which can be denoted as a sum of products of two triangular Toeplitz matrices. The stability of the inversion formula for a Sylvester matrix is also considered. The Sylvester matrix is numerically forward stable if it is nonsingular and well conditioned.

1. Introduction

Let 𝑅[𝑥] be the space of polynomials over the real numbers. Given univariate polynomials 𝑓(𝑥),𝑔(𝑥)∈𝑅[𝑥], 𝑎₁≠0, where𝑓(𝑥)=𝑎₁𝑥^𝑛+𝑎₂𝑥^𝑛−1+⋯+𝑎_𝑛+1,𝑎₁≠0,𝑔(𝑥)=𝑏₁𝑥^𝑚+𝑏₂𝑥^𝑚−1+⋯+𝑏_𝑚+1,𝑏₁≠0.(1.1)

Let 𝑆 denote the Sylvester matrix of 𝑓(𝑥) and 𝑔(𝑥)∶𝑆=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝𝑎¹𝑎₂⋯⋯𝑎_𝑛+1𝑎₁𝑎₂⋯⋯𝑎_𝑛+1⋱⋱𝑎⋱₁𝑎₂⋯⋯𝑎_𝑛+1𝑏₁𝑏₂⋯𝑏_𝑚𝑏_𝑚+1𝑏₁𝑏₂⋯𝑏_𝑚𝑏_𝑚+1⋱⋱𝑏⋱⋱₁𝑏₂⋯𝑏_𝑚𝑏_𝑚+1⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠⎫⎪⎬⎪⎭⎫⎪⎬⎪⎭𝑚𝑛rrooww(1.2)

Sylvester matrix is applied in many science and technology fields. The solutions of Sylvester matrix equations and matrix inequations play an important role in the analysis and design of control systems. In determining the greatest common divisor of two polynomials, the Sylvester matrix plays a vital role, and the magnitude of the inverse of the Sylvester matrix is important in determining the distance to the closest polynomials which have a common root. Assuming that all principal submatrices of the matrix are nonsingular, in [1], Jing Yang et al. have given a fast algorithm for the inverse of Sylvester matrix by using the displacement structure of 𝑚+𝑛-order Sylvester matrix.

By using the displacement structure of the Sylvester matrix, in this paper, we give a sufficient condition (the solvability of two standard equations) of Toeplitz matrix, and, according to the sufficient condition, we derive a new fast algorithm for the inversion of a Sylvester matrix, which can be denoted as a sum of products of two triangular Toeplitz matrices. At last, the stability of the inversion formula for a Sylvester matrix is also considered. The Sylvester matrix is numerically forward stable if it is nonsingular and well conditioned.

In this paper, ‖⋅‖₂ always denotes the Euclidean or spectral norm and ‖⋅‖_𝐹 the Frobenius norm.

2. Preliminary Notes

In this section, we present a lemma that is important to our main results.

Lemma 2.1 (see [2, Section 2.4.8]). Let 𝐴,𝐵∈ℂ^𝑛,𝑛 and 𝛼∈ℂ. Then for any floating-point arithmetic with machine precision 𝜀, one has that 𝑓𝑙(𝑓𝐴𝑙+(𝛼𝐵𝐴))==𝐴𝛼𝐴++𝐵𝐸+,𝐸,‖𝐸‖‖𝐸_𝐹‖≤_𝐹𝜀|𝛼|‖𝐴‖_𝐹≤𝜀√𝑛|𝛼|‖𝐴‖₂,𝑓𝑙(𝐴𝐵)=𝐴𝐵+𝐸,‖𝐸‖_𝐹≤𝜀‖𝐴+𝐵‖_𝐹≤𝜀√𝑛‖𝐴+𝐵‖₂,≤𝜀𝑛‖𝐴‖_𝐹‖𝐵‖_𝐹.(2.1) As usual, one neglects the errors of O(𝜀²), O(𝜀̃²), and O(𝜀𝜀̃).

3. Sylvester Inversion Formula

In this section, we present our main results.

Theorem 3.1. Let matrix 𝑆 be a Sylvester matrix; then it satisfies the formula 𝐾𝑆−𝑆𝐾=𝑒_𝑚𝑓^𝑇−𝑒_𝑚+𝑛𝑔^𝑇,(3.1) where 𝐾=(0𝑒₁𝑒₂⋯𝑒_{𝑚+𝑛−1}) is a (𝑚+𝑛)×(𝑚+𝑛) shift matrix, 𝑓=𝑏₁⋯𝑏_𝑚𝑏_𝑚+1−𝑎₁⋯−𝑎_𝑛^𝑇,𝑔=0⋯0𝑏₁⋯𝑏_𝑚^𝑇.(3.2)

Proof. We have that 𝐾𝑆−𝑆𝐾=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝0𝑎¹⋯𝑎𝑎_𝑛+10⋯0₁𝑏⋯⋮𝑎_𝑛+1⋯0₁⋯𝑏𝑏_𝑚+10⋯0₁0⋯𝑏⋯⋮_𝑚+10⋯00⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠−⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝0𝑎¹0⋯⋯𝑎₁𝑎⋮_𝑛+10⋯0⋯𝑎_𝑛𝑏₁⋯𝑏_𝑚+10⋯0𝑏⋮₁⋯𝑏_𝑚=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝𝑏00⋯⋯⋮00₁⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠⋯𝑏_𝑚𝑏_𝑚+1−𝑎₁−𝑎₂⋯−𝑎_𝑛000⋯−𝑏₁⋯⋮00⋯−𝑏_𝑚𝑒_𝑚⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠𝑓^𝑇−𝑒_𝑚+𝑛𝑔^𝑇=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝𝑏00⋯𝑏_𝑚𝑏⋯⋯_𝑚+1⋮−𝑎₁−𝑎₂⋯−00₁𝑎_𝑛000⋯−𝑏₁⋯⋮00⋯−𝑏_𝑚⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠.(3.3) So 𝐾𝑆−𝑆𝐾=𝑒_𝑚𝑓^𝑇−𝑒_𝑚+𝑛𝑔^𝑇.

Theorem 3.2. Let matrix 𝑆 be a Sylvester matrix and 𝑥=(𝑥₁𝑥₂⋯𝑥_𝑚+𝑛)^𝑇, 𝑦=(𝑦₁𝑦₂⋯𝑦_𝑚+𝑛)^𝑇, 𝜇=(𝜇₁𝜇₂⋯𝜇_𝑚+𝑛)^T, and 𝑉=(𝑉₁𝑉₂⋯𝑉_𝑚+𝑛)^T the solutions of the systems of equations 𝑆𝑥=𝑒_𝑚, 𝑆𝑦=𝑒_𝑚+𝑛, 𝑆^𝑇𝜇=𝑓, and 𝑆^𝑇𝑉=𝑔, respectively, where 𝑒_𝑚 and 𝑒_𝑚+𝑛 are both (𝑚+𝑛)×1 vectors; then (a)𝑆 is invertible, and the column vector 𝑤_𝑗(𝑗=1,2,…,𝑚+𝑛) of 𝑆⁻¹ satisfies the recurrence relation𝑤𝑤_𝑚+𝑛=𝑦,_𝑖−1=𝐾𝑤_𝑖+𝜇_𝑖𝑥−𝑉_𝑖𝑦,𝑖=𝑚+𝑛,…,3,2,(3.4)(b)  𝑆⁻¹=𝑆₁𝑈₁+𝑆₂𝑈₂=⎛⎜⎜⎜⎜⎜⎝𝑦^𝑚+𝑛𝑦_{𝑚+𝑛−1}⋯𝑦₂𝑦₁𝑦_𝑚+𝑛⋱𝑦₃𝑦₂⋱𝑦⋱⋮_𝑚+𝑛𝑦_{𝑚+𝑛−1}𝑦_𝑚+𝑛⎞⎟⎟⎟⎟⎟⎠⎛⎜⎜⎜⎜⎜⎝−𝑉1_𝑚+𝑛−⋮𝑉₃⋮1⋱−𝑉₂−𝑉₄⋯1+⎛⎜⎜⎜⎜⎜⎝𝑥_𝑚+𝑛−𝑉₃⋯−𝑉_𝑚+𝑛1⎞⎟⎟⎟⎟⎟⎠𝑥_{𝑚+𝑛−1}⋯𝑥₂𝑥₁𝑥_𝑚+𝑛⋱𝑥₃𝑥₂⋱𝑥⋱⋮_𝑚+𝑛𝑥_{𝑚+𝑛−1}𝑥_𝑚+𝑛⎞⎟⎟⎟⎟⎟⎠⎛⎜⎜⎜⎜⎜⎝𝜇0_𝑚+𝑛𝜇⋮⋮0⋱₃𝜇𝜇₄⋯0₂𝜇₃⋯𝜇_𝑚+𝑛0⎞⎟⎟⎟⎟⎟⎠.(3.5)

Proof. By Theorem 3.1  𝑆𝑥=𝑒_𝑚 and 𝑆𝑦=𝑒_𝑚+𝑛, we have that 𝐾𝑆=𝑆𝐾+𝑒_𝑚𝑓^𝑇−𝑒_𝑚+𝑛𝑔^𝑇=𝑆𝐾+𝑆𝑥𝑓^𝑇−𝑆𝑦𝑔^𝑇=𝑆𝐾+𝑥𝑓^𝑇−𝑦𝑔^𝑇,(3.6) so 𝐾^𝑖𝑆==𝐾𝐾^𝑖−1_𝑖−2𝑆𝐾+𝑥𝑓^𝑇−𝑦𝑔^𝑇=𝐾_𝑖−2𝐾𝑆𝐾+𝑥𝑓^𝑇−𝑦𝑔^𝑇𝑆𝐾+𝑥𝑓^𝑇−𝑦𝑔^𝑇²=⋮𝑆𝐾+𝑥𝑓_𝑇−𝑦𝑔^𝑇^𝑖.(3.7) Hence, we have that 𝐾^𝑖𝑒_𝑚+𝑛=𝐾^𝑖𝑆𝑦=𝑆𝐾+𝑥𝑓^𝑇−𝑦𝑔^𝑇^𝑖𝑦,𝑖=0,1,…,𝑚+𝑛−1.(3.8) Let 𝑤_{𝑚+𝑛−𝑖}=𝐾+𝑥𝑓^𝑇−𝑦𝑔^𝑇^𝑖𝑦,𝑖=0,1,…,𝑚+𝑛−1.(3.9) It is easy to see that 𝐾^𝑖𝑒_𝑚+𝑛=𝑒_{𝑚+𝑛−𝑖}, and by (3.8) 𝑆𝑤_{𝑚+𝑛−𝑖}=𝐾^𝑖𝑒_𝑚+𝑛=𝑒_{𝑚+𝑛−𝑖},𝑖=0,1,…,𝑚+𝑛−1.(3.10) From 𝑆𝑦=𝑒_𝑚+𝑛, we have that 𝑤_𝑚+𝑛=𝑦. Let 𝑋=(𝑤₁𝑤₂⋯𝑤_𝑚+𝑛); then 𝑆𝑋=𝑆𝑤₁𝑤₂⋯𝑤_𝑚+𝑛=𝑆𝑤₁𝑆𝑤₂⋯𝑆𝑤_𝑚+𝑛=𝐼_𝑚+𝑛,(3.11) so the matrix 𝑆 is invertible and the inverse of 𝑆 is the matrix 𝑋.
From (3.1), we have that𝑆⁻¹(𝐾𝑆−𝑆𝐾)𝑆⁻¹=𝑆⁻¹𝑒_𝑚𝑓^𝑇−𝑒_𝑚+𝑛𝑔^𝑇𝑆⁻¹,(3.12) and thus 𝑆⁻¹𝐾−𝐾𝑆⁻¹=𝑥𝜇^𝑇−𝑦𝑉^𝑇.(3.13)
So𝑆⁻¹𝐾−𝐾𝑆⁻¹𝑒_𝑖=𝑥𝜇^𝑇−𝑦𝑉^𝑇𝑒_𝑖.(3.14)
Since  𝐾𝑒_𝑖=𝑒_𝑖−1, we have that𝑤_𝑖−1=𝐾𝑠_𝑖+𝜇_𝑖𝑥−𝑉_𝑖𝑦,𝑖=𝑚+𝑛,…,3,2,(3.15) and hence 𝑤𝑤_𝑚+𝑛=𝑦,_𝑖−1=𝐾𝑤_𝑖+𝜇_𝑖𝑥−𝑉_𝑖𝑦,𝑖=𝑚+𝑛,…,3,2.(3.16)
For (b), by  (3.4)𝑤𝑤_𝑚+𝑛=𝑦,_{𝑚+𝑛−1}𝑤=𝐾𝑦+𝜇_𝑚+𝑛𝑥−𝑉_𝑚+𝑛𝑦,_{𝑚+𝑛−2}𝑤=𝐾²𝑦+𝐾𝜇_𝑚+𝑛𝑥−⋯𝐾𝑉_𝑚+𝑛𝑦+𝜇_{𝑚+𝑛−1}𝑥−𝑉_{𝑚+𝑛−1}𝑦,₁=𝐾^{𝑚+𝑛−1}𝑦+𝐾^{𝑚+𝑛−2}𝜇_𝑚+𝑛𝑥−𝐾^{𝑚+𝑛−2}𝑉_𝑚+𝑛𝑦+⋯+𝜇₂𝑥−𝑉₂𝑦.(3.17) So 𝑆⁻¹==𝑤𝐾¹𝑤_{2𝑚+𝑛−1}⋯𝑤_𝑚+𝑛𝑦𝐾^{𝑚+𝑛−2}𝑦⋯𝐾𝑦𝑦⎛⎜⎜⎜⎜⎜⎝−𝑉1^𝑚+𝑛−⋮𝑉₃⋮1⋱−𝑉₂−𝑉₄⋯1+𝐾_{𝑚+𝑛−1}−𝑉₃⋯−𝑉_𝑚+𝑛1⎞⎟⎟⎟⎟⎟⎠𝑥𝐾^{𝑚+𝑛−2}𝑥⋯𝐾𝑥𝑥⎛⎜⎜⎜⎜⎜⎝𝜇0^𝑚+𝑛𝜇⋮⋮0⋱₃𝜇𝜇₄⋯0₂=⎛⎜⎜⎜⎜⎜⎝𝑦_𝑚+𝑛𝜇₃⋯𝜇_𝑚+𝑛0⎞⎟⎟⎟⎟⎟⎠𝑦_{𝑚+𝑛−1}⋯𝑦₂𝑦₁𝑦_𝑚+𝑛⋱𝑦₃𝑦₂⋱𝑦⋱⋮_𝑚+𝑛𝑦_{𝑚+𝑛−1}𝑦_𝑚+𝑛⎞⎟⎟⎟⎟⎟⎠⎛⎜⎜⎜⎜⎜⎝−𝑉1_𝑚+𝑛−⋮𝑉₃⋮1⋱−𝑉₂−𝑉₄⋯1+⎛⎜⎜⎜⎜⎜⎝𝑥_𝑚+𝑛−𝑉₃⋯−𝑉_𝑚+𝑛1⎞⎟⎟⎟⎟⎟⎠𝑥_{𝑚+𝑛−1}⋯𝑥₂𝑥₁𝑥_𝑚+𝑛⋱𝑥₃𝑥₂⋱𝑥⋱⋮_𝑚+𝑛𝑥_{𝑚+𝑛−1}𝑥_𝑚+𝑛⎞⎟⎟⎟⎟⎟⎠⎛⎜⎜⎜⎜⎜⎝𝜇0_𝑚+𝑛𝜇⋮⋮0⋱₃𝜇𝜇₄⋯0₂𝜇₃⋯𝜇_𝑚+𝑛0⎞⎟⎟⎟⎟⎟⎠.(3.18) This completes the proof.

4. Stability Analysis

In this section, we will show that the Sylvester inversion formula presented in this paper is evaluation forward stable.

If, for all well conditioned problems, the computed solution 𝑥̃ is close to the true solution 𝑥, in the sense that the relative error ‖𝑥−𝑥̃‖₂/‖𝑥‖₂ is small, then we call the algorithm forward stable (the author called this weakly in [3]). Round-off errors will occur in the matrix computation.

Theorem 4.1. Let matrix 𝑆 be a nonsingular Sylvester matrix and well conditioned; then the formula in Theorem 3.2 is forward stable.

Proof. Assume that we have computed the solutions 𝑥̃,𝑦,̃𝜇, and 𝑉 of 𝑆𝑥=𝑒_𝑚, 𝑆𝑦=𝑒_𝑚+𝑛, 𝑆^𝑇𝜇=𝑓, and 𝑆^𝑇𝑉=𝑔 in Theorem 3.2 which are perturbed by the normwise relative errors bounded by 𝜀̃, ‖𝑥̃‖₂≤‖𝑥‖₂(1+𝜀̃),‖𝑦‖̃₂≤‖𝑦‖₂(1+𝜀̃),‖‖𝜇‖‖₂≤‖𝜇‖₂(1+𝜀̃),‖‖𝑉‖‖₂≤‖𝑉‖₂(1+𝜀̃).(4.1)
Thus, we have that‖‖𝑈₁‖‖𝑆₁‖‖_𝐹≤√𝑚+𝑛‖𝑦‖₂,‖‖𝑆₂‖‖_𝐹≤√𝑚+𝑛‖𝑥‖₂,‖‖_𝐹≤√𝑚+𝑛1+‖𝑉‖²²,‖‖𝑈₂‖‖_𝐹≤√𝑚+𝑛‖𝜇‖₂.(4.2)
The inversion formula in Theorem 3.2, using the perturbed solutions 𝑥̃, 𝑦̃, 𝜇, and 𝑉, can be expressed as𝑆⁻¹=𝑓𝑙𝑆₁𝑈₁+𝑆₂𝑈₂=𝑓𝑙𝑆₁=𝑆₋₁+Δ𝑆₁𝑈₁+Δ𝑈₁+𝑆₂+Δ𝑆₂𝑈₂+Δ𝑈₂+Δ𝑆₁𝑈₁+𝑆₁Δ𝑈₁+Δ𝑆₂𝑈₂+𝑆₂Δ𝑈₂+𝐸+𝐹.(4.3) Here, and in the sequel, 𝐸 is the matrix containing the error which results from computing the matrix products and 𝐹 contains the error from subtracting the matrices. For the error matrices Δ𝑆₁, Δ𝑈₁, Δ𝑆₂, and Δ𝑈₂, we have that ‖‖Δ𝑈‖‖₁Δ𝑆₁‖‖_𝐹≤𝜀̃‖‖𝑆₁‖‖_𝐹≤𝜀̃√𝑚+𝑛‖𝑦‖₂,‖‖Δ𝑆₂‖‖_𝐹≤𝜀̃‖‖𝑆₂‖‖_𝐹≤𝜀̃√𝑚+𝑛‖𝑥‖₂,‖‖_𝐹≤𝜀̃‖‖𝑈₁‖‖_𝐹≤𝜀̃√𝑚+𝑛1+‖𝑉‖²²,‖‖Δ𝑈₂‖‖_𝐹≤𝜀̃‖‖𝑈₂‖‖_𝐹≤𝜀̃√𝑚+𝑛‖𝜇‖₂.(4.4)
By Lemma 2.1, we have the following bounds on 𝐸 and 𝐹:‖𝐸‖₂≤‖𝐸‖_𝐹≤≤𝜀𝜀((𝑚𝑚++𝑛𝑛))₂‖‖𝑆₁‖‖_𝐹‖‖𝑈₁‖‖_𝐹+‖‖𝑆₂‖‖_𝐹‖‖𝑈₂‖‖_𝐹‖𝑦‖²1+‖𝑉‖²²+‖𝑥‖₂‖𝜇‖₂≤𝜀(𝑚+𝑛)₂‖𝐹‖₂‖𝑦‖₂1+‖𝑉‖²²+‖𝑥‖₂‖𝜇‖₂≤√𝑚+𝑛𝜀‖‖𝑆⁻¹‖‖₂.(4.5)
Consequently, adding all these error bounds, by (4.3), we have that‖‖𝑆⁻¹−𝑆⁻¹‖‖₂≤(𝑚+𝑛)(2𝜀̃+(𝑚+𝑛)𝜀)‖𝑦‖₂1+‖𝑉‖²²+‖𝑥‖₂‖𝜇‖₂+√𝑚+𝑛𝜀‖‖𝑆⁻¹‖‖₂.(4.6)
From the equations 𝑆𝑥=𝑒_𝑚, 𝑆𝑦=𝑒_𝑚+𝑛𝑆^𝑇𝜇=𝑓, and 𝑆^𝑇𝑉=𝑔 in Theorem 3.2, we have that‖𝑦‖₂≤‖‖𝑆⁻¹‖‖₂,‖𝑉‖₂≤‖‖𝑆⁻¹‖‖₂‖𝑔‖₂,‖𝑥‖₂≤‖‖𝑆⁻¹‖‖₂,‖𝜇‖₂≤‖‖𝑆⁻¹‖‖₂‖𝑓‖₂.(4.7)
Thus, the relative error is‖‖𝑆⁻¹−𝑆⁻¹‖‖₂‖‖𝑆⁻¹‖‖₂≤(𝑚+𝑛)(2𝜀̃+(𝑚+𝑛)𝜀)‖𝑦‖²1+‖𝑉‖²²+‖𝑥‖₂‖𝜇‖₂₂‖‖𝑆⁻¹‖‖₂≤(𝑚+𝑛)(2𝜀̃+(𝑚+𝑛)𝜀)1+‖‖𝑆₋₁+𝜀√𝑚+𝑛‖‖₂‖𝑔‖₂+‖𝑓‖₂+𝜀√𝑚+𝑛.(4.8)
Since 𝑆 is well conditioned, ‖𝑆⁻¹‖₂ is finite. It is easy to see that ‖𝑔‖₂,‖𝑓‖₂ are finite. Therefore, the formula presented in Theorem 3.2 is forward stable.
This completes the proof.

5. Numerical Example

This section gives an example to illustrate our results. All the following tests are performed by MATLAB 7.0.

Example 5.1. Given 𝑓(𝑥)=𝑥+1,𝑔(𝑥)=𝑥²+𝑥+1, that is, 𝑎₁=1,𝑎₂=1,𝑏₁=1,𝑏₂=1,𝑏₃=1,𝑛=1,𝑚=2,𝑓=(1,1,0)^T,𝑔=(0,1,1)^T, and 𝑆=^110011111.
So𝑒𝐾𝑆−𝑆𝐾==⎛⎜⎜⎝⎛⎜⎜⎝000010100110010110⎞⎟⎟⎠⎞⎟⎟⎠⎛⎜⎜⎝−101⎛⎜⎜⎝000111101011⎞⎟⎟⎠111−⎞⎟⎟⎠⎛⎜⎜⎝=101⎛⎜⎜⎝111010011−01⎞⎟⎟⎠1⎛⎜⎜⎝000−001100⎞⎟⎟⎠010,⎞⎟⎟⎠(5.1)_𝑚𝑓^𝑇−𝑒_𝑚+𝑛𝑔^𝑇=⎛⎜⎜⎝010⎞⎟⎟⎠110−⎛⎜⎜⎝001⎞⎟⎟⎠011=⎛⎜⎜⎝010−011−001⎞⎟⎟⎠.(5.2) Therefore, 𝐾𝑆−𝑆𝐾=𝑒_𝑚𝑓^T−𝑒_𝑚+𝑛𝑔^T.
By the condition of Theorem 3.2, we can get𝑥=(−1,1,0)^T,𝑦=(1,−1,1)^T,𝜇=(1,0,0)^T,𝑉=(0,1,0)^T.(5.3) Obviously, 𝑆 is invertible and 𝑆⁻¹=^{0−1111−1−101}. And it is easy to see that 𝑤₃=𝑦, 𝑤₁===𝑆⎛⎜⎜⎝⎛⎜⎜⎝−000𝑤011₂100⎞⎟⎟⎠010⎞⎟⎟⎠⎛⎜⎜⎝−101⎞⎟⎟⎠+0⎛⎜⎜⎝−101⎞⎟⎟⎠−⎛⎜⎜⎝−111⎞⎟⎟⎠𝑤+𝜇₂𝑥−𝑉₂𝑦,₂===𝑆⎛⎜⎜⎝⎛⎜⎜⎝−000𝑤101₃100⎞⎟⎟⎠010⎞⎟⎟⎠⎛⎜⎜⎝−111⎞⎟⎟⎠+0⎛⎜⎜⎝−101⎞⎟⎟⎠−0⎛⎜⎜⎝−111⎞⎟⎟⎠𝑆+𝜇₃𝑥−𝑉₃𝑦,₋₁===𝑆⎛⎜⎜⎝⎛⎜⎜⎝−1₁011−11−011−11−111⎞⎟⎟⎠1⎛⎜⎜⎝⎞⎟⎟⎠−101101⎞⎟⎟⎠+⎛⎜⎜⎝010−101⎞⎟⎟⎠⎛⎜⎜⎝000000⎞⎟⎟⎠𝑈₁+𝑆₂𝑈₂.(5.4)

Acknowledgments

This work was supported financially by the Youth Foundation of Shandong Province, China (ZR2010EQ014), Shandong Province Natural Science Foundation (ZR2010AQ026), and National Natural Science Foundation of China (11026047).