Re: Symbolic solution of quadratic matrix equations

From: Gregg Kelly (greggk_at_ozemail.com.au)
Date: 11/30/04 

Date: 30 Nov 2004 12:28:44 -0800

In answer to Dave's question below, I believe I can shed a little more
light on this. Rewrite the matrix equation as

A.X^2 + B.X + C = 0 (1)

Let k and v be an eigenvalue and eigenvector of X so X.v = k.v
Then multiplying (1) by v gives

(A.k^2 + B.k + C).v = 0 (2)

so that

det(A.k^2 + B.k + C) = 0 (3)

This is then a polynomial equation of degree 2.n for the n eigenvalues
of X. Now we take n solutions k_i of (3) we can then solve (2) for
v_i.

Then X is given by the eigen decomposition formula

X = V.K.V^-1

where

K = diagonal matrix k_i
V = matrix of vectors v_i

Which n solutions of (3) should be chosen I don't know. Maybe they all
give valid solutions, maybe only some of them do ...

rusin@vesuvius.math.niu.edu (Dave Rusin) wrote in message news:<cog70j$l5u$1@news.math.niu.edu>...
> In article <codc2i$2p5$1@nntp.itservices.ubc.ca>,
> Robert Israel <israel@math.ubc.ca> wrote:
> >In article <b73a68dc.0411280909.64d5dbbe@posting.google.com>,
> >Ceylan <Ceylan.Yozgatligil@gmail.com> wrote:
> >
> >>I need to find a symbolic solution of a quadratic matrix equation
> >> X^2 * A + X * B + A' = 0
> >>where A' is the transpose of (kxk) nonsingular matrix A and B is a
> >>(kxk) symmetric matrix.
>
> >This is not going to be easy, except in the case
> >where all the matrices commute. You can treat it
> >as a system of k^2 quadratic polynomials in the entries
> >of X, and use try Groebner basis techniques, but I'd guess
> >it's likely to be rather complicated unless k is very small.
> >A 2 x 2 example should work pretty well, but even 3 x 3
> >might be very ugly.
>
> Of course Robert is correct; this is going to be a mess.
> But it looks like there IS some structure to the answer;
> I don't have an explanation for it but I offer my findings
> so someone else can explain the theory making this work.
>
> I tried some 2x2 and 3x3 examples, talking to Maple like this:
> X:=matrix(3,3,[z,y,x,w,v,u,t,s,r]);
> A:=matrix(3,3,[seq(rand() mod 100, i=1..9)]);
> B:=matrix(3,3,[seq(rand() mod 100, i=1..9)]);
> for i to 3 do for j to i-1 do B[i,j]:=B[j,i]:od:od:print(B);
> evalm(X&*X&*A + X&*B + transpose(A));
> eq:={seq(seq(%[i,j],j=1..3),i=1..3)}; lprint(%);
> and then to Magma like this:
> Q:=RationalField();
> P<z,y,x,w,v,u,t,s,r>:=PolynomialRing(Q,9);
> I:=ideal<P| ... nine quadratic polynomials not shown ... >;
> Groebner(I);
> Write("SeeMe",I);
>
> The 2x2 case is finished right away; the 3x3 case takes a minute or two.
>
> The 3x3 solution looks like this: each of z, y, x, ... can be expressed
> as a degree-19 polynomial in r with coefficients which are ratios of
> roughly 150-digit numbers. So there is one solution matrix X for
> each value of r. And there are 20 possible values of r, namely the
> roots of a polynomial of degree 20 with integer coefficients.
> (In the 2x2 case, replace "20" by "6".)
>
> The only surprise (to me) is that the polynomial always seems to factor:
> in the 2x2 case it's the product of a quadratic and a quartic with
> coefficients in the 7- and 15-digit range, respectively.
> In the 3x3 case there's a degree-8 factor (60-digit coefficients)
> and a degree-12 factor (90-digit coefficients). Moreover, these
> polynomials are "special": the quartic has a dihedral Galois group,
> and the two factors in the 3x3 case have galois groups of order just 48.
>
> So it appears there are two "types" of solutions to the original
> quadratic equation, and that at least for these low degrees,
> the different solutions can be obtained by some easy root extractions.
> I admit I don't understand why there is this much of a pattern to
> what should be a more random-looking polynomial.
>
> dave