A Homotopy Method for Computing All Isolated Solvents of the - - PowerPoint PPT Presentation

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

A Homotopy Method for Computing All Isolated Solvents of the Quadratic Matrix Equation AX 2 + BX + C = 0

Yongwen Hou Bo Yu

School of Mathematical Sciences, Dalian University of Technology

ASCM Beijing Oct. 26-28, 2012

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Outline

1

Introduction

2

Homotopy method for quadratic matrix equation

3

More results for special problems

4

Numerical results and conclusion

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Quadratic matrix equation

The unilateral quadratic matrix equation P(X) = AX 2 + BX + C = 0, (QME) is considered, whose coefficients A, B, and C ∈ Cn×n, and the matrix solution X is called a solvent. The corresponding quadratic eigenvalue problem is: P(λ)v = (λ2A + λB + C)v = 0, λ ∈ C and v = 0 ∈ Cn, (QEP) where λ and v are eigenvalue and eigenvector respectively.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Relations between solvents and eigenvalues.

Given a solvent X, satisfying P(X) = AX 2 + BX + C = 0, then, P(λ) can be divided by the linear term X − λI on the right: λ2A + λB + C = −(B + AX + λA)(X − λI), and thus, eigenvalues of P(λ) are those of X and those of a generalized eigenvalue problem (B + AX)v = −λAv.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Relations between solvents and eigenvalues.

For n eigenpairs {λi, vi}n

i=1 of P(λ), where the eigenvectors

v1, ..., vn are linearly independent, denoting V = [v1, ..., vn] and Λ = diag(λ1, ..., λn), then AV Λ2 + BV Λ + CV = 0 is satisfied. Multiplied by V −1 on the right, we have AV Λ2V −1 + BV ΛV −1 + C = 0, which indicates V ΛV −1 is a solvent of P(X).

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Solvents of quadratic matrix equation

P(X) can have no solvents, a finite positive number, or infinitely

• many. For example, the equation

X 2 − Jn(λ), where Jn(λ) (n > 1) is a Jordan block with the eigenvalue λ, has no solvents when λ = 0, and precisely two solvents when λ = 0. And the equation X 2 − In has infinite many solvents, including two isolated solvents and 2n − 2 families of solvents.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Existing algorithms

Newton’s method is attractive for its local quadratic convergence, if a good approximation Z0 of the desired solvent Z ∗ is given. It generates a sequence of matrices converging to Z ∗: Solve P′

Zk(Ek) = −P(Zk) for Ek

Update Zk+1 = Zk + Ek

• k = 0, 1, 2, ...

where P′

X(E) : Cn×n → Cn×n is the Fr´

echet derivative of P at X in the direction E along X.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Existing algorithms

Bernoulli’s method, used to find the dominant or minimal solvent (if there is one) of P(X), is generalized from the case of the quadratic scalar equations. Dominant solvent: (AZk + B)Zk−1 + C = 0, for k = 1, 2, ... lim

k→∞ Zk = Zdom.

Minimal solvent: (AZk−1 + B)Zk + C = 0, for k = 1, 2, ... lim

k→∞ Zk = Zmin.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Existing algorithms

Definition Suppose P(λ) has exactly 2n eigenvalues |λ1| ≥ |λ2| ≥ ... ≥ |λ2n|, and denote the set of eigenvalues of a matrix Z by λ(Z), a solvent Z1 of P(X) is a dominant solvent if λ(Z1) = {λ1, ..., λn}, and a solvent Z2 is a minimal solvent if λ(Z1) = {λn+1, ..., λ2n}.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Eigenvalue technique, used to construct all diagonalizable solvents. Theorem (Higham et al., 2001) Suppose Q(λ) = Mλ2 + Lλ + K has p distinct eigenvalues {λi}p

i=1, with n ≤ p ≤ 2n, and that the corresponding set of p

eigenvectors {λi}p

i=1 satisfies the Haar condition (that is, every

subset of n of them is linearly independent). Then there are at least p

n

• different solvents of Q(X), and exactly this many if

p = 2n, which are given by X = W diag(µi)W −1, W = [w1, ..., wn], where the eigenpaires {µi, wi}n

i=1 are chosen among the

eigenvpairs {λi, vi}p

i=1 of Q.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Existing algorithms

• N. J. Higham and H.-M. Kim, Numerical analysis of

