This is an interesting construction, but I am wondering if a uniform
distribution for all polynomials of specified degree < d, with a
specified number of terms, is the most natural one to give, and how
grave the impact is on efficiency. (Depending on the coefficient
ring, this goal may not even be achievable).
Also, rather than specifying maximum degree/number of terms, it might
make more sense (and be much after to implement) to use a
distribution with an expected degree and/or number of terms.
- Robert
On Oct 25, 2007, at 2:04 PM, Martin Albrecht wrote:
>
> Hi,
>
> after realizing that Steffen and I share an office at RHUL I
> discussed this
> thing with him today for some time.
>
> The improved implementation (#980) still doesn't produce random
> polynomials
> for two reasons: The code does not seem to produce monomials of
> degree up to
> $d$ uniformly at random.
>
> deg = degree
> for i in range(n):
> exp = ZZ.random_element(0,deg+1)
> if exp <= deg:
> exps.append(exp)
> deg -= exp
> else:
> exps.append(0)
>
> E.g. in this implementation the probability that a (0,0,0....,0,0,0,0)
> exponent tuple is produced is (1/(d+1))^n because we have to hit 0
> out of d+1
> possibilities n times.
>
> However, if we write down a list M of all possible monomials of
> degree <= d
> this L has length
>
> #M = \sum_{i=0}^d ( binomial(n + i -1, i ) )
>
> and the probability to hit the constant monomial is (1/#M).
> Plugging in some
> numbers: n = 10, d = 2
>
> implemented: (1 / 3)^10
> uniformly at random: 1 / ( 1 + 10 + 55)
>
> To fix this we need to consider choosing monomials of a degree D at
> random
> which can be done by finding a bijection between the integers up to
> some N =
> binomial( n + d -1 , d ) and the set of compositions. Mike (mhansen)
> provided this code to do it:
>
> def random_monomial(n, degree):
> total = binomial(n+degree-1,degree)
> random_index = randint(0, total)
> comb = choose_nk.from_rank(random_index, n+degree-1, n-1)
> monomial = [ comb[0] ]
> monomial += map(lambda i: comb[i+1]-comb[i]-1, range(n-2))
> monomial.append( n+degree-1-comb[-1]-1 )
>
> and this explanation:
>
> http://www.cse.ucsd.edu/classes/sp07/cse21/Ron_Notes/starbin.pdf
>
> Side note: This code does not correctly for me as it can run into
> an infinite
> loop which probably can be tracked down to _comb_largest with a
> parameter 'w'
> < 0.
>
> This construction is random if the random number generator
> ("randint") is
> random. Btw. how random is randint?
>
> Now we need to choose a degree d at random. We cannot do this
> uniformly at
> random because the classes (associated to the degree) are not of
> the same
> size (e.g., 1, 10, 55 above) We don't want to choose classes
> uniformly but
> elements and again mhansen beat me to an implementation:
>
> def random_monomial_less_than_degree(n, degree):
> #Get the counts of the total number of monomials
> #of degree d
> counts = [1] #degree 0
> for d in range(1, degree):
> counts.append(binomial(n+d-1, d))
> total = sum(counts)
> #Select a random one
> random_index = randint(0, total-1)
> #Figure out which degree it corresponds to
> d = 0
> while random_index >= counts[d]:
> random_index -= counts[d]
> d += 1
>
> #Generate the corresponding monomial
> comb = choose_nk.from_rank(random_index, n+d-1, n-1)
> monomial = [ comb[0] ]
> monomial += map(lambda i: comb[i+1]-comb[i]-1, range(n-2))
> monomial.append( n+d-1-comb[-1]-1 )
> return monomial
>
> However, an implementation without the need to explicitly construct
> the counts
> array would be good.
>
> With this in place one can start thinking about constructing random
> polynomials from random monomials. For this we don't have an
> elegant solution
> yet. We (=Steffen and me) distinguish two cases:
>
> Let #T be the number of terms requested.
>
> If $MM is the number of all possible monomials of all allowed
> degree <= d and
> #T < $MM/2 then we generate #T _distinct_ monomials as above and
> assign
> random coefficients (possibly zero). The overhead of requesting
> distinct
> monomials should be at most 2 because if we are looking for 1/2 of
> the search
> space, every second choice should be a double.
>
> If #T > $MM/2 we generate all possible monomials, pick #T and assign
> coefficients. This should be more effective for this case because
> we don't
> generate any doubles (but memory is more precious than cpu cycles,
> so we
> might have to move that barrier up a bit).
>
> Does this sound like a sensible approach?
>
> Martin
>
> --
> name: Martin Albrecht
> _pgp: http://pgp.mit.edu:11371/pks/lookup?op=get&search=0x8EF0DC99
> _www: http://www.informatik.uni-bremen.de/~malb
> _jab: [EMAIL PROTECTED]
>
>
>
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~----------~----~----~----~------~----~------~--~---
