source: git/ncpreim.lib @ 6b4d43

//////////////////////////////////////////////////////////////////////
version="$Id$";
category="Noncommutative";
info="
LIBRARY: ncpreim.lib    Non-commutative elimination and preimage computations
AUTHORS: Daniel Andres, daniel.andres@math.rwth-aachen.de

Support: DFG Graduiertenkolleg 1632 'Experimentelle und konstruktive Algebra'


OVERVIEW:
In G-algebras, elimination of variables is more involved than in the
commutative case.
One, not every subset of variables generates an algebra, which is again a
G-algebra.
Two, even if the subset of variables in question generates an admissible
subalgebra, there might be no admissible elimination ordering, i.e. an
elimination ordering which also satisfies the ordering condition for
G-algebras.

The difference between the procedure eliminateNC provided in this library and
the procedure eliminate (plural) from the kernel is that eliminateNC will
always find an admissible elimination if such one exists. Moreover, the use
of slimgb for performing Groebner basis computations is possible.

As an application of the theory of elimination, the procedure preimageNC is
provided, which computes the preimage of an ideal under a homomorphism
f: A -> B between G-algebras A and B. In contrast to the kernel procedure
preimage (plural), the assumption that A is commutative is not required.


References:
@* (BGL) J.L. Bueso, J. Gomez-Torrecillas, F.J. Lobillo:
         'Re-filtering and exactness of the Gelfand-Kirillov dimension',
         Bull. Sci. math. 125, 8, 689-715, 2001.
@* (GML) J.I. Garcia Garcia, J. Garcia Miranda, F.J. Lobillo:
         'Elimination orderings and localization in PBW algebras',
         Linear Algebra and its Applications 430(8-9), 2133-2148, 2009.
@* (Lev) Viktor Levandovskyy: 'Intersection of ideals with non-commutative
         subalgebras', ISSAC'06, 212-219, ACM, 2006.


PROCEDURES:
admissibleSub(v);          checks whether subalgebra is admissible
isUpperTriangular(M,k)  ;  checks whether matrix is (strictly) upper triangular
appendWeight2Ord(w);       appends weight to ordering
elimWeight(v);             computes elimination weight
eliminateNC(I,v,eng);      elimination in G-algebras
extendedTensor(A,I);       tensor product of rings with additional relations
preimageNC(A,f,J[,P,eng]); preimage of ideals under homomorphisms of G-algebras


KEYWORDS: 
preimage; elimination

SEE ALSO: eliminate_lib, preimage (plural)
";


LIB "elim.lib";    // for nselect
LIB "nctools.lib"; // for makeWeyl etc.
LIB "dmodapp.lib"; // for sortIntvec
LIB "ncalg.lib";   // for makeUgl
LIB "dmodloc.lib"; // for commRing


/*
CHANGELOG


*/


// -- Tools ------------------------------------------------


proc admissibleSub (intvec v)
"
USAGE:    admissibleSub(v);  v intvec
ASSUME:   The entries of v are in the range 1..nvars(basering).
RETURN:   int, 1 if the variables indexed by the entries of v form an 
          admissible subalgebra, 0 otherwise
EXAMPLE:  example admissibleSub; shows examples
"
{
  v = checkIntvec(v);
  int i,j;
  list RL = ringlist(basering);
  if (size(RL) == 4)
  {
    return(int(1));
  }
  matrix D = RL[6];
  ideal I;
  for (i=1; i<=size(v); i++)
  {
    for (j=i+1; j<=size(v); j++)
    {
      I[size(I)+1] = D[v[j],v[i]];
    }
  }
  ideal M = maxideal(1);
  ideal J = M[v];
  attrib(J,"isSB",1);
  M = NF(M,J);
  M = simplify(M,2); // get rid of double entries in v
  intvec opt = option(get);
  attrib(M,"isSB",1);
  option("redSB");
  J = NF(I,M);
  option(set,opt);
  for (i=1; i<=ncols(I); i++)
  {
    if (J[i]<>I[i])
    {
      return(int(0));
    }
  }
  return(int(1));
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(e,f,h),dp;
  matrix d[3][3];
  d[1,2] = -h; d[1,3] = 2*e; d[2,3] = -2*f;
  def A = nc_algebra(1,d);
  setring A; A; // A is U(sl_2)
  // the subalgebra generated by e,f is not admissible since [e,f]=h
  admissibleSub(1..2);
  // but the subalgebra generated by f,h is admissible since [f,h]=2f
  admissibleSub(2..3);
}


