parallelizing log-sum-exp

I have some approximate likelihoods: $$L_1, \ldots, L_n$$. Each is quite expensive to calculate. They're approximate because they use random numbers. Each of them is being calculated on the same data set of size $$N$$. Actually, what's really being calculated are the log-likelihoods: $$\ell_i$$ for $$i=1,\ldots,n$$, not the likelihoods.

I'd like to average together these likelihoods. If each likelihood is unbiased, then so is the sample mean. This is an important requirement when using the Pseudo-marginal Metropolis-Hastings sampler, for instance.

To avoid numerical underflow, from exponentiating numbers close to $$-\infty$$, I'd like to use the log-sum-exp trick:

$$\log \left( \frac{1}{n}\sum_{i} L_i\right) = m + \log \sum_{i}\exp[\ell_i - m] - \log n$$

Question

Typically $$m$$ is set to $$\max(\ell_1, \ldots, \ell_n)$$. However,

• if one is using multithreading, calculating the maximum will require more thread communication or memory or locking, and
• I don't want to waste time calculating some example/pilot likelihoods (remember, they're expensive to compute).

Is there a good rule of thumb to use as a substitute? Perhaps $$m := -.5N$$? I reason that, if the data are iid, and each data point likelihood is approx. $$r \approx .6$$, then $$L_i \approx r^N$$ and $$\ell_i \approx N \log r \approx -.5 N$$. Another piece of relevant information: I expect these $$\ell_i$$ to be quite left-skewed.

The perfect answer is probably problem-specific, but I was wondering if anyone had some nice anecdotes or resources.

• There's no way to compute the sum without inter-thread communication or shared memory of the same sort as computing the per-thread max, so it's not clear based on the notation how avoiding inter-thread communication to find the maximum is an optimization in the end -- you having to go from N pieces of data for N threads to one number in the end. What's a typical value for N? Jul 31, 2022 at 17:03
• @jwimberley right the question might be too narrow here. Every time I update the running sum/average (anything that is the result of a commutative op) I have to hold a lock/mutex. I'd prefer to minimize doing this because this is what slows things down. Another logical solution would be to, instead of sharing an average or sum, to share a container of all the likelihoods, well then the maximum would be available like you say. I'm considering that as well. However, that would require more extensive refactoring at this point in time and I have a deadline :) Jul 31, 2022 at 17:08
• $\max$ is a commutative operation as well, so maybe I'll hold on to the max in addition to the sum, but no that won't work because you can't start summing together until you have the max, which requires two passes over the entire collection of log-likelihoods Jul 31, 2022 at 17:10
• If you're using C++ you can probably calculate both the max and then the sum with compare exchange operations on atomics to avoid some mutex overhead (e.g. stackoverflow.com/questions/16190078/…) Jul 31, 2022 at 18:07
• Your requirements are unclear. For example, why is not not effective to compute the $l_i$ and in a second pass perform the averaging? It's also possible that undisclosed information about your specific problem could solve it. For instance, if the random variation is not great, so that the $l_i$ are close to one another, then it wouldn't be necessary to find the maximum: you could use any of the values of the $l_i$ as the offset $m$ in the calculation.
– whuber
Aug 1, 2022 at 12:47

• uses a strict lock (i.e. a std::unique_lock<std::shared_mutex>)whenever writing individual approximations to the shared std::vector (e.g. here)
• uses a strict lock (i.e. a std::unique_lock<std::shared_mutex>) whenever writing the new input parameters to the log-likelihood function or refreshing the sum to begin accumulating again (e.g. here)
Because the bottleneck is the sampling-based approximations---not the intra-thread communication required for summing and maximum calculations---the reading from shared memory is only loosely protected using a shared lock (i.e. a std::shared_lock<std::shared_mutex>).
When all the work is done, the the standard log-sum exp trick is used for the final calculation. This uses the maximum $$m$$. See here.
If you wanted to visualize the whole process as a graph: first, one thread will write the new parameter input for the likelihood. Second, in parallel all the workers will calculate their approximation to the likelihood. When each thread is finished with a likelihood approximation, they grab the mutex and write to the output std::vector. Finally, when all threads are done, only one thread processes all these numbers with log-sum-exp and spits out the one floating point.