Better API for generating random elements

[fingolfin]

(I am aware of issue #100 but I think it is about something else)

I think we need a better interface for rand on e.g. polynomial rings and other of our datastructures. The problem with the current system is that it is kinda hard to discover how to use it, and even if you know, it requires a lot of guesswork. At least in my experience.

Here is an example of what we have right now:

julia> R, x = polynomial_ring(QQ, "x" => 1:3);

julia> rand(R, 1:3, 4:6, 7:9)
x[1]^6*x[2]^6*x[3]^4 + x[1]^6*x[2]^5*x[3]^6 + 9//8*x[1]^4*x[2]^4*x[3]^4

julia> rand(R, 1:3, 4:6, 7:9)
x[1]^5*x[2]^6*x[3]^5

julia> rand(R, 1:3, 4:6, 7:9)
x[1]^4*x[2]^4*x[3]^5 + 7//9*x[1]^4*x[2]^4*x[3]^4

So rand has to be called with three ranges (or vectors): the first indicates the acceptable number of terms, the second indicates the allowed range of degree of each variable, the final one is used to determine which denominators and numerators the coefficients may have (I think).

The logic is clear in principle: the first two ranges control what is specific to polynomial rings; everything after controls the randomness of the coefficient ring. This can be nested then:

julia> S, y = R[:y1, :y2];

julia> rand(S, 1:3, 4:6, 7:9, 10:12, 13:15)
(x[1]^12*x[2]^12*x[3]^12 + 15//13*x[1]^12*x[2]^11*x[3]^10 + 14//15*x[1]^11*x[2]^12*x[3]^12 + 14//15*x[1]^11*x[2]^10*x[3]^11 + 405//182*x[1]^10*x[2]^12*x[3]^11 + 14//15*x[1]^10*x[2]^10*x[3]^12 + 421//210*x[1]^10*x[2]^10*x[3]^11)*y1^6*y2^6 + (x[1]^12*x[2]^11*x[3]^11 + 379//182*x[1]^12*x[2]^11*x[3]^10 + x[1]^11*x[2]^12*x[3]^12 + x[1]^11*x[2]^10*x[3]^12 + 14//15*x[1]^11*x[2]^10*x[3]^11 + x[1]^10*x[2]^11*x[3]^12 + 14//15*x[1]^10*x[2]^10*x[3]^11)*y1^5*y2^6 + (15//13*x[1]^12*x[2]^12*x[3]^12 + 391//210*x[1]^12*x[2]^12*x[3]^10 + 15//14*x[1]^12*x[2]^11*x[3]^11 + x[1]^11*x[2]^11*x[3]^12 + x[1]^11*x[2]^11*x[3]^11 + 13//15*x[1]^10*x[2]^11*x[3]^10 + 13//14*x[1]^10*x[2]^10*x[3]^12)*y1^4*y2^6

(In case you wonder about 405//182: well if some terms are generated multiple times with different coefficients, they are of course merged...)

While this is "logical", I still find it difficult to use, and worse: to read in code using it without trying it out to verify.

It is also somewhat limited, e.g. there is no way to e.g. "get a random polynomial where all terms have total degree 4 to 6, and coefficients in the range -10:10".

I propose that we design and switch to a new system that relies on kwargs to become more "self documenting". Here is a first example of how it might look like in the above examples; it's just a first draft to convey the idea and get the discussion going, so I am happy to hear if others have different and better ideas :-)

I'll give these in the form of examples, without explaining much what they mean, in the hopes that they are mostly self-documenting. These lists are not meant to be exhaustive nor do I insist all variants are needed (or even sensible), I am more concerned about the "flavor" here.

For QQ we might offer:

• rand(QQ; denominator=1:10, numerator=-5:5)
• rand(QQ; denominator=10, value=-5:5)
• rand(QQ; nbits=1:10)
• rand(QQ; nbits=1:10)
• rand(QQ, 7:9) -- for backwards compatibility??
• rand(QQ) -- either returns random values wrt some documented default values, or produces an error explaining how to use it / pointing to the docs

For GF(p) we might offer

• rand(GF(p))
• rand(GF(p), non_zero=true)
• rand(GF(p), value=1:p-1)

For a polynomial ring, specify using coefficients the parameters to pass on to the rand function on the coefficient ring. I envision this is either a named tuple (passed as kwargs), or a tuple (passed as regular args) or some other value (passes as a single arg).

So e.g. if the coefficient ring is QQ then we might offer:

• rand(R; total_degree=4:6, coefficients=(denominator=1:10, numerator=-5:5))
• rand(R; nterms=1:3, degree=4:6, coefficients=7:9) -- invokes rand(QQ, 7:9)
• rand(R) -- either uses some default settings, or raises a helpful error

If the coefficient ring is another polynomial ring over QQ, one can recurse this, e.g.:

• rand(R; total_degree=4:6, coefficients=(nterms=1:3, degree=4:6, coefficients=7:9))

Similar strategies should work for matrix spaces, matrix algebras, group rings, Lie algebras, matrix groups, power series, ...

Perhaps there are structures where this is not a good fit, in that case it would be good to know about them, so please list them.

[fingolfin]

Let me say that a related but somewhat separate issue is that in many cases (for experiments but also algorithms) what one wants is not just a fixed way to get random elements, but really one needs a "random process" which is a stateful object that can be asked to produce additional random elements, but which can uses its state to e.g. increase the parameters over time: e.g. maybe the first 10 random polynomials should have total degree <= 5, and after that this is raised to <= 10, or whatever. But I think this could be mostly discussed in a separate issue (if someone wants to open one...?), here we should only discuss it if it has implications for the API design...

[fingolfin]

Also, of course this can also be made to work with even more flexible "sampler" objects as discussed in the Random documentation, e.g. via rand(QQ, sampler=some_object) respectively. rand(poly_ring, sampler=s1, coefficients=(sampler=s2,)).

[lgoettgens]

Do we wanna get consistent with julia Base again in the sense that one can pass dims to generate multiple elements at once? For more details, see https://docs.julialang.org/en/v1/stdlib/Random/#Base.rand

[thofma]

I like it. But I would be in favour of keeping rand(ZZ, -10:10) (forever, and not just as a deprecated method).

[fieker]

On Sun, Oct 08, 2023 at 03:52:22PM -0700, Lars Göttgens wrote: Do we wanna get consistent with julia Base again in the sense that one can pass `dims` to generate multiple elements at once? For more details, see https://docs.julialang.org/en/v1/stdlib/Random/#Base.rand -- Reply to this email directly or view it on GitHub: #2896 (comment) You are receiving this because you are subscribed to this thread. Message ID: ***@***.***>

I think we need to think carefully about what we want here. RandomExtensions was partly, in the end, supported even by us, but seems to be dead(?). Furthermore, we receive too many complains about rand and friends not being useable by "normal" users or not giving what they want. Julia is cheating: mostly they are working in a world where - the type contains "everything" - most "domains" are finite (Float, Int) So the complexity is in array dimensions and, maybe, esoteric distributions. For us, both assumptions are wrong. How do I create a random polynomial with entries in a number field? Or worse: a matrix of those? A random prime(ideal)? A random ideal? A random (sub)group? I don't know a "good" answer, but looking at some of the complaints, our current approach of adding parameters to rand following the construction of the parent is very hard to use, understand and predict. Any good suggestions anyone?

[fieker]

Max is going in the right direction (I read it only afterwards), but as this requires some serious work to roll out, we should think about it carefully.
E.g. named tuples vs Dicts, kwargs vs pairs, some "small random slowly increasing" thing, how to deal with recursion? With generics?
E.g. a polynomial aver a matrix algebra would need nterms (as a poly) and degrees (as a poly), but coefficients here has to be about entries (of the matrices) va a polynomial over the integers where only a range would be required for the coeffs...

[fieker]

Or an interface: random_params(some parent)? To return a dict, as recursive as neccessary, with all(?) parameters set to s.th. sensible(?) Possibly, for the connosieur with the possibility to specify some sizes directly?

And some random_grow(Dict from above) that increases all(?) paramters sensibly(?)?

[fingolfin]

some notes from discussion today:

• @jankoboehm requests random homogenous polynomials specified via (total) degree plus a sparsity (0-1) value
• request for random graded homomorphism of degree $d$ between e.g. graded polynomial rings (or more general) -> I think that should be done via rand(homSpace) but then homSpace needs to be constructed lazily (and of course if the coefficient rings are infinite, then one will have to restrict those, too)
• provide "sampler" objects (to use Julia lingo) that retain state and can use this to e.g. produce polyomials with slowly growing degrees and coefficients "over time" (or anything else, whatever we fancy -- and we can add more different kinds of these samplers over time)
  • samplers can be composed: one can e.g. specify a sampler as argument for coefficients
• immediate use case for random polynomials (suggested by @HereAround): generic section of a line bundle on a toric space; code for this already exists
  • more or less the same in F-Theory tools
• Claus points out it could be useful to query a ring to know "what are all the parameters one needs to specify to get "something random"
• also need things like "random symmetric matrix of graded polynomials" -> the "symmetric" is not an aspect of the parent...
• ...

[JohnAAbbott]

Some comments:

1. the interface for generating a random rational number (from a
   finite set specified by some parameters, somehow): do we desire
   that the values in this finite set are selected with uniform
   probability (or "near uniform")? Including for QQ?

2. I like the idea of specifying the randomization parameters via
   a Dict because:
   (-) this is an (single) object, so can easily be passed around
   (-) it can have a recursive structure (i.e. some field could be
   another such Dict)

3. Claus Fieker mentioned also the desire to allow the "range" for
   the random values to grow. Doing this in a fully general way
   seems to be tricky. We could consider an alternative: namely,
   state the maximum number of values which may be generated using
   the parameters in the Dict -- an exception is thrown if more are
   requested. I assume that Dicts are "passed by reference" in Julia.

[JohnAAbbott]

Some more thoughts about (pseudo-)random generators.

I wish to revise my previous endorsement of Dict.

I now think that a better approach is to use "iterators"
(similar to Julia's "samplers"): that is, objects which
produce a stream of values sampled from some population
with a specified distribution. The advantage of such
an approach is that all information about the type of
randomness desired is captured in a single object which
can be passed as parameter, stored in a variable etc.

I think we should avoid the Julia name "samplers" because,
while notionally similar, our objects are not compatible.

We can also have "iterators" which act as filters: i.e.
they obtain random values from some other "iterator",
and then filter them in some way. A very simple filter
could simply throw an error after k values have been
emitted. A more sophisticated example could be one to
insist that only distinct random values are produced:
this entails memorizing the values already emitted, and
emiting only future values which are not already in the
table. We may want to use a heuristic to detect "value
exhaustion": e.g. if we have already emitted k different
values, and no new, distinct value has been generated in
5*k attempts then we declare exhaustion (e.g. throw error).

A constructor for an iterator (random generator) may need
other iterators as argument. For instance, a random dense
matrix with entries being sampled from an iterator which
generates appropriate ring elements. Thus the user must
creat first the desired random-ring-element-iterator, and
then pass it as argument to the random-matrix-iterator.

I envisage different implementations of random-matrix-iterator
for "each entry is i.i.d", "matrix is symmetric", "matrix
is anti-symmetric" -- these would most likely all expect an
arument specifying random-ring-element-iterator. CoCoA
has functions for generating random unimodular matrices,
and for generating random sparse non-singular {0,1} matrices;
the latter would not require a random-ring-element-iterator
(internally it uses a biased boolean iterator where the bias
is computed inside the function).

For integers we need to specify a range: this could be done
either via bit-size, or by a given limit value (inclusive?).
We may want to specify "non-negative" or "symmetric range"
(e.g. via two different types of iterator, or maybe a flag?)
I would not initially cater for the case of ranges other than
[0,max] or [-max,+max].

For rationals it becomes trickier. I believe it is important
to keep the interface simple. As for integers there could be
"symmetric" and "non-negative" options/variants. Range bounds
could be either [max_numer, max_denom] or [max_value, max_denom];
the latter is a subset of [max_value*max_denom, max_denom].
How important is "near uniformity" here? Generating random
numerator & denominator then cancelling common factors seems to
favour simpler rationals!

A quick comment about multi-threading: as far as I know
pseudo-random generators depend on some hidden values of
global lifetime; I have no idea how to achieve reproducible
behaviour in a genuine multi-threaded context.

[thofma]

I think we already have the sampler/iterator and I am afraid no one ever uses them. It is very cumbersome to write r = RandomGenerator(QQ, numerator = bla, denominator = blub); x = sample(r). Most people will want to just write x = rand(...).

The generator/iterator thingy will probably be useful for the idea from #2896 (comment). But as I said, this is in principle already how it is implemented, since we use https://github.com/JuliaRandom/RandomExtensions.jl (but we did not implement all the things that you mentioned).

[JohnAAbbott]

I'm convinced about having a short form for random integers -- I would expect that random integers are much more widely used than other types of random value. I can see strong arguments in favour of "symmetric" and "non-negative" random integers; personally I'm not convinced by random integers from any other sort of range -- I know we could easily implement such functions, but that doesn't necessary mean that it is a good idea. In CoCoA when we had doubts about easy-to-implement extensions, we would often implement them, and then comment out the implementation with a note saying that it should be activated if there is sufficient demand for it (or presumably deleted if no demand appears before some deadline...).
I remain to be convinced that it is a good idea to offer a short form for any other type of random value.

[JohnAAbbott]

How useful are random rationals?

Consider generating a random element in QQ[x]; what do we expect to obtain?
I generated a random polynomial of degree 10 with random QQ coefficients
whose denominators are uniform in 1:1000 and numerators are in 1:2d, where
d is the chosen denominator. The factorization of the polynomial has large
coefficients:
205421476497987921600*x^10 + 638715825226667485440*x^9 +...+ 284366318688030205440
with numerical factor 1/372918988104039611520 whose denominator is the LCM
of the denominators generated.

I suggest that this is not the desired outcome, and that it is likely
more desirable to create a polynomial with integer coefficients in the
given range, and then a single common denominator for the whole polynomial.

Here's another "contrived" example:
Let M_ZZ be a 100x100 matrix with random ZZ entries in -2^99 to 2^99.
Let M_QQ be a 100x100 matrix with random QQ entries with 50-bit numer & denom.
So each matrix has 100 random bits in each entry.
Let V_ZZ be the inverse of M_ZZ, and V_QQ be the inverse of M_QQ.
The entries in V_ZZ have about 10000 bit numerators, and a single 10000 bit
common denominator for the entire matrix.
The entries in V_QQ are around 40-50 times "bigger": e.g. each entry has
numerator with about 450000 bits, and denominator with also about 450000 bits.
So even though the matrices appear to contain the same amount of randomness,
M_QQ is somehow "much more complicated".

Summary: we need to decide what we really want when we talk of random rationals!

[YueRen]

If it is difficult to agree on what rand(QQ) (and others) should return, I vote for @fingolfin suggestion that we over rand(QQ; ....) with optional arguments and that rand(QQ) should raise an error.