quadratic matrix equation, IMA Journal of Numerical Analysis, 20 (2000), pp. 499-519.

• N. J. Higham and H.-M. Kim, Solving a quadratic matrix

equation by Newton’s method with exact line searches, SIAM Journal on Matrix Analysis and Applications, 23 (2001), pp. 303-316.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

The aim–finding all isolated solvents

We focus on locating all isolated solvents of P(X), of which there exists a neighborhood not containing other solvents. It is necessary to consider the element-wise form of P(X), which is a polynomial system of n2 equations and n2 variables, and denoted by: p(x) = 0. Here, we set x = vec(X) = (x11, ..., xn1, ..., x1n, ..., xnn)T, and p = vec(P) = (p11, ..., pn1, ..., p1n, ..., pnn)T.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Example 1:

The element-wise form of the equation X 2 + −1 −6 2 −9

• X +
• 12

−2 14

• = 0,

can be written as p(x) =     x2

11 + x12x21 − x11 − 6x21

x11x21 + x12x22 − x12 − 6x22 + 12 x21x11 + x22x21 + 2x11 − 9x21 − 2 x21x12 + x2

22 + 2x12 − 9x22 + 14

    = 0.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Homotopy method for solving polynomial systems

Given a system f (z) : Cr → Cr, we design a start system g(z) = 0, and construct a homotopy: h(z, t) = (1 − t)g(z) + tf (z), satisfying Triviality: The solutions of start system g(x) = 0 are known. Smoothness: The solution set of h(x, t) = 0 with t ∈ [0, 1) consists of a finite number of smooth paths emanating from isolated solutions of g(x) = 0, each parameterized by t. Accessibility: Every isolated solution of h(x, 1) = p(x) = 0 is reached by a path.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Homotopy method for solving polynomial systems

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

• Step1. Given a start solution z0 at t0 = 0, a stopping criteria

ǫ > 0, and an initial step h0.

• Step2. For i = 1, 2, ...

Predict: Predict solution (˜ zi, ti) such that: ti = ti−1 + hi−1 ˜ zi = zi−1 − hz(zi−1, ti−)−1ht(zi−1, ti−1)hi−1 Correct: Hold t constant, and correct ˜ zi = ˜ zi − h(˜ zi, ti)−1h(˜ zi, ti) until h(˜ zi, ti) < ǫ, then let zi = ˜ zi. Adjust: Cut the step length in half on failure of the corrector, and double it if several corrections at the current step size have been successful. Terminate At t = 1, refine z to high accuracy by Newton iteration.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Example 1.

There are different classes of homotopy methods, depending on different upper bounds estimation of the number of isolated

• solutions. To show the effect on Example 1, three of them have

been applied, in which we need to track 10 or more solution paths, while P(X) has exactly just 5 isolated solutions.

Homotopy Upper bound No.paths Classical homotopy B´ ezout number 16 m-Homogeneous homotopy m-Homogeneous B´ ezout number 16 Polyhedral homotopy BKK bound 10

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Two questions

1. Can we implement without transforming to a large polynomial system? 2. Can solution paths be reduced more in the path-following procedure?

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Homotopy method for quadratic matrix equation

We propose a homotopy method to solve P(X) = 0 in matrix form: H(X, t) = (1 − t)γQ(X) + tP(X), where Q(X) = MX 2 + LX + K = 0 is the start quadratic matrix equation.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Homotopy method for quadratic matrix equation

Theorem Let N(P) denote the number of isolated solvents of P(X) = AX 2 + BX + C, (1)N(P) is finite and it is the same, say N, for almost all A, B, and C ∈ Cn×n, where N = 2n n

• ;

(2)For all A, B, and C ∈ Cn×n, we have N(P) ≤ N.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Homotopy method for quadratic matrix equation

Theorem For the homotopy H(X, t) = γ(1 − t)Q(X) + tP(X), t ∈ [0, 1) and γ ∈ C where Q(X) = MX 2 + LX + K = 0 has exactly N mutually different isolated solvents, then, for almost all γ ∈ C, (1) The homotopy generates N smooth solution paths starting from isolated solvents of Q(X). (2) As t → 1, the limits of the solution paths include all isolated solvents of P(X).

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Homotopy method for quadratic matrix equation

