What is the best way to apply the log-sum-exp trick in this situation? I am aware of the "log-sum-exp" trick for calculating the logarithm of sums that handles overflow and underflow issues. However, I would like to know more about how it works. In particular, I am wondering: what is the most numerically precise way to evaluate the following expression using the "log-sum-exp" trick?
The expression I am concerned with this the following:  
$$
\log\left[\sum_i \left(\frac{a^i }{\sum_j a^j}\right) \left(\frac{b^i }{\sum_k b^k} \right)\right],
$$
where $i$ and $j$ are indices, not powers.
Instead of having $\{a^i\}$ and $\{b^i\}$, I have their logarithms. As you can see, there are several ways to apply this rule-of-thumb. Here are a few.
 1. The all-together way
$$
m + \log\left[\sum_i \exp\left\{ \log a^i + \log b^i - \log\sum_ja^j - \log\sum_k b^k - m \right\}\right]
$$
where $m = \max_i\{ \log a^i + \log b^i - \log\sum_ja^j - \log\sum_k b^k  \}$ and we could apply the log-sum-exp technique to the inner sums as well, which would be used for both calculating the maximum $m_1$, and for calculating the outer sum over $i$.
 2. The split-up way:
You could also split up the thing into products like $\log\sum_i a^ib^i - \log \sum_j a^j - \log \sum_k b^k$, and use the trick on each of the three terms:
\begin{align*}
&m_1 + \log\sum_i \exp [\log a^i + \log b^i - m_1] \\
&- m_2 - \log\sum_j\exp[\log a^j - m_2] \\
&- m_3 -\log \sum_k\exp[ \log b^k - m_3].
\end{align*}
 3. The straightforward way:
Don't do any of this. Just exponentiate everything, get $\{a^i\}$ and $\{b^i\}$, and compute the original expression.
 A: So, the sum terms in the denominator are invariant, over each of the terms in the outer sum, since you're always summing over all $j$ in the denominator sums.
Then, what you want to do probably is make sure that the largest of the terms in the outer sum is reasonably representable, e.g. if it equals $1$, that probably meets our requirements.
The numerator is not combinatorial in the number of terms, but is linear in values of $i$, i.e. it's:
$$
\exp(a_1)\exp(b_1) + \exp(a_2)\exp(b_2) + \dots \\
= \exp(a_1 + b_1) + \exp(a_2 + b_2) + \dots
$$
(where I'm writing $\exp(a_i)$ to mean your $a^i$, since I want to make the exponentials explicit).
So, to make sure the outer sum has $1$ as its maximum numerator, we can determine $\max_i(a_i +  b_i)$, and divide by the exponential of that.
Let:
$$
E = \max_i(a_i + b_i)
$$
Then we'd calculate:
$$
\log \left(
    \sum_i \left(
          \frac{\exp(a_i)\exp(b_i)/\exp(E)}
             {\left(\sum_j \exp(a_j) \sum_j(\exp(b_j))\right) / \exp(E)}
    \right)
\right)
$$
$$
= \log \left(
    \sum_i \left(
          \frac{\exp(a_i + b_i - \max_i(a_i + b_i))}
             {\left(\sum_j \exp(a_j) \sum_j(\exp(b_j))\right)/\exp(\max_i(a_i + b_i))}
    \right)
\right)
$$
What to do with the $\exp(\max_i(a_i + b_i))$ term in the denominator?
So, let's determine:
$$
i' = \arg \max_i(a_i + b_i)
$$
Then $E = a_{i'} + b_{i'}$. And our expression becomes:
$$
= \log \left(
    \sum_i \left(
          \frac{\exp(a_i + b_i - \max_i(a_i + b_i))}
             {\left(\sum_j \exp(a_j) \sum_j(\exp(b_j))\right)/\exp(a_{i'} + b_{i'})}
    \right)
\right)
$$
$$
= \log \left(
    \sum_i \left(
          \frac{\exp(a_i + b_i - \max_i(a_i + b_i))}
             {\left(\sum_j \exp(a_j) \sum_j(\exp(b_j))\right)\exp(-a_{i'})\exp(-b_{i'})}
    \right)
\right)
$$
and then we can just multiply the exponential sum of $a_j$ in the denominator by $\exp(-a_{i'})$, and the same for the exponential sum of $b_j$. So we have:
$$
\log \left(
    \sum_i \left(
          \frac{\exp(a_i + b_i - \max_i(a_i + b_i))}
             {\left(\sum_j \exp(a_j - a_{i'}) \sum_j(\exp(b_j - b_{i'}))\right)}
    \right)
\right)
$$
and since we've determined $i'$ anyway, we could if we want rewrite the numerator slightly:
$$
= \log \left(
    \sum_i \left(
          \frac{\exp(a_i - a_{i'} + b_i - b_{i'}))}
             {\left(\sum_j \exp(a_j - a_{i'}) \sum_j(\exp(b_j - b_{i'}))\right)}
    \right)
\right)
$$
This expression will give small numerical error, since the largest numerator will be exactly 1, neither under nor overflowing. Other numerators will all be smaller than this, and not overflow.
