[
https://issues.apache.org/jira/browse/RNG-195?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=18106798#comment-18106798
]
Alex Herbert commented on RNG-195:
----------------------------------
The DirichletSampler uses gamma samples with concentration alpha and the other
parameter theta=1.0. I have investigated the values of alpha that can cause an
issue with the Dirichlet sample.
Note that the Dirichlet distribution has at least 2 gamma samples that are
added together. The maximum number of samples to add is limited by an integer
(2^31 - 1).
h2. Small Alpha
Uses the Ahrens Dieter algorithm.
Statistics from 100000 random gamma samples counting the number of samples
equal to zero:
||log2(alpha)||alpha||p(x=0)||
|-1|0.500|0.0|
|-2|0.250|0.0|
|-3|0.125|0.0|
|-4|0.0625|0.0|
|-5|0.0313|0.0|
|-6|0.0156|1.0E-5|
|-7|0.00781|0.00294|
|-8|0.00391|0.05539|
|-9|0.00195|0.23584|
|-10|0.000977|0.48614|
|-11|0.000488|0.69803|
|-12|0.000244|0.83681|
|-13|0.000122|0.91541|
|-14|6.10e-05|0.95657|
|-15|3.05e-05|0.97797|
|-16|1.53e-05|0.98897|
|-17|7.63e-06|0.99446|
|-18|3.81e-06|0.9971|
|-19|1.91e-06|0.99848|
|-20|9.54e-07|0.99924|
|-21|4.77e-07|0.99967|
|-22|2.38e-07|0.99986|
|-23|1.19e-07|0.99994|
|-24|5.96e-08|0.99997|
|-25|2.98e-08|0.99998|
|-26|1.49e-08|0.99998|
|-27|7.45e-09|0.99999|
|-28|3.73e-09|1.0|
|-29|1.86e-09|1.0|
|-30|9.31e-10|1.0|
For example as the alpha parameter drops below 1e-4 there is above 90% chance
of generating a zero. For k concentration parameters the chance of all zeros is
p^k. This will always be true when the concentration parameter is very small.
The algorithm internally uses some predefined constants:
{noformat}
1 / alpha
1 + alpha / Math.E{noformat}
The second constant is multiplied by u in [0, 1) to generate p and compared to
1 to make a branch decision. When alpha / Math.E is less than machine precision
this branch is always chosen as p <= 1. Inside the branch the sampled p is
raised to the power of 1/alpha. This will tend towards zero for small alpha,
e.g.
{noformat}
Math.pow(Math.nextDown(1.0), 1e19) == 0{noformat}
An alpha value of 1e-19 cannot possibly generate a non-zero gamma sample.
However higher values can and there is no definitive cut-off.
h2. Large Alpha
Uses the Marsaglia Tsang algorithm.
Statistics from 100000 random gamma samples counting the number of samples
equal to alpha:
||log2(alpha)||alpha||p(x=alpha)||
|90| 1.24e+27| 0.00681|
|91| 2.48e+27| 0.00984|
|92| 4.95e+27| 0.01395|
|93| 9.90e+27| 0.02002|
|94| 1.98e+28| 0.02858|
|95| 3.96e+28| 0.04056|
|96| 7.92e+28| 0.05705|
|97| 1.58e+29| 0.08035|
|98| 3.17e+29| 0.11332|
|99| 6.34e+29| 0.15954|
|100| 1.27e+30| 0.22434|
|101| 2.54e+30| 0.31304|
|102| 5.07e+30| 0.42942|
|103| 1.01e+31| 0.57461|
|104| 2.03e+31| 0.73523|
|105| 4.06e+31| 0.86778|
|106| 8.11e+31| 0.95787|
|107| 1.62e+32| 0.99845|
|108| 3.25e+32| 1.0|
|109| 6.49e+32| 1.0|
|110| 1.30e+33| 1.0|
Internally the algorithm uses 2 constants:
{noformat}
dOptim = alpha - 1/3
cOptim = 1/3 / sqrt(dOptim)
{noformat}
When alpha -> 3e32 then
* dOptim = alpha
* cOptim = (1 / (3 * sqrt(alpha)) -> 0
The sample is generated from a random Gaussian sample using:
{noformat}
oPcTx = 1 + cOptim * Gauss(X)
v = oPcTx * oPcTx * oPcTx{noformat}
If conditions are satisfied the sample is:
{noformat}
theta * dOptim * v{noformat}
When alpha is large v=1, dOptim=alpha and theta=1 results in all samples being
alpha.
Note that in the DirichletSampler as at least 2 gamma samples are added. So if
alpha=2^1023 then the sum will overflow to infinity and the sampler enters
unbounded recursion. This situation can be predicted: if all concentration
parameters are large then if they sum to infinity then the sampler will
generate stack overflow error for every sample. However the practical threshold
for a valid sampler is much lower. A sampler created using a large but not
excessive alpha will return constant samples, e.g.:
{code:java}
UniformRandomProvider rng = ...
int k = 2;
double alpha = 1e35;
// Generates constant output
DirichletSampler.symmetric(rng, k, alpha) {code}
h2. Conclusion
The DirichletSampler uses an underlying gamma sampler.
The gamma sampler will output all zeros when the alpha parameter is small (<
1e-19), and is very likely to output zeros with alpha < 1e-10. This will
invalidate the Dirichlet sample algorithm. The current behaviour is a stack
overflow error on sample generation.
The gamma sampler will output all samples equal to the alpha parameter when
alpha > 1e33 (approx). The gamma sampler, and consequently the
DirichletSampler, is thus invalid with large alpha as it does not sample a
distribution but a constant. When the sum of Dirichlet concentration parameters
is infinite then the sampler will generate a stack overflow error, otherwise it
will output a constant.
Since the sampler uses a source of randomness it is possible that as these
alpha limits are approached the rejection algorithm can still generate samples.
Adding a hard limit to the upper bound for alpha would be possible to avoid the
sampler emitting constant values. But this could still occur below the hard
limit and provides little value beyond catching a grossly misconfigured
sampler. The same goes for the lower limit. It may avoid an obviously broken
sampler from unbounded recursion, but this would still be possible above the
hard limit.
There are at least two options:
# Impose hard limits where the sampler is known to be non-functional. Document
it could still be non-functional within these limits.
# Leave the behaviour as is; document the known limits and advise the user to
expect stack overflow errors if they approach these extremes.
> DirichletSampler has unbounded recursion and stack overflow for small alpha
> parameters
> --------------------------------------------------------------------------------------
>
> Key: RNG-195
> URL: https://issues.apache.org/jira/browse/RNG-195
> Project: Commons RNG
> Issue Type: Bug
> Components: sampling
> Affects Versions: 1.3
> Reporter: Alex Herbert
> Priority: Minor
>
> The DirichletSampler uses a rejection method to sample. A gamma sample is
> obtained for each concentration parameter alpha. The vector of k samples is
> normalised to unit length. When the sum of the samples is zero or infinite
> the sample cannot be normalised and recursion occurs to generate another
> sample.
> If the alpha parameters are very small the gamma samples are always zero and
> a stack overflow error occurs.
> Note that recursion is expected to be used to indicate a non-functional RNG
> underlying the sampler. It was not the original intention to use a stack
> overflow error to indicate a badly parameterised sampler where samples are
> impossible.
> Issue identified using a security scan.
--
This message was sent by Atlassian Jira
(v8.20.10#820010)
