Normal Sets and Multiplication Matrices

Function: normset --- compute the normal set of a 0-dim Groebner basis

Function: mulmatrix --- compute a multiplication matrix from a normal set

References :

Robert M. Corless, ``Groebner Bases and Matrix Eigenproblems'',

Sigsam Bulletin (Communications in Computer Algebra) Volume 30,

number 4, Issue 118, December 1996, pp. 26--32.

David Cox, John Little, and Donal O'Shea, Ideals, Varieties and Algorithms, Springer-Verlag, 1992.

Calling Sequences:

> readlib(normset): # Get both routines

> ns := normset(gb, ordering):

> ns := normset(gb):

The default ordering is tdeg(op(indets(gb)), but you may specify any ordering

in the third argument. You must make sure that the ordering is the same as

that used to compute the Groebner basis gb. Only the case of zero-dimensional

varieties (finite-dimensional normal sets) is handled.

> mx := mulmatrix(x, ns, gb, ordering):

> my := mulmatrix(y, ns, gb, ordering):

Output:

The output of normset is an expression sequence with two elements. The first element,

> ns[1];

is a list containing the elements of the normal set of the Groebner basis. It will usually look like , where the are the terms of the normal set. The second element, , contains exactly the same information, but reversed in a table, for efficient access by the routine mulmatrix . You will not normally want to look at .

The output of mulmatrix is the multiplication matrix corresponding to the variable that was input. The eigenvalues of the multiplication matrices give the roots of the original polynomial system.

> restart:

Examples:

Example 1: A Lagrange multiplier problem.

First load in the routines.

> f := x^3 + 2*x*y*z - z^2:

> g := x^2 + y^2 + z^2 - 1:

> F := f + lambda*g;

> gF := convert(linalg[grad](F, [lambda, x, y, z]),set);

> with(Groebner):

> gb := gbasis(gF, tdeg(lambda,x,y,z)):

Now call normset to get the normal set of this Groebner basis.

> ns := normset(gb, tdeg(lambda,x,y,z)):

> ns[1];

We don't wish to look at these 12 by 12 matrices here, but we will

verify that the matrices commute.

> M[x] := mulmatrix(x, ns, gb, tdeg(lambda,x,y,z)):

> M[y] := mulmatrix(y, ns, gb, tdeg(lambda,x,y,z)):

> M[z] := mulmatrix(z, ns, gb, tdeg(lambda,x,y,z)):

> M[lambda] := mulmatrix(lambda, ns, gb, tdeg(lambda,x,y,z)):

> linalg[norm](evalm( M[x]&*M[y]-M[y]&*M[x] ));

Example 2: A geometric intersection problem.

> restart:

> f := 3*x^2*y + 9*x^2 +2*x*y+5*x +y-3;

> g := 2*x^3*y+6*x^3-2*x^2-x*y-3*x-y+3;

> h := x^3*y+3*x^3+x^2*y+2*x^2;

> with(Groebner):

> gb := gbasis({f,g,h}, tdeg(x,y)):

> ns := normset(gb, tdeg(x,y)):

> ns[1];

> M[x] := mulmatrix(x, ns, gb, tdeg(x,y));

> M[y] := mulmatrix(y, ns, gb, tdeg(x,y));

We now use the eigenvectors of a generic combination of these matrices to find all the possible roots. We find that this problem has only simple roots and so we may use eigenvectors to find out everything.

> alpha := rand(): beta := rand():

> M[generic] := evalm(alpha*M[x] + beta*M[y]):

> evs := linalg[eigenvects](M[generic]):

> evs := map(t -> (t[3][1]), [evs]);

> xev := linalg[vector](3,0):

> yev := linalg[vector](3,0):

> for i to 3 do

> xev[i] := evalm(M[x] &* evs[i] ):

> yev[i] := evalm(M[y] &* evs[i] ):

> od:

Now we use the eigenvectors to compute the roots.

> xvals := linalg[vector](3,0):

> yvals := linalg[vector](3,0):

> for i to 3 do

> xvals[i] := expand( radsimp(
xev[i][1]/evs[i][1], ratdenom )):

> yvals[i] := expand( radsimp(
yev[i][1]/evs[i][1], ratdenom )):

> od:

> answers := seq( {x=xvals[i],y=yvals[i]}, i=1..3);

So one root is x=0 and y=3. The others are almost as uncomplicated.

Example 3 . A problem with free parameters. This is taken from

the paper ``The Cheater's Homotopy'' by T. Y. Li, Tim Sauer, and J. A.

Yorke, from the SIAM J. Numerical Analysis, Volume 26, No. 5, pp. 1241-1251, October 1989. This problem cannot be done by a pure lexicographic ordered Groebner basis (the answer is just too complicated to be of any use even if we could calculate it in under two hours and 267 Mbytes).

But the approach here works in under a minute on a 486, with only a small

amount of memory.

> restart:

> LSY[1] := x^3*y^2+c1*x^3*y+y^2+c2*x+c3;

> LSY[2] := c4*x^4*y^2-x^2*y+y+c5;

> with(Groebner):

The following step does not succeed if you ask for a plex order basis.

> gb := gbasis({LSY[1],LSY[2]},tdeg(x,y)):

> ns := normset(gb, tdeg(x,y)):

> ns[1];

> nops(ns[1]);

> M[x] := mulmatrix(x, ns, gb, tdeg(x,y)):

> M[y] := mulmatrix(y, ns, gb, tdeg(x,y)):

At this point, one may insert numerical values for c1 through c5 and find eigenvalues of a generic (convex random) combination of these two matrices, cluster any multiple roots, and use the Schur vectors to find the roots of the system. We see that there are generically 10 roots. We take random values for the parameters below, as an example, ignoring the possibility of multiple

roots.

> c1 := rand()/10.^12: c2 := rand()/10.^12: c3 := rand()/10.^12: c4 := rand()/10.^12: c5 := rand()/10.^12;

> Digits := 20: V := 'V': readlib(forget)(evalf):

> xvals := evalf(Eigenvals(M[x],'V'));

Those eigenvectors stored in V are complex, and stored in a funny way, with real and imaginary parts separately.

We construct the correct matrix of eigenvectors from V as follows.

> Vc := array(1..10,1..10):

> for i to 10 do

> for j from 1 by 2 to 9 do

> Vc[i,j] := V[i,j] + I*V[i,j+1]:

> od:

> for j from 2 by 2 to 10 do

> Vc[i,j] := V[i,j-1] - I*V[i,j]:

> od:

> od:

> yvalmtx := map(fnormal,evalm( Vc^(-1) &* M[y] &* Vc)):

> yvals := [seq(yvalmtx[i,i],i=1..10)];

> residuals := [seq(subs(x=xvals[i],y=yvals[i],{LSY[1],LSY[2]}),i=1..10)];

>

Alternatively one can use tdeg bases to explore the bifurcation behaviour of these equations; one can quite easily determine parameter sets that force double, triple, or even quadruple roots.

See Also: solve , resultant , Groebner