Does machine learning on random situations require a cryptographically-secure random number generator?

I'm working on a project involving using machine learning to figure out optimal strategies to play board games, and after a few rolls of the virtual dice a thought hit me:

For the games that involve randomness, such as those with dice, there can be significant value in predicting future dice rolls. And the rolls themselves are not, of course, truly random. Given this, do I reasonably have to fear that after learning through all the low-hanging fruit as far as gameplay is concerned, my ML models may try to optimize for the specific random number generator I am using, such as by predicting future rolls?

If so, would a cryptographically-secure random number generator be necessary to obviate this, or are there computationally-cheaper alternatives that offer the same protection from ML "gaming the system" in this way?

• A sloppy pRNG for online poker undermined the security of the game. Importantly, you only have to know how to deduce the pRNG seed to get a game-breaking edge; you don't have to use machine learning to approximate the modulus function itself. pdfs.semanticscholar.org/b0b5/… – Sycorax Aug 11 '20 at 20:56

Edit: My original answer below is mostly informal, but I want to address some of the comments in a more technical and hopefully convincing manner. Please see the technical appendix for these details.

Does machine learning on random situations require a cryptographically secure random number generator, or in other words, is it reasonable to fear that your machine learning algorithm will learn how to predict the output of your pseudo-random number generator (PRNG)? Generally, no.

Could a machine learning model such as a neural network emulate a PRNG? By this, I mean: could the function $$f$$ that produces the sequence of pseudo-random numbers be in the class of functions $$V$$ that the machine learning model is capable of representing. Possibly, depending on the model in question.

Could a capable machine learning model accidentally be trained from data generated by the PRNG to predict the output. Almost certainly not, though the probability of this is non-zero. Could we successfully create and train a custom machine learning model with the sole purpose of predicting the output of a PRNG? Also probably not, at least not without a great deal of "cheating."

The key point is that even if a machine learning model is capable of representing the PRNG, it has to be capable of finding the right parameters to predict the output of the PRNG. Training a machine learning model to predict the output of a PRNG is an extremely difficult task, bordering on the impossible. To understand why, let's first talk about how a PRNG works.

Pseudo-Random Number Generation

Most PRNGs use some form of congruential algorithm which involves starting with a positive integer $$X_0$$ called the seed and then making a recursive sequence according to a rule similar to

$$X_{n + 1} = g(X_n) \text{ mod } m$$