Both of the theorems are based on the parameter continuation of isolated solutions of polynomial systems, that is, for a class of systems whose coefficients are continuous functions of some parameters, a continuous curve in the parameter space determines a continuous solution path accordingly.

• A. J. Sommese and C. W. Wampler, II, The numerical

solution of systems of polynomials, Arising in engineering and

• science. World Scientific Publishing Co. Pte. Ltd.,

Hackensack, NJ, 2005.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Construct the start equation

• Step1. Give 2n eigenpairs

{λi, vi}2n

i=1

in which λ′

is are different from each other, and the eigenvectors

satisfy the Haar condition.

• Step2. Construct a quadratic lambda polynomial

Q(λ) = Mλ2 + Lλ + K as well as Q(X), satisfying Q(λ)vi = λivi, for i = 1, ..., 2n.

• Step3. Compute all start solvents

X = W diag(µi)W −1, W = [w1, ..., wn], (1) where the eigenpaires (µi, wi)n

i=1 are chosen from among the

eigenvpairs {λi, vi}2n

i=1 of Q(λ).

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Construct the start equation

Our first approach is based on the method for solving a special kind of quadratic inverse eigenvalue problem, give by Chu and Xu in 2009:

• M. T. Chu and S.-F. Xu, Spectral decomposition of real

symmetric quadratic λ-matrices and its applications., Math. Comput., 78 (2009), pp. 293-313.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Corollary Choose eigenvalues λ1 > λ2 > ... > λ2n in R, and the columns of V := [In, Un] as the eigenvectors, where In and Un are n × n identity matrix and random orthogonal matrix respectively. If we let Λ := diag{λ1, ..., λ2n} and S = diag{In, −In}, then, the matrix coefficients M, L and K of Q(X) can be given by M = (V ΛSV T)−1, L = −MV Λ2SV TM, K = −MV Λ3SV TM + CM−1C.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Construct the start equation

Another method allows us to reconstruct a monic matrix polynomial Q(X) = X 2 + LX + K, along with its block companion matrix T of order n2 T =

• In

−K −L

• from the information of the associated eigenpairs.
• E. Pereira, On solvents of matrix polynomials, Applied

Numerical Mathematics, 47 (2003), pp. 197-208.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Corollary Randomly choose λ1, ..., λ2n ∈ C as the eigenvalues of Q(λ), and v1, ..., v2n ∈ Cn as the eigenvectors, then the 2n × 2n nonsingular matrix S =

• v1

· · · vn · · · v2n λ1v1 · · · λnvn · · · λ2nv2n

• is a similarity matrix of a block companion matrix associated with

Q(X), that is, T = Sdiag(λ1...λ2n)S−1. Thus, partitioning T as a 2 × 2 block matrix, we have M = In, L = −T22, K = −T12

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Fr´ echet derivative of the homotopy equation

Consider the homotopy equation H(X, t) = ((1−t)M +tA)X 2+((1−t)L+tB)X +((1−t)K +tC)X . = AtX 2 + BtX + CtX, and the expansion H(X + E, t) = H(X) + AtEX + (AtX + Bt)E + AtE 2, the Fr´ echet derivative of H at X in the direction E along X is H′

X(X, t; E) = AtEX + (AtX + Bt)E,

and the partial derivative of H to t is easy to get: H′

t(X, t) = (A − M)X 2 + (B − L)X + (C − K)X

.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Predictor-corrector scheme in matrix form

Step 1. Construct a start equation Q(X) = 0, along with a given start solvent X0, satisfying Q(X0) = 0; Step 2. Choose a tolerance ǫ > 0, and a step length h0; Step 3. Loop: for k = 0, 1, 2, ... Predict: Solve (Euler-Predictor) AtiEXi + (AtiXi + Bti)E = −H′

t(Xi, ti),

for E, and let X = Xi + hiE, and ti+1 = ti + hi

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Predictor-corrector scheme in matrix form

Correct: Solve (Newton-Corrector) Ati+1EX + (Ati+1X + Bti+1)E = −H(X, ti+1), for E and correct X in the subspace of t = ti+1 by X = X + E repeatedly until the condition H(X, ti+1) ≤ ǫ is satisfied, and let Xi+1 = X.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Predictor-corrector scheme in matrix form

Adjust: Double the step length if several successive corrections at the current step length have been successful, and cut the step length in half on failure of the corrector. Refine: Terminate and refine the solution at t = 1 using the Newton iteration. Next: Choose another start solvent and go to step 2.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Generalized sylvester equation