proc isUpperTriangular(matrix M, list #)
"
USAGE:    isUpperTriangular(M[,k]);  M a matrix, k an optional int
RETURN:   int, 1 if the given matrix is upper triangular, 
          0 otherwise.
NOTE:     If k<>0 is given, it is checked whether M is strictly upper
          triangular.
EXAMPLE:  example isUpperTriangular; shows examples
"
{
  int strict;
  if (size(#)>0)
  {
    if ((typeof(#[1])=="int") || (typeof(#[1])=="number"))
    {
      strict = (0<>int(#[1]));
    }
  }
  int m = Min(intvec(nrows(M),ncols(M)));
  int j;
  ideal I;
  for (j=1; j<=m; j++)
  {
    I = M[j..nrows(M),j];
    if (!strict)
    {
      I[1] = 0;
    }
    if (size(I)>0)
    {
      return(int(0));
    }
  }
  return(int(1));
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,x,dp;
  matrix M[2][3] = 
    0,1,2,
    0,0,3;
  isUpperTriangular(M);
  isUpperTriangular(M,1);
  M[2,2] = 4;
  isUpperTriangular(M);
  isUpperTriangular(M,1);
}


proc appendWeight2Ord (intvec w)
"
USAGE:    appendWeight2Ord(w);  w an intvec
RETURN:   ring, the basering equipped with the ordering (a(w),<), where < is
          the ordering of the basering.
EXAMPLE:  example appendWeight2Ord; shows examples
"
{
  list RL = ringlist(basering);
  RL[3] = insert(RL[3],list("a",w),0);
  def A = ring(RL);
  return(A);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(a,b,x,d),Dp;
  intvec w = 1,2,3,4;
  def r2 = appendWeight2Ord(w); // for a commutative ring
  r2;
  matrix D[4][4];
  D[1,2] = 3*a;  D[1,4] = 3*x^2;  D[2,3] = -x;
  D[2,4] = d;    D[3,4] = 1;
  def A = nc_algebra(1,D);
  setring A; A;
  w = 2,1,1,1;
  def B = appendWeight2Ord(w);  // for a non-commutative ring
  setring B; B;
}


static proc checkIntvec (intvec v)
"
USAGE:    checkIntvec(v);  v intvec
RETURN:   intvec consisting of entries of v in ascending order
NOTE:     Purpose of this proc: check if all entries of v are in the range
          1..nvars(basering).
"
{
  if (size(v)>1)
  {
    v = sortIntvec(v)[1];
  }
  int n = nvars(basering);
  if ( (v[1]<1) || v[size(v)]>n)
  {
    ERROR("Entries of intvec must be in the range 1.." + string(n));
  }
  return(v);
}



// -- Elimination ------------------------------------------


/*
// this is the same as Gweights@nctools.lib
// 
// proc orderingCondition (matrix D)
// "
// USAGE:    orderingCondition(D);  D a matrix
// ASSUME:   The matrix D is a strictly upper triangular square matrix.
// RETURN:   intvec, say w, such that the ordering (a(w),<), where < is
//           any global ordering, satisfies the ordering condition for
//           all G-algebras induced by D.
// NOTE:     If no such ordering exists, the zero intvec is returned.
// REMARK:   Reference: (BGL)
// EXAMPLE:  example orderingCondition; shows examples
// "
// {
//   if (ncols(D) <> nrows(D))
//   {
//     ERROR("Expected square matrix.");
//   }
//   if (isUpperTriangular(D,1)==0)
//   {
//     ERROR("Expected strictly upper triangular matrix.");
//   }
//   intvec v = 1..nvars(basering);
//   intvec w = orderingConditionEngine(D,v,0);
//   return(w);
// }
// example
// {
//   "EXAMPLE:"; echo = 2;
//   // (Lev): Example 2
//   ring r = 0,(a,b,x,d),dp;
//   matrix D[4][4];
//   D[1,2] = 3*a;  D[1,4] = 3*x^2;  D[2,3] = -x;
//   D[2,4] = d;    D[3,4] = 1;
//   // To create a G-algebra, the ordering condition implies 
//   // that x^2<a*d must hold (see D[1,4]), which is not fulfilled:
//   x^2 < a*d;
//   // Hence, we look for an appropriate weight vector
//   intwec w = orderingCondition(D); w;
//   // and use it accordingly.
//   ring r2 = 0,(a,b,x,d),(a(w),dp);
//   x^2 < a*d;
//   matrix D = imap(r,D);
//   def A = nc_algebra(1,D);
//   setring A; A;
// }
*/


proc elimWeight (intvec v)
"
USAGE:    elimWeight(v);  v an intvec
ASSUME:   The basering is a G-algebra.
@*        The entries of v are in the range 1..nvars(basering) and the
          corresponding variables generate an admissible subalgebra.
RETURN:   intvec, say w, such that the ordering (a(w),<), where < is
          any admissible global ordering, is an elimination ordering
          for the subalgebra generated by the variables indexed by the
          entries of the given intvec.
NOTE:     If no such ordering exists, the zero intvec is returned.
REMARK:   Reference: (BGL), (GML)
EXAMPLE:  example elimWeight; shows examples
"
{
  list RL = ringlist(basering);
  if (size(RL)==4)
  {
    ERROR("Expected non-commutative basering.");
  }
  matrix D = RL[6];
  intvec w = orderingConditionEngine(D,v,1);
  return(w);
}
example
{
  "EXAMPLE:"; echo = 2;
  // (Lev): Example 2
  ring r = 0,(a,b,x,d),Dp;
  matrix D[4][4];
  D[1,2] = 3*a;  D[1,4] = 3*x^2;  D[2,3] = -x;
  D[2,4] = d;    D[3,4] = 1;
  def A = nc_algebra(1,D);
  setring A; A;
  // Since d*a-a*d = 3*x^2, any admissible ordering has to satisfy
  // x^2 < a*d, while any elimination ordering for {x,d} additionally
  // has to fulfil a << x and a << d.
  // Hence neither a block ordering with weights
  // (1,1,1,1) nor a weighted ordering with weight (0,0,1,1) will do.
  intvec v = 3,4;
  elimWeight(v);
}


static proc orderingConditionEngine (matrix D, intvec v, int elimweight)
{
  // algorithm from (BGL) and (GML), respectively
  // solving an LPP via simplex
  int ppl = printlevel - voice + 1;
  def save = basering;
  int n = nvars(save);
  ideal EV = maxideal(1);
  EV = EV[v]; // also assumption check for v
  attrib(EV,"isSB",1);
  ideal NEV = maxideal(1);
  NEV = NF(NEV,EV);
  intmat V1[n-size(NEV)][n+1];
  if (elimweight)
  {
    intmat V2[size(NEV)][n+1];
  }
  int rowV1,rowV2;
  intmat M[1][n];
  intmat M2,oldM;
  int i,j,k;
  for (i=1; i<=n; i++)
  {
    if (elimweight)
    {
      if (NEV[i]<>0)
      {
        V2[rowV2+1,i+1] = 1; // xj == 0
        rowV2++;
      }
      else
      {
        V1[rowV1+1,1] = 1; // 1-xi <= 0
        V1[rowV1+1,i+1] = -1;
        rowV1++;
      }
    }
    else
    {
      V1[i,1] = 1; // 1-xi <= 0
      V1[i,i+1] = -1;
      rowV1++;
    }
    for (j=i+1; j<=n; j++)
    {
      if (deg(D[i,j])>0)
      {
        M2 = newtonDiag(D[i,j]);
        for (k=1; k<=nrows(M2); k++)
        {
          M2[k,i] = M2[k,i] - 1; // <beta,x> >= 0
          M2[k,j] = M2[k,j] - 1;
        }
        oldM = M;
        M = intmat(M,nrows(M)+nrows(M2),n);
        M = oldM,M2;
      }
    }
  }
  intvec eq = 0,(-1:n);
  ring r = 0,x,dp; // to avoid problems with pars or char>0
  module MM = module(transpose(matrix(M)));
  MM = simplify(MM,2+4);
  matrix A;
  if (MM[1]<>0)
  {
    if (elimweight)
    {
      MM = 0,transpose(MM);
    }
    else
    {
      MM = module(matrix(1:ncols(MM)))[1],transpose(MM);
    }
    A = transpose(concat(matrix(eq),transpose(-MM)));
  }
  else
  {
    A = transpose(eq);
  }
  A = transpose(concat(transpose(A),matrix(transpose(V1))));
  if (elimweight)
  {
    A = transpose(concat(transpose(A),matrix(transpose(V2))));
  }
  int m = nrows(A)-1;
  ring realr = (real,10),x,lp;
  matrix A = imap(r,A);
  dbprint(ppl,"// Calling simplex...");
  dbprint(ppl-1,"// with the matrix " + print(A));
  dbprint(ppl-1,"// and parameters " 
          + string(intvec(m,n,m-rowV1-rowV2,rowV1,rowV2)));
  list L = simplex(A,m,n,m-rowV1-rowV2,rowV1,rowV2);
  int se = L[2];
  if (se==-2)
  {
    ERROR("simplex yielded an error. Please inform the authors.");
  }
  intvec w = 0:n;
  if (se==0)
  {
    matrix S = L[1];
    intvec s = L[3];
    for (i=2; i<=nrows(S); i++)
    {
      if (s[i-1]<=n)
      {
        w[s[i-1]] = int(S[i,1]);
      }
    }
  }
  setring save;
  return(w);
}


proc eliminateNC (ideal I, intvec v, list #)
"
USAGE:    eliminateNC(I,v,eng);  I ideal, v intvec, eng optional int
RETURN:   ideal, I intersected with the subring defined by the variables not
          index by the entries of v
ASSUME:   The entries of v are in the range 1..nvars(basering) and the
          corresponding variables generate an admissible subalgebra.
REMARKS:  In order to determine the required elimination ordering, a linear
          programming problem is solved with the simplex algorithm.
@*        Reference: (GML)
@*        Unlike eliminate, this procedure will always find an elimination
          ordering, if such exists.
NOTE:     If eng<>0, @code{std} is used for Groebner basis computations,
          otherwise (and by default) @code{slimgb} is used.
@*        If printlevel=1, progress debug messages will be printed,
          if printlevel>=2, all the debug messages will be printed.
SEE ALSO: eliminate (plural)
EXAMPLE:  example eliminateNC; shows examples
"
{
  int ppl = printlevel - voice + 2;
  v = checkIntvec(v);
  if (!admissibleSub(v))
  {
    ERROR("Subalgebra is not admissible: no elimination is possible.");
  }
  dbprint(ppl,"// Subalgebra is admissible.");
  int eng;
  if (size(#)>0)
  {
    if (typeof(#[1])=="int" || typeof(#[1])=="number")
    {
      eng = int(#[1]);
    }
  }
  def save = basering;
  int n = nvars(save);
  dbprint(ppl,"// Computing elimination weight...");
  intvec w = elimWeight(v);
  if (w==(0:n))
  {
    ERROR("No elimination ordering exists.");
  }
  dbprint(ppl,"// ...done.");
  dbprint(ppl-1,"// Using elimination weight " + string(w) + ".");
  def r = appendWeight2Ord(w);
  setring r;
  ideal I = imap(save,I);
  dbprint(ppl,"// Computing Groebner basis with engine " + string(eng)+"...");
  I = engine(I,eng);
  dbprint(ppl,"// ...done.");
  dbprint(ppl-1,string(I));
  I = nselect(I,v);
  setring save;
  I = imap(r,I);
  return(I);
}
example
{
  "EXAMPLE:"; echo = 2;
  // (Lev): Example 2
  ring r = 0,(a,b,x,d),Dp;
  matrix D[4][4];
  D[1,2] = 3*a; D[1,4] = 3*x^2;
  D[2,3] = -x;  D[2,4] = d;     D[3,4] = 1;
  def A = nc_algebra(1,D);
  setring A; A;
  ideal I = a,x;
  // Since d*a-a*d = 3*x^2, any admissible ordering has to satisfy
  // x^2 < a*d, while any elimination ordering for {x,d} additionally
  // has to fulfil a << x and a << d.
  // Hence, the weight (0,0,1,1) is not an elimination weight for
  // (x,d) and the call eliminate(I,x*d); will produce an error.
  eliminateNC(I,3..4);
  // This call uses the elimination weight (0,0,1,2), which works.
}



// -- Preimages ------------------------------------------------

// TODO A or B commutative
proc extendedTensor(def A, ideal I)
"
USAGE:    extendedTensor(A,I);  A ring, I ideal
RETURN:   ring, A+B (where B denotes the basering) extended with non-
          commutative relations between the vars of A and B, which arise from
          the homomorphism A -> B induced by I in the usual sense, i.e. if the
          vars of A are named x(i) and the vars of B y(j), then putting
          q(i)(j) = leadcoef(y(j)*I[i])/leadcoef(I[i]*y(j)) and
          r(i)(j) = y(j)*I[i] - q(i)(j)*I[i]*y(j) yields the relation
          y(j)*x(i) = q(i)(j)*x(i)*y(j)+r(i)(j).
REMARK:   Reference: (Lev)
EXAMPLE:  example extendedTensor; shows examples
"
{
  def B = basering;
  setring A;
  int nA = nvars(A);
  string varA = "," + charstr(A) + "," + varstr(A) + ",";
  setring B;
  int nB = nvars(B);
  list RL = ringlist(B);
  list L = RL[2];
  string vB;
  int i,j;
  for (i=1; i<=nB; i++)
  {
    vB = "," + L[i] + ",";
    while (find(varA,vB)<>0)
    {
      vB[1] = "@";
      vB = "," + vB;
    }
    vB = vB[2..size(vB)-1];
    L[i] = vB;
  }
  RL[2] = L;
  def @B = ring(RL);
  kill L,RL;
  setring @B;
  ideal I = fetch(B,I);
  def E = A+@B;
  setring E;
  ideal I = imap(@B,I);
  matrix C = ringlist(E)[5];
  matrix D = ringlist(E)[6];
  poly p,q;
  for (i=1; i<=nA; i++)
  {
    for (j=nA+1; j<=nA+nB; j++)
    {
      // upper right block: new relations
      p = var(j)*I[i];
      q = I[i]*var(j);
      C[i,j] = leadcoef(p)/leadcoef(q);
      D[i,j] = p - C[i,j]*q;
    }
  }
  def @EE = commRing();
  setring @EE;
  matrix C = imap(E,C);
  matrix D = imap(E,D);
  def EE = nc_algebra(C,D);
  setring B;
  return(EE);
}
example
{
  "EXAMPLE:"; echo = 2;
  def A = makeWeyl(2);
  setring A; A;
  def B = makeUgl(2);
  setring B; B;
  ideal I = var(1)*var(3), var(1)*var(4), var(2)*var(3), var(2)*var(4);
  I;
  def C = extendedTensor(A,I);
  setring C; C;
}


proc preimageNC (list #)
"
USAGE:    preimageNC(A,f,J[,P,eng]);  A ring, f map or ideal, J ideal, 
                                      P optional string, eng optional int
ASSUME:   f defines a map from A to the basering.
RETURN:   nothing, instead exports an object 'preim' of type ideal to ring A,
          being the preimage of J under f.
NOTE:     If P is given and not equal to the empty string, the preimage is
          exported to A under the name specified by P.
          Otherwise (and by default), P is set to 'preim'.
@*        If eng<>0, @code{std} is used for Groebner basis computations,
          otherwise (and by default) @code{slimgb} is used.
@*        If printlevel=1, progress debug messages will be printed,
          if printlevel>=2, all the debug messages will be printed.
REMARK:   Reference: (Lev)
SEE ALSO: preimage (plural)
EXAMPLE:  example preimageNC; shows examples
"
{
  int ppl = printlevel - voice + 2;
  if (size(#) <3)
  {
    ERROR("Expected 3 arguments.")
  }
  def B = basering;
  if (typeof(#[1])<>"ring")
  {
    ERROR("First argument must be a ring.");
  }
  def A = #[1];
  setring A;
  ideal mm = maxideal(1);
  setring B;
  if (typeof(#[2])=="map")
  {
    map phi = #[2];
  }
  else
  {
    if (typeof(#[2])=="ideal")
    {
      map phi = A,#[2];
    }
    else
    {
      ERROR("Second argument must define a map from the specified ring to the basering.");
    }
  }
  if (typeof(#[3])<>"ideal")
  {
    ERROR("Third argument must be an ideal in the specified ring");
  }
  ideal J = #[3];
  string str = "preim";
  int eng;
  if (size(#)>3)
  {
    if (typeof(#[4])=="string")
    {
      if (#[4]<>"")
      {
        str = #[4];
      }
    }
    if (size(#)>4)
    {
      if (typeof(#[5])=="int")
      {
        eng = #[5];
      }
    }
  }
  setring B;
  ideal I = phi(mm);
  def E = extendedTensor(A,I);
  setring E;
  dbprint(ppl,"// Computing in ring");
  dbprint(ppl,E);
  int nA = nvars(A);
  int nB = nvars(B);
  ideal @B2E = maxideal(1);
  @B2E = @B2E[(nA+1)..(nA+nB)];
  map B2E = B,@B2E;
  ideal I = B2E(I);
  ideal Iphi;
  int i,j;
  for (i=1; i<=nA; i++)
  {
    Iphi[size(Iphi)+1] = var(i) - I[i];
  }
  dbprint(ppl,"// I_{phi} is  " + string(Iphi));
  ideal J = imap(B,J);
  J = J + Iphi;
  intvec v = (nA+1)..(nA+nB);
  dbprint(ppl,"// Starting elimination...");
  dbprint(ppl-1,string(J));
  J = eliminateNC(J,v,eng);
  dbprint(ppl,"// ...done.");
  dbprint(ppl-1,string(J));
  J = nselect(J,v);
  attrib(J,"isSB",1);
  setring A;
  dbprint(ppl,"// Writing output to specified ring under the name '" 
          + str + "'.");
  str = "ideal " + str + " = imap(E,J); export(" + str + ");";
  execute(str);
  setring B;
  return();
}
example
{
  "EXAMPLE:"; echo = 2;
  def A = makeUgl(3); setring A; A; // universal enveloping algebra of gl_3
  ring r3 = 0,(x,y,z,Dx,Dy,Dz),dp;
  def B = Weyl(); setring B; B;     // third Weyl algebra
  ideal ff = x*Dx,x*Dy,x*Dz,y*Dx,y*Dy,y*Dz,z*Dx,z*Dy,z*Dz;
  map f = A,ff;                     // f: A -> B, e(i,j) |-> x(i)D(j)
  ideal J = 0;
  preimageNC(A,f,J,"K");            // compute K := ker(f)
  setring A;
  K;
}


// -- Examples ---------------------------------------------

static proc ex1 ()
{
  ring r1 = 0,(a,b),dp;
  int t = 7;
  def St = nc_algebra(1,t*a);
  ring r2 = 0,(x,D),dp;
  def W = nc_algebra(1,1); // W is the first Weyl algebra
  setring W;
  map psit = St, x^t,x*D+t;
  int p = 3;
  ideal Ip = x^p, x*D+p;
  preimageNC(St,psit,Ip);
  setring St; preim;
}


static proc ex2 ()
{
  ring r1 = 0,(e,f,h),dp;
  matrix D1[3][3]; D1[1,2] = -h; D1[1,3] = 2*e; D1[2,3] = -2*f;
  def U = nc_algebra(1,D1); // D is U(sl_2)
  ring r2 = 0,(x,D),dp;
  def W = nc_algebra(1,1); // W is the first Weyl algebra
  setring W;
  ideal tau = x,-x*D^2,2*x*D;
  def E = extendedTensor(U,tau);
  setring E; E;
  elimWeight(4..5); 
  // zero, since there is no elimination ordering for x,D in E
}


static proc ex3 ()
{
  ring r1 = 0,(x,d,s),dp;
  matrix D1[3][3]; D1[1,2] = 1;
  def A = nc_algebra(1,D1);
  ring r2 = 0,(X,DX,T,DT),dp;
  matrix D2[4][4]; D2[1,2] = 1; D2[3,4] = 1;
  def B = nc_algebra(1,D2);
  setring B;
  map phi = A, X,DX,-DT*T;
  ideal J = T-X^2, DX+2*X*DT;
  preimageNC(A,phi,J);
  setring A;
  preim;
}