for some function $$g$$ and constant $$m \in \mathbb{N}$$. There are some slight variations on the method, and of course some methods which are completely different, such as cellular automata-based methods (like Wolfram's Mathematica uses by default). To answer your question, I'm going to focus on one of the simplest PRNGs: the linear congruential method, which uses the function $$g(x) = ax + c$$ for some integer constants $$a$$ and $$c$$. This method is used by the Java programming language, even though it has relatively poor statistically properties. I'm then going to appeal to intuition to claim that, if we don't have to worry about a machine learning algorithm learning how to predict the output of a very simple PRNG with poor statistical properties, we probably don't have to worry about it learning how to predict the output of a fancy PRNG with better statistical properties.

Now, let's consider the actual constants $$a$$, $$c$$, and $$m$$ to use. There are various properties that these need to satisfy to make a good PRNG that I won't discuss (see Donald Knuth's The Art of Computer Programming vol. 2 which is an authoritative treatment of the topic). Let's just consider the constants that Java's PRNG uses as a real-world example. From the source code (on line 173), the values it uses are $$a = 25214903917$$, $$c = 11$$, and $$m = 2^{48} = 281474976710656$$. We also can't forget that in trying to learn the output of the PRNG, the machine learning model will also have to learn the seed $$X_0$$.

Learning the $$x$$ mod $$m$$ Function is Hard

This is the first difficulty that our machine learning model has to surmount. There is already an excellent discussion of this problem on this stackoverflow post which you should read before continuing this post any further.

Hopefully you aren't reading this unless you poked through the linked post. Note that the best solutions use recurrent neural networks (RNN), with the motivation explained in the accepted answer:

Please understand that this solution is a bit of a tongue-in-cheek thing, because it is based on the task-domain knowledge that our target function can be defined by a simple recurring formula on the sequence of input bits.

In reality, if we aren't using domain knowledge for this problem (for instance, if you're designing your model to play a dice game) then the model might not be able to learn the $$x$$ mod $$m$$ function. You could test this by using your model architecture and applying it to this problem directly to see if you can get good results.

Cost Functions and Convexity

Okay, so maybe learning $$x$$ mod $$m$$ is difficult, but as the stackoverflow answer above demonstrates, it is doable. So what's the next hurdle? Let's talk about the training of a model, i.e. the finding of the parameters that best fit the data.

The "magic" of modern machine learning is very much reliant on the fact that convex optimization techniques like gradient descent seem to "just work" even when applied to non-convex optimization problems. They don't work perfectly, and often require a fair amount of tinkering to train properly, but they can still get good results. One of the reasons for this "magic" is that lots of cost functions, while non-convex, aren't that non-convex. For example, your cost function might look something like:

This cost function might look bad at first glance, but notice that it has some degree of regularity/smoothness. You can still tell that the underlying function is continuous because "small" movements along the $$x$$ or $$y$$-axis result in "small" changes in height. You can also pick out a general basin-shaped structure, and it's believable that a convex optimization algorithm with some random perturbations could eventually find the global minimum.

Essentially, a cost function with some regularity might not be convex, but can still be "locally convex" in some sense. This means that gradient descent can find a local minimum if the initial point is within a locally convex "basin." In other words, being close the minimum counts for something, so "partial" correctness can be rewarded.

Indeed, this is the idea behind transfer learning. Finding a good minimum for one task that is sufficiently similar to another task can provide the second task with a good initial point and then convex optimization can fine-tune the result to find a nearby minimum for the second task.

An Experiment

However, the cost function for trying to learn a PRNG has virtually no regularity whatsoever. It shouldn't come as a shock, but the cost function behaves like noise. But don't take my word for it: let's do an experiment to try to predict the output of Java's PRNG.

For this experiment, we're going to cheat as much as possible and still lose. To start with, instead of using some kind of neural network or other machine learning model with a large number of parameters, we're going to use the exact functional form that we know Java's PRNG takes:

$$X_{n + 1} = (a X_n + c) \text{ mod } m$$

which has parameters $$a$$, $$c$$, $$m$$, and $$X_0$$. This completely sidesteps the difficulty of learning $$x$$ mod $$m$$ discussed above. And our model has only four parameters! Modern machine learning algorithms can have hundreds of millions of parameters that require training, so just four should be a piece of cake, right?

Let's make it even easier though: suppose that an oracle (no pun intended) tells us three of four correct parameters for Java's PRNG, and our task is simply to learn the value of the fourth. Learning one parameter can't be that hard, right?

Here's some Julia code to emulate Java's PRNG and to plot an $$\ell_2$$ cost function over each of the four slices of the four slices we get from not knowing one of the four parameters:

using LinearAlgebra: norm
using Plots
theme(:dark)
seed = 12150615 # Date the Magna Carta was signed

# Constants used by Java's linear congruential PRNG
a = 25214903917
c = 11
m = 2^48

"""Generates the next integer in a sequence of pseudo-random_sequence
numbers in a linear congruential sequence."""
function next(x, a, c, m)
return mod(a*x + c, m)
end

"""Generates a random sequence of M random integers from a linear
congruential sequence with the parameters a, c, m, and seed."""
function random_sequence(a, c, m, seed, M)
nums = zeros(Int, M)
nums[1] = seed
for i = 2:M
nums[i] = next(nums[i-1], a, c, m)
end
return nums
end

# Generate Java's random sequence
y = random_sequence(a, c, m, seed, M)

i_values = -200:200 # Range around the correct parameter to test
n_trials = length(i_values)

# Test a neighborhood of the a-values
as = [a + i for i = i_values]
avg_errors = []
for i = 1:n_trials
# Generate another random sequence using random constants a, b, c, and a random seed
y_test = random_sequence(as[i], c, m, seed, M)
avg_error = norm(y_test - y) / M
push!(avg_errors, avg_error)
end
plot(avg_errors, size=(400, 400), legend=false)

# Test a neighborhood of the c-values
cs = [c + i for i = i_values]
avg_errors = []
for i = 1:n_trials
# Generate another random sequence using random constants a, b, c, and a random seed
y_test = random_sequence(a, cs[i], m, seed, M)
avg_error = norm(y_test - y) / M
push!(avg_errors, avg_error)
end
plot(avg_errors, size=(400, 400), legend=false, ylim=(1.145e11, 1.151e11))

# Test a neighborhood of the m-values
ms = [m + i for i = i_values]
avg_errors = []
for i = 1:n_trials
# Generate another random sequence using random constants a, b, c, and a random seed
y_test = random_sequence(a, c, ms[i], seed, M)
avg_error = norm(y_test - y) / M
push!(avg_errors, avg_error)
end
plot(avg_errors, size=(400, 400), legend=false, ylim=(1.145e11, 1.151e11))

# Test a neighborhood of the seed-values
seeds = [seed + i for i = i_values]
avg_errors = []
for i = 1:n_trials
# Generate another random sequence using random constants a, b, c, and a random seed
y_test = random_sequence(a, c, m, seeds[i], M)
avg_error = norm(y_test - y) / M
push!(avg_errors, avg_error)
end
plot(avg_errors, size=(400, 400), legend=false, ylim=(1.147e11, 1.151e11))


So you can clearly see that even with three of the four parameters and the exact functional form known, the cost function still has the form $$c + (\text{noise})$$ where $$c$$ is a constant. In this case, a gradient-descent-type algorithm would compute a gradient of $$0 + (\text{noise})$$. Then gradient descent is simply performing a random walk. While possible that a random walk could converge to correct parameters, it is extremely unlikely given that the size of the space is $$10^{77}$$ (see below). Without any regularity, convex optimization tools are no better than a random walk looking for that one "valley" in the middle of each graph where the correct parameter lies.

Conclusion

It turns out that even with all of this simplification, the last step is still virtually impossible. "Learning" the last parameter boils down to a brute force search over the entire range of possible values for the parameters, because the "magic" of applying convex optimization techniques to train a machine learning model does not help solve a search problem when the cost function does not have any information whatsoever about the direction of even a good local minimum. If you wanted to try every possible 64-bit integer for the four parameters, this would mean searching through $$(2^{64})^4 = 2^{256} \approx 10^{77}$$ combinations. And this is just for a very simple PRNG.

Ultimately, if you really want to alleviate any worries you might have about your particular task, you could always drop the board game aspect and see if your model can learn the output of the pseudo-random dice roll using your programming language's PRNG. Good luck (you're going need it)!

