Asymptotic critical value algorithm - Lagrangian


An implementation of the algorithms given in above paper in the computer algebra system Maple. This implementation allows one to choose between several variants. Two use a Lagrangian-based geometric characterisation of the set of asymptotic critical values while the other uses a matrix a variables to construct a basis for a kernel of a certain Jacobian matrix. Then, one can combine these characterisations with elements of randomisation to efficiently compute a superset of said values.

AsympCritVals := proc(F, G, vars, out, char, slice, change, method)

### Input ###
# F - a list of polynomials that define a dominant polynomial mapping.
# G - a reduced regular sequence that defines a smooth algebraic set.
# vars - the variables of the polynomials of the input F.
# out - the variables of the output system.
# char - the characteristic of the field in which the computation is performed.
# slice - 1 (No randomisation in Lagrangian equations), 2 (Randomisation in Lagrangian equations), 3 (Alternate randomisation in Lagrangian equations)
# change - 0 (No generic linear change of variables), 1 (Generic linear change of variables)
# method - 1 (Lagrangian method), 2 (Kernel method)

### Output ###
# acv - a list a polynomials in the variables out whose simultaneous
# vanishing set contains the asymptotic critical values of the
# polynomial mapping defined by F restricted to the algebraic set defined by G.

local n, m, p, cs, FA, VA, varend, A, Asubs, J, i, j;

randomize():
n := nops(vars);
m := nops(V);
p := nops(F);
if change = 0 then

FA := F;
VA := V;
varend := m;

else

A := Matrix([seq([seq(rand(),i=1..n)],j=1..n)]);
Asubs := {seq(vars[i] = add(A(i,j)*vars[j],j=1..n),i=1..n)};
FA := subs(Asubs, F);
VA := subs(Asubs, V);
varend := 1;

fi;
J := VectorCalculus[Jacobian]([op(VA),op(FA)], vars);
if method = 2 then

return KernelAlgo(FA, VA, J, vars, char, n, m, p, varend, out, slice, change);

else

return LagrangeAlgo(FA, VA, J, vars, char, n, m, p, varend, out, slice, change);

fi;
end:

LagrangeAlgo(FA, VA, J, vars, char, n, m, p, varend, out, slice, change);

### Input ###
# FA - a list of polynomials that define a dominant polynomial mapping.
# VA - a reduced regular sequence that defines a smooth algebraic set.
# J - Jacobian of FA and VA
# vars - the variables of the polynomials of the input F.
# char - the characteristic of the field in which the computation is performed.
# n - number of vars
# m - number of VA
# p - number of FA
# varend - number of loops required
# out - the variables of the output system.
# slice - Level of randomisation
# change - Generic linear change of coordinates

### Output ###
# acv - a list a polynomials in the variables out whose simultaneous
# vanishing set contains the asymptotic critical values of the
# polynomial mapping defined by F restricted to the algebraic set defined by G.

local acv, lms, rlist, es, s, i, newvars, tau, varres, k, Vk, G, Gsat, cvals, l, j;

acv := [];
lms := [seq(cat(lm,i), i=1..m+p-1)];
rlist := [seq(rand(),i=1..n)];
es := [seq(cat(e,i), i=1..n)];
for s from 1 to varend do

newvars := subsop(s=NULL, vars);
tau := {vars[s]=1/vars[s], seq(newvars[j] = newvars[j]/vars[s], j=1..n-1)};
varres := [];
for k from 1 to p do

Vk := LinearAlgebra:-DeleteRow(J, m+k);
if slice = 1 then

G := numer(subs(tau, [seq(FA[i]-out[i],i=1..p),op(VA),seq(vars[s]*J[m+k][i] - add(lms[j]*Vk[j][i],j=1..m+p-1) - es[i],i=1..n)]));
Gsat := remove(has, Groebner:-Basis([op(G), l*vars[s]-1], lexdeg([l], [op(vars), op(es), op(lms), op(out)]), characteristic=char), {l});
cvals := remove(has, Groebner:-Basis([op(Gsat), op(es), vars[s]], lexdeg([op(vars),op(es),op(lms)],out), characteristic=char), {op(vars),op(es),op(lms)});

if slice = 2 then

G := numer(subs(tau, [seq(FA[i]-out[i],i=1..p),op(VA),seq(vars[s]*J[m+k][i] - rand()*es[1] - rand()*es[2]),i=1..n)]));
Gsat := remove(has, Groebner:-Basis([op(G), l*vars[s]-1], lexdeg([l], [op(vars), op(lms), es[1], es[2], op(out)]), characteristic=char), {l});
cvals := remove(has, Groebner:-Basis([op(Gsat), es[1], es[2], vars[s]], lexdeg([op(vars),es[1], es[2],op(lms)],out), characteristic=char), {op(vars),es[1], es[2], op(lms)});

else

