Ideal-Theoretic Strategies for Asymptotic Approximation of Marginal Likelihood Integrals
Ideal-Theoretic Strategies for
Asymptotic Approximation of
Marginal Likelihood Integrals
Shaowei Lin
Research Paper
The research paper "Ideal-Theoretic Strategies for Asymptotic Approximation of Marginal Likelihood Integrals" can be found here.
Macaulay2 Libraries
MonomialMultiplierIdeals (Zach Teitler)
Zach Teitler also has a very fast implementation of finding T-distances of
Newton polyhedra. First, install Macaulay2 on your computer. Next, go to
Zach Teitler's website to download his MonomialMultiplierIdeals
package and follow his instructions for installing it.
Here are some examples for computing the T-distance, which should be
self explanatory.
R = QQ[x,y2,y3,y4,y5];
I = monomialIdeal(x*y2,y2^2*y3,y2^2*y4,y2^2*y5);
monomialThreshold(I,y2^4);
You should get that the RLCT = 3, and there is a 1x6 matrix following
that output which means the multiplicity equals 1. But one must be
careful: the matrix is a list of hyperplanes meeting at the
intersection point of the tau-diagonal with the Newton polyhedron, so
it is the codimension of this intersection, not the number of
hyperplanes, that gives us the multiplicity.
Here is another example where the multiplicity is not 1.
R = QQ[x,y];
I = monomialIdeal(x*y);
monomialThreshold(I,1_R);
Here, 1_R stands for the unit in the ring R and the output has a 2x3
matrix, which means that the multiplicity is 2. The RLCT is 1.
asymptotics.m2 (Shaowei Lin)
Download the library asymptotics.m2. The
RLCT (exact for monomial ideals, upper bound otherwise) can be
computed using the following code.
load "asymptotics.m2";
R = QQ[x,y2,y3,y4,y5];
I = ideal(x*y2,y2^2*y3,y2^2*y4,y2^2*y5);
RLCT(I,1)
RLCT(I,y2^4)
The RLCT(I,m) function computes the reciprocal of the tau-distance
of the Newton polyhedron of the ideal I, where tau is the exponents of the monomial m. Thus, if the input ideal is a monomial, the
output is precisely the RLCT; otherwise, it is only an upper bound. To
check if your input ideal can be generated by monomials, try the
following command to compute the Grobner basis. If the ideal is
monomial, then the Grobner basis will be a list of monomials.
gens gb I
The function newtonPolyhedron(I) returns the Newton polyhedron
of the input ideal. It uses the object type Polyhedron from
the Macaulay2 package Polyhedra. For polynomials f,
the function newtonPolyhedron(f) returns the Newton polyhedron
of f.
P = newtonPolyhedron I
The function tauDistance(P,t) computes the tau-distance of
the polyhedron P where tau is given by the list t of
integers. Here, the tau-distance is defined to be the smallest
rational number s such that the point s(t+1) lies in
P. The function returns a pair (l,Q) where l is the
tau-distance and Q is the tau-cone, i.e. the cone dual to the
face of the polyhedron P at the intersection of the
tau-diagonal with P. The dimension of Q is the
multiplicity of the tau-distance.
(l,Q) = tauDistance(P,{0,0,0,0,0})
dim Q
The function exponentsIdeal(I) returns the MonomialIdeal
generated by terms appearing in polynomials f in the input ideal.
exponentsIdeal(I)
The function faceIdeal(I,Q) returns
the face ideal of the ideal I at a face polyhedron Q, while
the function facePolynomial(f,Q) returns
the face polynomial of the polynomial f at the face
polyhedron Q.
Q = (faces(1,P))_0
f = x*y2+y2^2*y3
facePolynomial(f,Q)
faceIdeal(I,Q)
The function sosNondegenerate(I) checks if the ideal
I is sos-nondegenerate. It returns a pair (b,L)
where b is boolean and L is a list of faces of the
Newton polyhedron. If b==true, then the input ideal is
sos-nondegenerate. If b==false, then the function was not
able to determinine for sure if the ideal is sos-nondegenerate. Recall that an ideal I is
sos-nondegenerate if for all compact faces Q of the Newton
polyhedron P of I, the real variety of the face
ideal I_Q does not intersect the real torus (R^*)^d
where d is the number of variables in the base ring. Because
we do not have efficient algorithms for checking if a real variety is
empty, the implemented function first checks if the ideal is
monomial. If it is monomial, it returns that the ideal is
sos-nondegenerate. Otherwise, for each compact face Q of the
Newton polyhedron, it saturates the face ideal I_Q with the
product of the ring variables to remove components of the (complex)
variety that lies completely in the coordinate hyperplanes. If all the
resulting saturation ideals are the unit ideal, then the function
returns that the ideal is sos-nondegenerate. Otherwise, it returns a
list L of "bad faces" Q where the function could not
determinine if the real variety of I_Q intersects the
torus. This function can be slow, even for simple
examples, because saturating ideals is computationally intensive.
(b,L) = sosNondegenerate I
To find the face ideal corresponding to the 0-th bad face, we can use
the following code where we saturate out the coordinate hyperplanes to
find roots in the torus.
J = faceIdeal(I, L#0)
saturate (J, product flatten entries vars ring I)
The function nondegenerate(f) checks if the polynomial
f is nondegenerate. That means that it runs through all the
compact faces of the Newton polyhedron of f, computes the
face polynomial of f at that face, and checks if the singular
locus is empty in the torus. If all these singular loci are empty,
then f is nondegenerate. As with sosNondegenerate,
it returns a pair (b,L) where the boolean b states
if f is nondegenerate and the list L contains all
the "bad faces".
(b,L) = nondegenerate f
The function slocus(I,e) computes the ideal of the singular locus of the
ideal I with respect to the exceptional divisor defined by
e = 0 where e is a ring variable. This singular locus ideal is generated by generators
f of I, partial derivatives df/dx for
each non-exceptional variable x and e(df/de) for the
exceptional variable e.
The function blowupMap(R,c,L) returns a ring map that blows
up the origin of the ring R with respect to the
full-dimensional simplicial cone c. The list L is
a list of strings that will be used as names for the new variables in
the blowup ring.
The function strictTransform(I,L) saturates out the
coordinate hyperplanes corresponding to the variables in the list
L. The list L can be a list of variables or a list
of strings that are names of the variables.
The function jacobianDeterminant(p) returns the Jacobian
determinant of the ring map p.
The following code computes the Newton polyhedron of
the ideal I and finds a smooth refinement of its normal fan. It then
selects the 0-th cone in that smooth fan, and blows up the ideal using
the corresponding blowup map. It also computes the Jacobian
determinant of the map and the strict transform of the ideal.
P = newtonPolyhedron I;
T = normalToricVariety P;
S = makeSmooth T;
F = fan S;
C = cones (dim F,F);
U = {"e1","e2","e3","e4","e5","e6"};
p = blowupMap(ring I, C#0, U);
pI = p I
jacobianDeterminant(p)
strictTransform(pI, U)
The function simplifyRegularParameters(I) computes an ideal
whose RLCT is equal to that of the input ideal I. What it
does is that it evaluates at the origin the Jacobian matrix of the generators of the
ideal, and looks for ring variables that has a nonzero row in the
matrix. We can then project the variety of the ideal to the hyperplane
defined by this ring variable, and the image will have
the same RLCT modulo the dimension lost from the projection. In terms of ideals, this corresponds to eliminating the
above ring variable from the ideal. After running through all the
variables, we add the eliminated variables to the resulting ideal to
get an ideal that has the same RLCT as the original ideal. The
following example demonstrates this idea.
restart;
load "asymptotics.m2";
R = QQ[t,a1,a2,b1,b2,c1,c2,d1,d2];
A = matrix {{a1,a2,1-a1-a2}};
B = matrix {{b1,b2,1-b1-b2}};
C = matrix {{c1,c2,1-c1-c2}};
D = matrix {{d1,d2,1-d1-d2}};
P = t*(transpose A)*B + (1-t)*(transpose C)*D;
shift = map(R,R,{t+1/2,a1+1/3,a2+1/3,b1+1/3,b2+1/3,c1+1/3,c2+1/3,d1+1/3,d2+1/3});
eval = map(R,R,{1/2,1/3,1/3,1/3,1/3,1/3,1/3,1/3,1/3});
I = ideal (shift P - eval P)
I1 = simplifyRegularParameters I
The function removeUnitComponents(I) is another function that
can help simplify the ideal for which we are trying to compute the
RLCT at the origin. It removes components of the variety of the
ideal which are not passing through the origin. In terms of ideals, it
looks for associated primes of the input ideal which do not evaluate
to zero at the origin, and saturates the ideal
by these primes. Below, we continue from the above example. Note that
at the end, we get a monomial ideal for which we can compute the RLCT.
associatedPrimes I1
I2 = removeUnitComponents I1
associatedPrimes I2
RLCT(I2,1)
Alternatively, one can use the software Singular to remove units, such as
1-2*t in the previous example, from the
ideal in the local ring at the origin. Mora's algorithm for standard
bases, which is implemented in Singular, should do the
job. Here is the sample code in Singular for the previous
example.
ring R = 0,(t,a1,a2,b1,b2,c1,c2,d1,d2),ds;
ideal I = a1,a2,b1,b2,
2*t*c2*d2-c2*d2,2*t*c1*d2-c1*d2,
2*t*c2*d1-c2*d1,2*t*c1*d1-c1*d1;
groebner(I);
The code produces the following output.
_[1]=a1
_[2]=a2
_[3]=b1
_[4]=b2
_[5]=c1*d1
_[6]=c2*d1
_[7]=c1*d2
_[8]=c2*d2
Singular Library
You should have GFan
and Polymake installed on your
system.
The Singular library can be downloaded here: rlct.lib. You will also need the library polymake.lib.
Put both files in your work directory and type the command
LIB "rlct.lib";
when you start your Singular session.
The above library should work with Polymake version 2.9 and
above. If you are using Polymake version 2.3, use this instead: rlct0.lib. After downloading the file, remember to change its name to rlct.lib.
Functions
The command
list L = ndg(f);
determines if a
polynomial f is nondegenerate. The list L contains information
about the Newton polyhedron of f.
-
integer with value 1 if f is nondegenerate and 0 otherwise *
-
exponents of highest common monomial dividing f
-
f-vector of Newton polyhedron G of f
-
weight vectors for faces (non-vertices) of G
-
inequalities defining G
-
faces of G and the facets which meet them
-
polynomial Lf of lowest terms of f
-
dimension of Newton polytope of Lf
-
vertices of Newton polytope of Lf
-
if f is degenerate, a weight vector w giving a bad face
-
if f is degenerate, the initial form w.r.t. w
* The software currently checks nondegeneracy by checking the existence of complex roots, rather than the existence of real roots. This is because we do not know of any good libraries for real root computations. Thus, if L[1] = 1, then f is definitely nondegenerate, but if L[1] = 0, we are unsure of the nondegeneracy of f.
For an ideal I, use the command
list L = sosndg(I);
It returns a list L with 11 items as before, and it will contain information about the Newton polyhedron of I. Items 7 to 11 will be information about the sum of squares f of the given generators of I. If I is a monomial ideal, the software detects that and sets L[1] = 1. If I is not a monomial ideal, the software checks if the sum of squares f is nondegenerate. Right now, the software does not check if the ideal is sos-nondegenerate.
The command
rlct(L);
returns (1/l, t) where l is the distance of the
Newton polyhedron of f or of the ideal I and t its multiplicity. If
f is nondegenerate or if I is monomial or nondegenerate, this is precisely the real log canonical
threshold of f or I.
To find the distance in the direction of some non-negative vector T+1, use the command
rlct(L, T);
The default value of T is the zero vector. For instance, if the ring has 3 variables and we want to find the distance in the direction (1+1, 2+1, 3+1), we type
rlct(L, intvec(1,2,3));
Examples
LIB "rlct.lib";
ring RR = 0, (x,y,z), ds;
poly f = (x*y)^2+(y*z)^2+(x*z)^2;
list L = ndg(f);
if (L[1]==1)
{
print ("The polynomial is nondegenerate.");
print ("The RLCT is");
print (rlct(L));
}
else
{
print ("The polynomial may be degenerate.");
print ("One face polynomial that is singular in the torus is");
print (L[11]);
}
Suggestions and Bugs
If you find any bugs in the code or have any comments or suggestions, please do not hesitate to contact Shaowei Lin. His email address is
shaowei [underscore] lin [dot] sutd [dot] edu [dot] sg