Technical Appendix

First, I want to point out that the function $$x$$ mod $$m$$ being difficult to approximate is more of an interesting side note, relevant mostly for the concern in the original question that a machine learning algorithm might incidentally discover how to predict the output of the PRNG while being trained for some other purpose. The fact is that it is difficult even when this is one's sole purpose. Difficult, but not unreasonably difficult. You don't need to appeal to universal approximation theorem to claim this is possible, because in the linked stackoverflow post from above there are several examples of models that successful approximated $$x$$ mod $$m$$ (albeit with the input in binary-vector representation). So not only was it possible to represent the function $$x$$ mod $$m$$ by a neural network (which is all the UAT guarantees), they were also able to successfully find weights that worked (which is not guaranteed).

Second, what is the technical meaning of the claim that the cost function has the form $$C(x) = \begin{cases} \text{constant} + \text{noise}, & \text{ if } x \ne x^* \\ 0, & \text{ if } x = x^* \end{cases}$$ where $$x$$ denotes the parameters $$x = (a, c, m, \text{seed})$$ and $$x^*$$ denotes the correct parameters? This can be defined technically by picking a radius $$\delta > 0$$ and then computing the average value $$\text{Avg} (\delta, t) = \frac{1}{m(B_\delta (t))}\int_{B_\delta (t)} C(x) dx$$ where $$K$$ can represent either $$\mathbb{Z}^4$$ or $$\mathbb{R}^4$$, $$m$$ is either the Lebesgue measure or the counting measure respectively, and $$B_\delta (t) = \{ x \in K: \| x - t \| < \delta \}$$ is the ball of radius $$\delta$$ centered at $$t \in K$$. Now the claim that $$C = \text{constant} + \text{noise}$$ means that as $$\delta$$ increases, the local average $$\text{Avg} (\delta, t)$$ converges quickly to a constant $$L$$, so long as the true parameters $$x^* \notin B_\delta (t)$$. Here, I say "quickly" to rule out the fact that eventually this limit would be constant after surpassing the bounds of the domain. This definition makes sense even though the "noise" term is technically deterministic. In other words, local averages of $$C$$ are globally constant. Local averaging smooths out the noise term, and what remains is a constant.

Plotted below is a much larger scale experiment on the interval $$[-9 \times 10^{12}, 9 \times 10^{12}]$$ that shows essentially the same phenomenon as before. For this experiment, I only tested the case where the seed is unknown as this experiment took a much longer. Each point here is not the cost function, but the local average $$\text{Avg} (100, t)$$ of the cost function, which smooths out some of the noise:

Here I've plotted the constant as well, which turns out to be roughly $$\text{constant} = 1.150 \times 10^{12}$$

Ultimately, this is a problem for gradient-based optimization methods not because of the noise term per se, but because the cost function is "flat." Of course, when you do add in the noise term, a flat function plus noise makes an extremely large number of local minima that certainly doesn't help the convergence of any gradient-based optimization algorithm.

Moreover, I am well aware that this is an empirical claim, and I can not prove it analytically. I just wanted to demonstrate empirically that the gradient for this function is essentially 0 on average, and contains no information about the direction of $$x^*$$. In Experiment 1, the neighborhood was purposefully small to demonstrate that even if you started close to $$x^*$$, there is no visible gradient pointing in that direction. The four slices of the neighborhood $$B_{200} (x^*)$$ are small, but still don't show a local "basin" (locally approximately convex region) of the sort that gradient-based optimization is good at minimizing in. Experiment 2 demonstrates this same phenomenon on a much larger scale.

The last technical detail I want to touch on is the fact that I am only analyzing the the model and the cost function as functions over a subset of the domain $$\mathbb{Z}^4$$, not over $$\mathbb{R}^4$$. This means that the gradient/derivative is not defined. So how can I claim something about the convergence or non-convergence of a gradient-based method when gradients aren't defined? Well, we can of course try fitting a differentiable model defined on $$\mathbb{R}^4$$ to the data, and compute its derivative, but if the data is already "flat" a model that fits it well will be "flat" as well.

This is not something that I can prove, but I can prove that it is unprovable by constructing a continuously differentiable ($$\mathcal{C}^1$$) interpolation function $$f : \mathbb{R} \to \mathbb{R}$$ to the cost function data $$C(x)$$ that would cause gradient descent to converge to the true global minimizer $$x^*$$ in one step with high probability. This is an absurd example, but it demonstrates that trying to prove that gradient-based algorithms couldn't conceivably work here is impossible. To construct the interpolating function, consider two adjacent points $$n, n+1 \in \mathbb{Z}$$ with cost function values $$C(n)$$ and $$C(n+1)$$. Pick a threshold $$\epsilon > 0$$. Now, on the interval $$[n + \epsilon, n + 1 - \epsilon]$$, we can construct $$f$$ so that a regular gradient-descent step will reach $$x^*$$ in one step, i.e. $$x^* = x - f'(x)$$. This defines an easy differential equation that we can solve as follows: \begin{align} x^* & = x - f'(x) \\ \int x^* dx & = \int x - f'(x) dx \\ x x^* & = \frac{1}{2} x^2 - f(x) + D\\ f(x) & = \frac{1}{2} x^2 - x x^* + D \end{align} for any constant $$D$$. The constant is irrelevant, because regardless of its value, we can still define $$f$$ in such a way on the intevals $$[n, n + \epsilon)$$ and $$(n+1-\epsilon, n+1]$$ to make $$f \in \mathcal{C}^1$$ and such that $$C(n)$$ and $$C(n+1)$$ are the correct values, using splines for instance. This construction can be repeated on all intervals and the results can be stitched together in a $$\mathcal{C}^1$$ manner (using splines again, as one particular method). The result will be a $$\mathcal{C}^1$$ function that interpolates the cost function at all $$n \in \mathbb{Z}$$ (so it fits the data here perfectly well), and one which will converge to $$x^*$$ in one step of the gradient descent algorithm with probability $$1 - 2\epsilon$$. Take $$\epsilon > 0$$ to be as small as desired.

• Except all of these functions are continuous because they are functions on $\mathbb{Z}$. And I don't have a good term for "degrees of convexity" other than to say something about the extent to which a convex optimization algorithm works or does not work on a function. But I think it is clear that a functions like the cost function here cannot be optimized using anything other than brute force. There's no local regularity to take advantage of. – Eric Perkerson Aug 11 '20 at 13:43
• This experiment does not demonstrate that $(aX_n + c) \mod m$ cannot be well-approximated using machine learning -- the deficiency you've found could plausibly be due to a deficiency of experimental design, and unrelated to the modulus function. You'll have to make an appeal to a relevant theorem to conclusively show why machine learning is not appropriate. – Sycorax Aug 11 '20 at 14:06
• Moreover, the universal approximation theorem does apply to discontinuous functions, so long as they are Lebesgue integrable. – Eric Perkerson Aug 11 '20 at 14:10
• Is the modulus function Lebesgue integrable? – Sycorax Aug 11 '20 at 14:14
• There are lots of "bumpy" loss functions. You show a graph of one, but say that it's not "too bumpy." How do we know this? Also, a neural network with dropout has a different bumpy loss every iteration, yet we can still make progress. For lcm modulus, how do we know that this loss function is "too bumpy"? How do we know when a problem has enough "regularity" to be learnable? What does "regularity" of a loss function mean, and how do we know when it's present? – Sycorax Aug 11 '20 at 14:41