randomSquareFreeStep -- A step in a random walk with uniform distribution over all monomial ideals

Synopsis

Usage:

M = randomSquareFreeStep(I)

M = randomSquareFreeStep(I, AlexanderProbability => p)

M = randomSquareFreeStep(L)

M = randomSquareFreeStep(L, AlexanderProbability => p)
Inputs:
- I, an ideal, square-free monomial Ideal or MonomialIdeal
- L, a list, {I,Igens,ISocgens} where I is a square-free MonomialIdeal, Igens is a List of its minimal generators, ISocgens is a List of the minimal generators of the socle mod I.
Optional inputs:
- AlexanderProbability => ..., default value .05
Outputs:
- M, a list, {J,Jgens,JSocgens} where J is a square-free MonomialIdeal, Jgens is a List of its minimal generators, JSocgens is a List of the minimal generators of the socle mod J.

Description

With probability p the routine takes the Alexander dual of I; the default value of p is .05, and it can be set with the option AlexanderProbility.

Otherwise uses the Metropolis algorithm to produce a random walk on the space of square-free ideals. Note that there are a LOT of square-free ideals; these are the Dedekind numbers, and the sequence (with 1,2,3,4,5,6,7,8 variables) begins 3,6,20,168,7581, 7828354, 2414682040998, 56130437228687557907788. (see the Online Encyclopedia of Integer Sequences for more information). Given I in a polynomial ring S, we make a list ISocgens of the square-free minimal monomial generators of the socle of S/(squares+I) and a list of minimal generators Igens of I. A candidate "next" ideal J is formed as follows: We choose randomly from the union of these lists; if a socle element is chosen, it's added to I; if a minimal generator is chosen, it's replaced by the square-free part of the maximal ideal times it. the chance of making the given move is then 1/(#ISocgens+#Igens), and the chance of making the move back would be the similar quantity for J, so we make the move or not depending on whether random RR < (nJ+nSocJ)/(nI+nSocI) or not; here random RR is a random number in [0,1].

i1 : setRandomSeed(currentTime())

o1 = 1711041345

i2 : S=ZZ/2[vars(0..3)]

o2 = S

o2 : PolynomialRing

i3 : J = monomialIdeal"ab,ad, bcd"

o3 = monomialIdeal (a*b, a*d, b*c*d)

o3 : MonomialIdeal of S

i4 : randomSquareFreeStep J

o4 = {monomialIdeal (a*b, a*c*d, b*c*d), {a*b, a*c*d, b*c*d}, {c*d, b*d, a*d,
     ------------------------------------------------------------------------
     b*c, a*c}}

o4 : List

With 4 variables and 168 possible monomial ideals, a run of 5000 takes less than 6 seconds on a reasonably fast machine. With 10 variables a run of 1000 takes about 2 seconds.

i5 : setRandomSeed(1)

o5 = 1

i6 : rsfs = randomSquareFreeStep

o6 = randomSquareFreeStep

o6 : MethodFunctionWithOptions

i7 : J = monomialIdeal 0_S

o7 = monomialIdeal ()

o7 : MonomialIdeal of S

i8 : time T=tally for t from 1 to 5000 list first (J=rsfs(J,AlexanderProbability => .01));
     -- used 2.32331 seconds

i9 : #T

o9 = 168

i10 : T

o10 = Tally{monomialIdeal () => 45                            }
            monomialIdeal (a*b*c, a*b*d) => 33
            monomialIdeal (a*b*c, a*b*d, a*c*d) => 22
            monomialIdeal (a*b*c, a*b*d, a*c*d, b*c*d) => 26
            monomialIdeal (a*b*c, a*b*d, b*c*d) => 37
            monomialIdeal (a*b*c, a*b*d, c*d) => 31
            monomialIdeal (a*b*c, a*c*d) => 25
            monomialIdeal (a*b*c, a*c*d, b*c*d) => 25
            monomialIdeal (a*b*c, a*d) => 31
            monomialIdeal (a*b*c, a*d, b*c*d) => 33
            monomialIdeal (a*b*c, a*d, b*d) => 36
            monomialIdeal (a*b*c, a*d, b*d, c*d) => 35
            monomialIdeal (a*b*c, a*d, c*d) => 47
            monomialIdeal (a*b*c, b*c*d) => 30
            monomialIdeal (a*b*c, b*d) => 39
            monomialIdeal (a*b*c, b*d, a*c*d) => 36
            monomialIdeal (a*b*c, b*d, c*d) => 39
            monomialIdeal (a*b*c, c*d) => 24
            monomialIdeal (a*b*c, d) => 20
            monomialIdeal (a*b*d, a*c*d) => 35
            monomialIdeal (a*b*d, a*c*d, b*c*d) => 32
            monomialIdeal (a*b*d, b*c*d) => 33
            monomialIdeal (a*b*d, c*d) => 45
            monomialIdeal (a*b, a*c) => 37
            monomialIdeal (a*b, a*c*d) => 42
            monomialIdeal (a*b, a*c*d, b*c*d) => 31
            monomialIdeal (a*b, a*c, a*d) => 29
            monomialIdeal (a*b, a*c, a*d, b*c*d) => 38
            monomialIdeal (a*b, a*c, a*d, b*d) => 28
            monomialIdeal (a*b, a*c, a*d, b*d, c*d) => 22
            monomialIdeal (a*b, a*c, a*d, c*d) => 30
            monomialIdeal (a*b, a*c, b*c) => 37
            monomialIdeal (a*b, a*c, b*c*d) => 35
            monomialIdeal (a*b, a*c, b*c, a*d) => 26
            monomialIdeal (a*b, a*c, b*c, a*d, b*d) => 33
            monomialIdeal (a*b, a*c, b*c, a*d, b*d, c*d) => 31
            monomialIdeal (a*b, a*c, b*c, a*d, c*d) => 28
            monomialIdeal (a*b, a*c, b*c, b*d) => 29
            monomialIdeal (a*b, a*c, b*c, b*d, c*d) => 23
            monomialIdeal (a*b, a*c, b*c, c*d) => 26
            monomialIdeal (a*b, a*c, b*c, d) => 29
            monomialIdeal (a*b, a*c, b*d) => 28
            monomialIdeal (a*b, a*c, b*d, c*d) => 25
            monomialIdeal (a*b, a*c, c*d) => 20
            monomialIdeal (a*b, a*c, d) => 27
            monomialIdeal (a*b, a*d) => 30
            monomialIdeal (a*b, a*d, b*c*d) => 31
            monomialIdeal (a*b, a*d, b*d) => 30
            monomialIdeal (a*b, a*d, b*d, c*d) => 25
            monomialIdeal (a*b, a*d, c*d) => 34
            monomialIdeal (a*b, b*c) => 25
            monomialIdeal (a*b, b*c*d) => 30
            monomialIdeal (a*b, b*c, a*c*d) => 30
            monomialIdeal (a*b, b*c, a*d) => 27
            monomialIdeal (a*b, b*c, a*d, b*d) => 32
            monomialIdeal (a*b, b*c, a*d, b*d, c*d) => 32
            monomialIdeal (a*b, b*c, a*d, c*d) => 32
            monomialIdeal (a*b, b*c, b*d) => 42
            monomialIdeal (a*b, b*c, b*d, a*c*d) => 35
            monomialIdeal (a*b, b*c, b*d, c*d) => 31
            monomialIdeal (a*b, b*c, c*d) => 24
            monomialIdeal (a*b, b*c, d) => 34
            monomialIdeal (a*b, b*d) => 41
            monomialIdeal (a*b, b*d, a*c*d) => 43
            monomialIdeal (a*b, b*d, c*d) => 31
            monomialIdeal (a*b, c) => 32
            monomialIdeal (a*b, c*d) => 26
            monomialIdeal (a*b, c, a*d) => 36
            monomialIdeal (a*b, c, a*d, b*d) => 36
            monomialIdeal (a*b, c, b*d) => 38
            monomialIdeal (a*b, c, d) => 38
            monomialIdeal (a*b, d) => 16
            monomialIdeal (a*c*d, b*c*d) => 27
            monomialIdeal (a*c, a*b*d) => 30
            monomialIdeal (a*c, a*b*d, b*c*d) => 34
            monomialIdeal (a*c, a*b*d, c*d) => 26
            monomialIdeal (a*c, a*d) => 32
            monomialIdeal (a*c, a*d, b*c*d) => 39
            monomialIdeal (a*c, a*d, b*d) => 16
            monomialIdeal (a*c, a*d, b*d, c*d) => 23
            monomialIdeal (a*c, a*d, c*d) => 27
            monomialIdeal (a*c, b*c) => 21
            monomialIdeal (a*c, b*c*d) => 24
            monomialIdeal (a*c, b*c, a*b*d) => 24
            monomialIdeal (a*c, b*c, a*b*d, c*d) => 25
            monomialIdeal (a*c, b*c, a*d) => 37
            monomialIdeal (a*c, b*c, a*d, b*d) => 28
            monomialIdeal (a*c, b*c, a*d, b*d, c*d) => 35
            monomialIdeal (a*c, b*c, a*d, c*d) => 28
            monomialIdeal (a*c, b*c, b*d) => 28
            monomialIdeal (a*c, b*c, b*d, c*d) => 18
            monomialIdeal (a*c, b*c, c*d) => 17
            monomialIdeal (a*c, b*c, d) => 26
            monomialIdeal (a*c, b*d) => 30
            monomialIdeal (a*c, b*d, c*d) => 23
            monomialIdeal (a*c, c*d) => 20
            monomialIdeal (a*c, d) => 27
            monomialIdeal (a*d, b*c*d) => 30
            monomialIdeal (a*d, b*d) => 32
            monomialIdeal (a*d, b*d, c*d) => 49
            monomialIdeal (a*d, c*d) => 52
            monomialIdeal (a, b) => 27
            monomialIdeal (a, b*c) => 13
            monomialIdeal (a, b*c*d) => 17
            monomialIdeal (a, b*c, b*d) => 11
            monomialIdeal (a, b*c, b*d, c*d) => 21
            monomialIdeal (a, b*c, c*d) => 12
            monomialIdeal (a, b*c, d) => 26
            monomialIdeal (a, b*d) => 10
            monomialIdeal (a, b*d, c*d) => 15
            monomialIdeal (a, b, c) => 38
            monomialIdeal (a, b, c*d) => 20
            monomialIdeal (a, b, c, d) => 31
            monomialIdeal (a, b, d) => 34
            monomialIdeal (a, c) => 30
            monomialIdeal (a, c*d) => 21
            monomialIdeal (a, c, b*d) => 26
            monomialIdeal (a, c, d) => 37
            monomialIdeal (a, d) => 15
            monomialIdeal (b*c, a*b*d) => 26
            monomialIdeal (b*c, a*b*d, a*c*d) => 26
            monomialIdeal (b*c, a*b*d, c*d) => 29
            monomialIdeal (b*c, a*c*d) => 21
            monomialIdeal (b*c, a*d) => 29
            monomialIdeal (b*c, a*d, b*d) => 31
            monomialIdeal (b*c, a*d, b*d, c*d) => 29
            monomialIdeal (b*c, a*d, c*d) => 47
            monomialIdeal (b*c, b*d) => 50
            monomialIdeal (b*c, b*d, a*c*d) => 35
            monomialIdeal (b*c, b*d, c*d) => 30
            monomialIdeal (b*c, c*d) => 14
            monomialIdeal (b*c, d) => 23
            monomialIdeal (b*d, a*c*d) => 33
            monomialIdeal (b*d, c*d) => 42
            monomialIdeal (b, a*c) => 35
            monomialIdeal (b, a*c*d) => 43
            monomialIdeal (b, a*c, a*d) => 39
            monomialIdeal (b, a*c, a*d, c*d) => 27
            monomialIdeal (b, a*c, c*d) => 22
            monomialIdeal (b, a*c, d) => 30
            monomialIdeal (b, a*d) => 37
            monomialIdeal (b, a*d, c*d) => 32
            monomialIdeal (b, c) => 26
            monomialIdeal (b, c*d) => 26
            monomialIdeal (b, c, a*d) => 35
            monomialIdeal (b, c, d) => 34
            monomialIdeal (b, d) => 33
            monomialIdeal (c, a*b*d) => 28
            monomialIdeal (c, a*d) => 35
            monomialIdeal (c, a*d, b*d) => 25
            monomialIdeal (c, b*d) => 23
            monomialIdeal (c, d) => 30
            monomialIdeal 1 => 21
            monomialIdeal a => 10
            monomialIdeal b => 39
            monomialIdeal c => 20
            monomialIdeal d => 40
            monomialIdeal(a*b) => 35
            monomialIdeal(a*b*c) => 26
            monomialIdeal(a*b*c*d) => 32
            monomialIdeal(a*b*d) => 32
            monomialIdeal(a*c) => 26
            monomialIdeal(a*c*d) => 38
            monomialIdeal(a*d) => 42
            monomialIdeal(b*c) => 14
            monomialIdeal(b*c*d) => 29
            monomialIdeal(b*d) => 22
            monomialIdeal(c*d) => 43

o10 : Tally

i11 : J

o11 = {monomialIdeal (a*b, a*c, b*c, c*d), {a*b, a*c, b*c, c*d}, {b*d, a*d,
      -----------------------------------------------------------------------
      c}}

o11 : List

Ways to use randomSquareFreeStep :

randomSquareFreeStep(Ideal)
randomSquareFreeStep(List)
randomSquareFreeStep(MonomialIdeal)

For the programmer

The object randomSquareFreeStep is a method function with options.