The generalized sylvester equation: AEX + (AX + B)E = F, involved in the algorithm can be solved by the generalized schur

• method. Apply generalized schur decomposition to A and AX + B,

and schur decomposition to X, we have W ∗AZ = T, W ∗(AX + B)Z = S, and U∗XU = R, and a substitution leads to T ˆ ER + S ˆ E = ˆ F, ˆ E = Z ∗EU and ˆ F = W ∗FU. Here, W , Z, and U are orthogonal matrices, while T, S, and S are upper diagonal matrices.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Generalized sylvester equation

Compute ˆ E column by column: (S + RkkT)ˆ Ek = ˆ Fk −

k−1

• i=1

RikT ˆ Ei, k = 1, 2, ..., n and the solution E = Z ˆ EU∗ is obtained. Hence, the predictor-corrector method respecting to the matrix form can be done.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Computing all isolated matrix square roots

Consider the equation X 2 − C = 0, which has at most 2n isolated solvents for general n × n complex matrix C, we present a homotopy map with a trivial start equation X 2 − D = 0, where D = diag(d1, · · · , dn) with nonzero diagonal elements different from each other, to find all isolated square roots of C if there exists.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Principal-preserving property

The principal square root, which is the unique square root for which every eigenvalue has nonnegative real part, can be located just by following one path, since the principal root goes to the principal along the homotopy path.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Overdamped problems

Definition The quadratic eigenvalue problem P(λ) = Aλ2 + Bλ + C is overdamped if A and B are symmetric positive definite, C is symmetric positive semidefinite and (xTBx)2 > 4(xTAx)(xTCx) for all x = 0.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Properties overdamped problems

The eigenvalues are real and nonpositive and there is a gap between the n largest and the n smallest; There are n linearly independent eigenvectors associated with the n largest eigenvalues and likewise for the n smallest ones; P(X) has at least two real solvents, having as their eigenvalues the n largest eigenvalues and for the n smallest

• nes respectively.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Order-preserving for overdamped problems

For overdamped problems, our algorithm guarantees the dominant solvent of the start equation goes to the dominant solvent of the target equation, and likewise for the minimal one, that means only

• ne solution curve needs to be tracked for computing a dominant
• r minimal solvent.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Numerical results

The algorithm is implemented in Matlab language, using Matlab’s qz to do the generalized Schur decomposition . All programmes are running on an ordinary computer, configured to Inter Core i5 2.76GHz processor, 8G RAM, Windows 7 Ultimate.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Example 2.

We present a comparison between our algorithm with brute-force method (here, applying the polyhedral homotopy method on the the polynomial system p(x) = 0 with n2 equation in n2 variables) for random matrix equations of the size n = 2, n = 3, and n = 4. n HOMOTOPY-QME HOM4PS2 PHC no.path time (s) no.path time (s) no.path time (s) 2 6 0.417 12 0.047 12 0.078 3 20 1.906 266 3.760 266 19.207 4 70 12.1486 26284 3037 26284 > 3600

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Example 3.

Consider the matrix       6 6 3 5 1 7 1 5 6 1 7 1 2 8 2 2 4 7 9 8 6 9 2 5 2       , all of its 32 isolated square roots can be found in 0.1s if the symmetry property is considered, since one of the roots X and −X is enough.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Recomputing Example 1

Again, consider the Example 1 X 2 + −1 −6 2 −9

• X +
• 12

−2 14

• = 0,

4

2

• = 6 solution paths are tracked in our homotopy, 5 of which

converged to its all isolated solvents: 1 2

• ,

1 2 3

• ,

3 1 2

• ,

1 4

• , and

4 2 2

• ,

and the other one diverged (converging to infinity).

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Conclusion

We have proposed an algorithm to locate all isolated solvents

• f the unilateral quadratic matrix equation, with the strategy

for choosing the start equation. In particular, it can be used to find all isolated square roots of a given matrix it there exists. We can locate a dominant or minimal solution for

• verdamped problems of great interest in applications by

tracing only one solution curve, and thus can be applied to solve larger problems.

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0

Introduction Homotopy method for quadratic matrix equation More results for special problems Numerical results and conclusion

Thanks for your attention!

Yongwen Hou, Bo Yu HM for Solving AX 2 + BX + C = 0