[
https://issues.apache.org/jira/browse/RNG-195?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=18126324#comment-18126324
]
Alex Herbert commented on RNG-195:
----------------------------------
I have fixed the issue in the Dirichlet sampler with small or large
concentration parameters. The issue occurs when gamma samples are all zero or
one or more is infinite. The sum cannot be used to normalise the samples to
unit length.
The solution is to modify the gamma sampler to return log samples when
requested. This is a small modification of the gamma sampler to change the
final output function. This either returns the sample x (unchanged behaviour)
or an equivalent log x. This is expected to have no impact on the existing
gamma sampler functionality.
In some cases the log sample can be returned with more accuracy than using log
on the sample value. Infinite real values can be returned as finite using a log
sample which solves the large alpha case. Small alpha are still limited.
When this modification is made the gamma sampler will output finite values
using an alpha parameter between 2^-1018 and double MAX_VALUE.
In order to extend support to all non-zero alpha requires some adjustment. The
small alpha gamma sampler will always return the sample for any alpha below
2^-52 using:
{noformat}
log(p) / alpha : p in [0, 1)
{noformat}
When p is 2^-53 the log value is -36.74. This sets a limit on alpha since
-36.74 / 2^-1019 is -infinity. Other samples with larger p may not overflow.
However as alpha -> 2^-1075 the number of samples that are not -infinity
reduces to almost none. Given log(1 - 2^-50) / 2^-1075 == -infinity, only 7 of
2^53 possible dyadic values of p will not cause overflow.
However since all small alpha use the same output function we can use the
identity:
{noformat}
log(p) / alpha == scale * log(p) / (scale * alpha)
{noformat}
The Dirichlet sampler can detect the largest alpha, and if far below 2^-52
apply scaling to all alpha parameters as long as they remain under 2^-52. The
log gamma samples can be rescaled by the Dirichlet sampler before normalization.
Normalization of log samples is performed by subtracting the maximum log value
from all samples before applying the exp function:
{noformat}
sum { exp( s_i - max(s_i) ) }
{noformat}
The maximum log sample will have a value of 1, all other samples will be
smaller and the sum is finite. This is an application of the log-sum-exp trick
for log domain calculations (see [Log-Sum-Exp
(Wikipedia)|https://en.wikipedia.org/wiki/LogSumExp#log-sum-exp_trick_for_log-domain_calculations]).
If scaling was used then the difference of logs can be scaled, preventing
overflow of any log sample: sa - sb = s * (a - b).
When scaling is applied to tiny alpha the Dirichlet sampler supports all finite
non-zero alpha.
> 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)