Mbase := Matrix([[seq(vars[s]*J[m+k][i]-add(lms[j]*Vk[j][i]-rand()*es[1],j=1..m+p-1),i=1..n)],rlist]);
G := numer(subs(tau, [seq(FA[i]-out[i],i=1..p),op(VA),seq(linalg:-det(Mbase[[1,2],[1,i]]),i=2..n)]));
Gsat := remove(has, Groebner:-Basis([op(G), l*vars[s]-1], lexdeg([l], [op(vars), es[1], op(lms), op(out)]), characteristic=char), {l});
cvals := remove(has, Groebner:-Basis([op(Gsat), vars[s], es[1]], lexdeg([op(vars), es[1], op(lms)] out), characteristic=char), {op(vars), es[1], op(lms)});

fi;
varres := [op(varres), cvals];

od;
acv := [op(gcv), varres];

od;
return acv;
end:

KernelAlgo(FA, VA, J, vars, char, n, m, p, varend, out, slice, change);

### Input ###
# FA - a list of polynomials that define a dominant polynomial mapping.
# VA - a reduced regular sequence that defines a smooth algebraic set.
# J - Jacobian of FA and VA
# vars - the variables of the polynomials of the input F.
# char - the characteristic of the field in which the computation is performed.
# n - number of vars
# m - number of VA
# p - number of FA
# varend - number of loops required
# out - the variables of the output system.
# slice - Level of randomisation
# change - Generic linear change of coordinates

### Output ###
# acv - a list a polynomials in the variables out whose simultaneous
# vanishing set contains the asymptotic critical values of the
# polynomial mapping defined by F restricted to the algebraic set defined by G.

local acv, rlist, es, K, Um ks, s, newvars, tau, varres, cvals, k, j, i, Vk, Jac_M, Rank_M, Jac_eqs, Rank_eqs, Grad_eqs, G, Gsat, Cvals, l, Mbase;

acv := [];
rlist := [seq(rand(),i=1..n-m-p+1)];
es := [seq(cat(e,i), i=1..n-m-p+1)];
K := Matrix([seq([seq(cat(k,i,r,j),i=1..n-m-p+1)],j=1..n)]);
U := Matrix([seq([seq(rand(),i=1..n)],j=1..n-m-p+1)]);
ks := indets(K);
for s from 1 to varend do

newvars := subsop(s=NULL, vars);
tau := {vars[s]=1/vars[s], seq(newvars[j] = newvars[j]/vars[s], j=1..n-1)};
varres := [];
for k from 1 to p do

Vk := LinearAlgebra:-DeleteRow(J, m+k);
Jac_M := LinearAlgebra:-Multiply(Vk,K);
Jac_eqs := convert(Jac_M,list);
Rank_M := LinearAlgebra:-Multiply(U,K)-LinearAlgebra:-IdentityMatrix(n-m-p+1);
Rank_eqs := convert(Rank_M,list);
Grad_eqs := LinearAlgebra:-Multiply(J[m+k,1..n],K); if slice = 1 then

G := numer(subs(tau, [seq(FA[i]-out[i],i=1..p),op(VA),op(Rank_eqs),op(Jac_eqs),seq(vars[s]*Grad_eqs[i] - es[i],i=1..n)]));
Gsat := remove(has, Groebner:-Basis([op(G), l*vars[s]-1], lexdeg([l], [op(vars), op(es), op(ks), op(out)]), characteristic=char), {l});
cvals := remove(has, Groebner:-Basis([op(Gsat), op(es), vars[s]], lexdeg([op(vars),op(es),op(ks)],out), characteristic=char), {op(vars),op(es),op(ks)});

if slice = 2 then

G := numer(subs(tau, [seq(FA[i]-out[i],i=1..p),op(VA),op(Rank_eqs),op(Jac_eqs),seq(vars[s]*Grad_eqs[i] - rand()*es[1]-rand()*es[2],i=1..n)]));
Gsat := remove(has, Groebner:-Basis([op(G), l*vars[s]-1], lexdeg([l], [op(vars), es[1], es[2], op(ks), op(out)]), characteristic=char), {l});
cvals := remove(has, Groebner:-Basis([op(Gsat), es[1], es[2], vars[s]], lexdeg([op(vars),es[1], es[2],op(ks)],out), characteristic=char), {op(vars),es[1], es[2],op(ks)});

else

Mbase := Matrix([[seq(vars[s]*Grad_eqs[i]-rand()*es[1],i=1..n-m-p+1)],rlist]);
G := numer(subs(tau, [seq(FA[i] - out[i], i=1..p), op(VA), seq(linalg:-det(Mbase[[1,2],[1,i]]),i=2..n-m-p+1), op(Rank_eqs), op(Jac_eqs)]));
Gsat := remove(has, Groebner:-Basis([op(G), l*vars[s]-1], lexdeg([l], [op(vars), es[1], op(ks), op(out)]), characteristic=char), {l});
cvals := remove(has, Groebner:-Basis([op(Gsat), vars[s], es[1]], lexdeg([op(vars), es[1], op(ks)],out), characteristic=char), {op(vars), es[1], op(ks)});

fi;
varres := [op(varres), cvals];

od;
acv := [op(gcv), varres];

od;
return acv;
end: