Code covered by the BSD License  

Highlights from
groebner

groebner

by

 

19 Jun 2009 (Updated )

manipulate and solve systems of multivariate polynomial equations by computing the groebner basis

gbasis=groebner(polyset,ord,varnames,tol)
function gbasis=groebner(polyset,ord,varnames,tol)
% GROEBNER - calculate the reduced Groebner basis of a set of polynomials
%  usage: gbasis = groebner(polyset,ord)
%         gbasis = groebner(polyset,ord,varnames)
%         gbasis = groebner(polyset,ord,varnames,tol)
%
% INPUTS: polyset, ord, varnames, tol
%  polyset is a cell array of multivariate polynomial coefficients
%    the cell elements can be all strings or all arrays.
%    If the polynomials are specified with arrays, each row of the array
%    represents a single term.  The first element is the coefficient and
%    the others are the degrees of each unknown in the term.
%    e.g. [-1,2,1;1,0,1] represents -x1^2*x2+x2
%    If the polynomials are specified as strings, they must use variable
%    names 'x1','x2',... (unless varnames is provided) and be valid inputs
%    for str2poly.m
%  ord is a string specifying the monomial ordering:
%    'lex': lexicographical, order by highest power of most significant
%      indeterminate, e.g. x1>x2^2, x1*x2>x1.
%    'grlex': graded lex, order by total degree then lex,
%      e.g. x1^2*x2 > x1*x2^2 > x1^2
%    'grevlex': graded reverse lex, order by total degree then by lowest
%      power of least significant indeterminate, e.g. x1*x2^2 > x1*x2*x3
%    (for all cases have x1>x2>...>xn)
%    note: revlex is not a well-ordering because 1>x1>x1^2>...
%  varnames (optional) is a cell array of variable names if polyset is a
%    cell array of strings and the variable names are not 'x1', 'x2', etc.
%    If specified, there must always be at least 2 variables.
%  tol (optional) is the zero tolerance for coefficients, otherwise
%    catastrophic cancellation may occur.  Default value is 0
%
% OUTPUTS: gbasis
%  gbasis is the reduced Groebner basis of the set of polynomials, in the
%    same format (strings or matrices) as polyset{1}.
%
% ALGORITHM:
%  A modified Buchberger's algorithm is used.  In Buchberger's algorithm,
%  at each step calculate Sij = (Lij/ai)*Pi - (Lij/aj)*Pj, for all pairs of
%  polynomials Pi, Pj, where ai and aj are the leading terms of Pi and Pj,
%  and Lij is the least common multiple of ai and aj; then reduce Sij with
%  respect to the polynomials and add it to the set (note that Sij always 
%  reduces to zero if ai and aj have no variables in common).  The
%  algorithm terminates when no new polynomials can be added to the set.
%  This function reduces the polynomials at each step in order to find the
%  reduced Groebner basis.
%
% EXAMPLE: simplify equations x^2+2xy^2=0, xy+2y^3=1
%  A1=[1,2,0;2,1,2]; % 1*x^2*y^0 + 2*x^1*y^2
%  A2=[-1,0,0;1,1,1;2,0,3]; % -1*x^0*y^0 + 1*x^1*y^1 + 2*x^0*y^3
%  y=groebner({A1,A2},'lex');
%  poly2str(y) % returns {'x2^3-0.5', 'x1'}
%
% EXAMPLE: same, but with equations specified as strings:
%  groebner({'x^2+2*x*y^2','x*y+2*y^3-1'},'lex',{'x','y'})
%  % returns {'y^3-0.5','x'}
%
% EXAMPLE: show that a set of equations is inconsistent
%  groebner({'x + y', 'x^2 - 1', 'y^2 - 2*x'}, 'lex', {'x', 'y'})
%  % returns {'1'}
%
%
% KNOWN ISSUES
%
%  0. GENERAL: while a significant effort has been made to ensure that the
%    program works as described, unidentified bugs may remain and the results
%    should be regarded with caution.  In particular, the suite of tests
%    performed was by no means comprehensive.
%
%  1. ACCURACY: the algorithm by default works in floating-point which
%    ultimately cannot perform exact arithmetic, so unexpected silent errors may
%    occur even when the coefficients are specified as integers.  A workaround
%    is to use other data types but this has its own issues (see below).
%
%  1a. EXAMPLE (roundoff - linearly independent polynomials)
%    eps % returns 2.2204e-16 (32-bit PC windows machine epsilon)
%    groebner({'x1+x2','x1+1.000000000000001*x2'},'lex') % ok, returns {'x2','x1'}
%    groebner({'x1+x2','x1+1.0000000000000001*x2'},'lex') % not ok, returns {'x2+x1'}
%
%  1b. EXAMPLE (roundoff - linearly dependent polynomials)
%    2/sqrt(2)-sqrt(2) % expect 0, returns -2.2204e-016
%    poly2str(groebner({[1,1,0;sqrt(2),0,1],[1/sqrt(2),1,0;1,0,1]},'lex')) % ok, returns {'1.41421*x2+x1'}
%    poly2str(groebner({[1,1,0;sqrt(2),0,1],[sqrt(2),1,0;2,0,1]},'lex')) % not ok, returns {'x2','x1'}
%    poly2str(groebner({[1,1,0;sqrt(2),0,1],[sqrt(2),1,0;2,0,1]},'lex',{},1e-14)) % ok, returns {'1.41421*x2+x1'}
%
%  1c. EXAMPLE (courtesy of Christophe Lauwerys - catastrophic cancellation & tolerance bug)
%    groebner({'t^3+x+y','t^2+0.5*x^2-x-z^2','t^2+y-z^2'},'lex',{'t','x','y','z'}) % ans{1} contains Inf
%
%  1d. EXAMPLE (from Mathematica documentation - catastrophic cancellation because of large intermediate coefficients)
%    % expect 3 polynomials returns, one of degree 8 in z; instead can get '1'
%    groebner({'x^2 + y^2 + z^2 - 1', 'x*y - z + 2', 'z^2 - 2*x + 3*y'}, 'lex', {'x', 'y','z'})
%
%  2. NON-DOUBLE INPUTS: use of non-double coefficients is not supported.  While
%    it may be tempting to introduce symbolic coefficients, it is not clear how
%    they should be consistently treated - for example 'A*x' could reduce to '0'
%    or 'x' or both, depending upon the assumptions on the coefficient A.
%
%  3. EFFICIENCY: The reduction algorithms aren't terribly efficient.  One of
%    the major bottlenecks at present is the calculation of the leading term,
%    which currently involves a call to sortrows each time.  It could be more
%    efficient to maintain an ordered list or even a priority queue (e.g.
%    implemented as a heap).  Also lattice reduction algorithms such as LLL
%    might speed things up.
%
%
% SEE ALSO:
%  poly2str, str2poly, polynsolve

% Author: Ben Petschel 19/6/2009
%
% Change history:
%  19/6/2009 - first release (ord='grlex' and ord='grevlex' not
%    implemented yet)
%  22/6/2009 - implemented ord='grlex' and ord='grevlex'
%    reduction algorithm now reduces lower-order terms as well
%  23/6/2009 - removed ord='revlex' because it is not a well-ordering
%  17/7/2009 - fixed bug in handling polynomials with >=3 variables
%  20/3/2010 - changed polynomial representation from multidimensional
%    arrays to 2d arrays (more memory efficient with large number of vars)
%  11/10/2010 - bugfix for handling error tolerance

if (numel(polyset)>0) && ischar(polyset{1}),
  charout = true;
  if nargin<3,
    polyset = str2poly(polyset);
  else
    polyset = str2poly(polyset,varnames);
  end;
else
  charout = false;
end;

if nargin<4,
  tol = 0;
end;

% determine the number of degrees
nd=max(cellfun(@(x)size(x,2)-1,polyset)); % error if = 0
if nd == 0
  error('number of variables must be >= 1');
end;

gbasis = fullreduce(polyset,ord,tol,nd);
oldgbasis = {};

while ~isequal(oldgbasis,gbasis),
  oldgbasis = gbasis;
  Qset = SPoly(gbasis,ord,tol,nd);
  gbasis = fullreduce([gbasis,Qset],ord,tol,nd);
end;

if charout,
  if nargin<3,
    gbasis = poly2str(gbasis);
  else
    gbasis = poly2str(gbasis,varnames);
  end;
end;

end % main function groebner(...)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Qset = SPoly(Pset,ord,tol,nd)
% calculates all reduced S-polynomials of Pset (Pset must be reduced)
%   Sij = (Lij/ai)*Pi - (Lij/aj)*Pj,
% for all pairs of polynomials Pi, Pj, where ai and aj are the leading
% terms of Pi and Pj, and Lij is the least common multiple of ai and aj
% then reduces Sij wrt Pset

Qset = {};
for i=1:numel(Pset)-1,
  for j=i+1:numel(Pset),
    [c1,d1] = leadterm(Pset{i},ord,tol,nd);
    [c2,d2] = leadterm(Pset{j},ord,tol,nd);
    if any(and(d1>0,d2>0)),
      % leading terms have a variable in common, so S-poly is nontrivial
      L = max(d1,d2); % LCM of leading terms
      S = addpolys(multiplyterm(Pset{i},L-d1),multconst(multiplyterm(Pset{j},L-d2),-1));
      S = reduceset(S,Pset,ord,tol,nd);
      %if any(S(:)~=0),
      if abs(leadterm(S,ord,tol,nd))>tol
        Qset{end+1}=S; % add S to the list
      end;
    end;
  end;
end;

end % helper function SPoly(...)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Qset = fullreduce(Pset,ord,tol,nd)
% reduces all polynomials wrt each other
Qset = Pset;
oldQset = {};

while ~isequal(oldQset,Qset),
  oldQset = Qset;
  for i=1:numel(Qset)-1,
    for j=i+1:numel(Qset),
      Qset{i}=reduce(Qset{i},Qset{j},ord,tol,nd);
      Qset{j}=reduce(Qset{j},Qset{i},ord,tol,nd);
    end;
  end;
end;

% keep only the nonzero results (to within zero tolerance tol)
Qset = Qset(cellfun(@(x)abs(leadterm(x,ord,tol,nd))>tol,Qset));

if numel(Qset)==1,
  % special case, a single equation can slip through without being reduced
  c = leadterm(Qset{1},ord,tol,nd);
  % coefficients occupy first column of the array
  Qset{1} = multconst(Qset{1},1/c); % know c~=0 because have already removed zero eqns
end;

end % helper function fullreduce(...)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q = reduceset(P,Pset,ord,tol,nd)
% reduces P wrt all polynomials in Pset
Q = P;
oldQ = [];

while ~isequal(oldQ,Q),
  oldQ = Q;
  for i=1:numel(Pset),
    Q=reduce(Q,Pset{i},ord,tol,nd);
  end;
end;

end % helper function reduceset(...)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q1=reduce(P1,P2,ord,tol,nd)
% reduces a polynomial P1 by subtracting multiples of P2 so that leading
% term coefficient is 1 and leading term of P2 does not divide that of P1

Q2=P2;
[c2,d2] = leadterm(Q2,ord,tol,nd);
if abs(c2)<=tol,
  % ignore terms less than tol
  keepgoing = false;
  Q1 = P1;
else
  keepgoing = true;
  Q2 = multconst(Q2,1/c2);
  Q1 = 0; % will successively add terms to Q
  % remove the leading terms, in order to reduce wrt lower terms
  remain = P1;
end;


while keepgoing,
  [c1,d1] = leadterm(remain,ord,tol,nd);
  if abs(c1)<=tol,
    % ignore terms less than tol
    keepgoing = false;
  else
    %reduce remain by Q2, if possible, or remove leading term
    if d1>=d2,
      % can reduce Q1 by Q2, so subtract a multiple of Q2 from Q1
      remain = addpolys(remain,multconst(multiplyterm(Q2,d1-d2),-c1));
    else
      % cannot reduce any further wrt leading term, but try lower terms
      lead = multiplyterm(c1,d1);
      remain = addpolys(remain,multconst(lead,-1));
      Q1 = addpolys(Q1,lead); % add irreducible terms to Q1
    end;
  end;
end;

% reduce leading coefficient, if possible
c1 = leadterm(Q1,ord,tol,nd);
if abs(c1)>tol,
  Q1 = multconst(Q1,1/c1);
end;

end % helper function reduce(...)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q=multconst(P,c)
% returns c*P where c is a constant

Q = P;
Q(:,1)=Q(:,1)*c;

end % helper function multconst(...)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q=addpolys(P1,P2)
% adds two polynomials

s1 = size(P1,2);
s2 = size(P2,2);

if s1>s2
  Q = [P1;P2,zeros(size(P2,1),s1-s2)];
else
  Q = [P1,zeros(size(P1,1),s2-s1);P2];
end;

Q = sortrows(Q,2:size(Q,2));
i=1;
while i<size(Q,1)
  % merge rows representing terms with the same exponents
  if Q(i,1)==0,
    Q(i,:)=[];
  elseif all(Q(i,2:end)==Q(i+1,2:end))
    Q(i,1)=Q(i,1)+Q(i+1,1);
    Q(i+1,:)=[];
    if Q(i,1)==0
      % remove terms that add to zero
      Q(i,:)=[];
    end;
  else
    i=i+1;
  end;
end;


end % helper function addpolys


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q=multiplyterm(P,d)
% multiplies a polynomial by a term (no coefficient)

% padding a matrix at the end is same as multiplying by x
d=d(:).';
sP = size(P,2);
sd = size(d,2)+1; % d doesn't include coefficient
if sP<sd,
  Q=[P,zeros(size(P,1),sd-sP)];
else
  Q=P;
  d=[(d(:)).',zeros(1,sP-sd)];
end;
Q=Q+repmat([0,d],size(Q,1),1);

end % helper function multiplyterm


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [coeff,deg]=leadterm(P,ord,tol,nd)
% returns the leading term of P
% deg is a row vector of length size(P,2)-1
% coeff is the leading coefficient, deg(i) is the degree of xi
% ignores terms with coefficients <= tol
% nd is the minimum number of degrees (in order to handle cases where deg=0)

d=size(P,2);

if all(abs(P(:,1))<=tol),
  % P=0
  coeff=0;
  deg=zeros(1,d-1);
else
  switch ord
    case 'lex'
      P=sortrows(P,2:d);

    case 'revlex'
      error('revlex is not a well-ordering');

    case 'grlex'
      P = [P,sum(P(:,2:d),2)];
      P = sortrows(P,[d+1,2:d]);
      P = P(:,1:d);

    case 'grevlex'
      P = [P,sum(P(:,2:d),2)];
      P = sortrows(P,[d+1,-(d:-1:2)]);
      P = P(:,1:d);
    
    otherwise
      error('ordering %s not recognized',ord);
      
  end; % switch ord

  ind = find(abs(P(:,1))>tol,1,'last');
  if isempty(ind),
    coeff = 0;
    deg = zeros(1,d-1);
  else
    coeff = P(ind,1);
    deg = P(ind,2:end);
  end;
  
end; % if all(P(:)==0)

if nd>d-1,
  deg = [deg,zeros(1,nd-d+1)];
end;

end % helper function leadterm(...)


Contact us