Operator notation - numerically stable forward recursions in linear-chain conditional random fields (CRFs) I would appreciate some clarity on my understanding of some notation for the mathematical specification of forward recursions for inference in linear-chain conditional random fields (CRFs), with a view to its implementation in the log-domain for numerical stability.
Query.
My query concerns the indexed operator denoted $\oplus_{i \in S}$ in equation (4.38) from the following excerpt from p329-330 of the Foundations and Trends in Machine Learning tutorial paper, "An Introduction to Conditional Random Fields" by Sutton, McCallum (2011):

A second approach to preventing underflow is to perform computations in the logarithmic domain e.g. the forward recursion (4.6) [which I have included for ease of reference]
$$\begin{equation*}
\alpha_t (j) = \sum_{i \in S} \Psi_t(j, i, x_t) \alpha_{t-1}(i)
\tag{4.6}
\end{equation*}$$
becomes
$$\begin{equation*}
\log \alpha_t (j) = \oplus_{i \in S} (\log \Psi_t(j, i, x_t) + \log \alpha_{t-1}(i))
\tag{4.38}
\end{equation*}$$
where $\oplus$ is the operator $a \oplus b = \log(e^a + e^b)$. At first this does not seem much of an improvement, since numerical precision is lost when computing $e^a$ and $e^b$. But $\oplus$ can be computed as
$$\begin{equation*}
a \oplus b = a + \log(1+e^{b-a}) = b + \log(1 + e^{a-b})
\tag{4.39}
\end{equation*}$$
which can be much more numerically stable if we pick the version of the identity with the smaller exponent.

Where $i$ ranges over the set of discrete states $S = \{1, ..., M\}$, $t$ indexes discrete-time so that equations (4.6) and (4.38) are recursions, and the idea is that we compute $\alpha_t(j)$ and $\log \alpha_t(j)$ for all $j \in S$ also.
I am comfortable with examples of operators in mathematics (e.g. product, summation, expectation etc.), but lack familiarity in the sense of never having studied their formal properties as mathematical objects of study. After doing some wiki-ing:
1. Is the operator $\oplus$ used in this particular context an example of the direct-sum operator in the sense of this wikipedia article? Or is it idiosyncratically chosen notation on the part of the author?
2. Which of the following expansions of the indexing $i \in S$ on this operator in equation (4.38) is correct?:
$$\begin{equation*}
\log \alpha_t(j) = (\log \Psi_t(j, 1, x_t) + \log \alpha_{t-1}(1)) \oplus \space  ... \space \oplus (\log \Psi_t(j, M, x_t) + \log \alpha_{t-1}(M))
\tag{A}
\end{equation*}$$
$$\begin{equation*}
\log \alpha_t(j) = (\log \Psi_t(j, 1, x_t) \oplus \log \alpha_{t-1}(1)) + \space  ... \space + (\log \Psi_t(j, M, x_t) \oplus \log \alpha_{t-1}(M))
\tag{B}
\end{equation*}$$
By analogy to the product and summation operators (which may not be appropriately justified), I would expect it to be $(A)$.
However, my confusion stems from wondering if it should be $(B)$, because the author has specifically chosen to present a numerically stable version of $a \oplus b$ in (4.39), known as the "LogSumExp" trick in machine learning, over two terms ($a$ and $b$), rather than over $M$ terms. Adding to my confusion somewhat is the fact that the clique potentials $\Psi_t$ are log-linear functions, with the form $e^{(\cdot)}$. Perhaps this is a case of overthinking leading me astray?
Further context.
The forward recursion described in equation (4.6), amongst other purposes, provides a means of efficiently computing the normalisation constant $Z(\mathbf{x}, \Theta)$ in conditional random fields (so that the computational complexity is not exponential). In the language of undirected directed graphical models, for state vectors $\mathbf{y} = (y_0, y_1, ... , y_T)$ and observation vectors $\mathbf{x} = (x_1, ..., x_T)$, CRFs specify the conditional probability $p(\mathbf{y} | \mathbf{x}; \Theta)$ as a product of clique potentials, or factors $\Psi_t(y_t, y_{t-1}, x_t)$, paramterised by $\Theta$:
$$\begin{align}
p(\mathbf{y} | \mathbf{x}; \Theta) &= \frac{1}{Z(\mathbf{x}, \Theta)} \prod^T_{t=1} \Psi_t(y_t, y_{t-1}, x_t)\end{align}$$
Where the clique potentials are also a weighted linear combination of feature functions $f_k(y_t, y_{t-1}, x_t)$, with each weight being an element $\theta_k$ of $\Theta = (\theta_1, ..., \theta_K)$:
$$\Psi_t(y_t, y_{t-1}, x_t) = \exp \left(\sum^K_{k=1} \theta_k f_k(y_t, y_{t-1}, x_t) \right)$$
And the normalisation constant $Z(\mathbf{x}, \Theta)$ is given by:
$$Z(\mathbf{x}, \Theta) = \sum_{y_1} ... \sum_{y_T} \prod^T_{t=1} \Psi_t(y_t, y_{t-1}, x_t)$$
 A: The correct expansion of $\oplus_{i \in S}(\log \Psi_t(j, i, x_t) + \log \alpha_{t-1}(i))$ you are looking for is $(A)$.
Taking the logs of the forward recursion in equation (4.6), we have:
$$\begin{align}
\log \alpha_t(j) &= \log \sum_{i \in S} \Psi_t(j, i, x_t) \alpha_{t-1}(i) \\
&= \log \left[\Psi_t(j, 1, x_t) \alpha_{t-1}(1) + ... +  \Psi_t(j, M, x_t) \alpha_{t-1}(M)\right] \\
&= \log \left[e^{\log \Psi_t(j, 1, x_t)} e^{\log \alpha_{t-1}(1)} + ... e^{\log \Psi_t(j, M, x_t)}e^{\log \alpha_{t-1}(M)} \right] \\
&= \log \left[e^{\log \Psi_t(j, 1, x_t) + \log \alpha_{t-1}(1)} + ... e^{\log \Psi_t(j, M, x_t) + \log\alpha_{t-1}(M)} \right] \\
&= \oplus_{i \in S}(\log \Psi_t(j, i, x_t) + \log \alpha_{t-1}(i))
\end{align}$$
Now, using the definition of the operator $\oplus_{i \in S} \space u_i = \log[e^{u_1} +  ... + e^{u_M}]$ where $u_i = \log \Psi_t(j, i, x_t) + \log \alpha_{t-1}(i)$ to compute $\log \alpha_t(j)$ will run into numerical precision errors. Furthermore, notice that $\oplus_{i \in S} u_i$ is of the form $\log \sum_i \exp(u_i)$.
To compute a numerically stable version of this, we use the Log-Sum-Exp trick and instead compute:
$$\begin{align}
\log \alpha_t(j) &= \oplus_{i \in S} \space u_i \\
&= \max \{u_1, ..., u_M \} + \log \left(\sum^M_{i=1} \exp(u_i - \max \{u_1, ..., u_M \}) \right) 
\end{align}$$
Which is the same as $a \oplus b$ generalised to $M$ terms rather than two terms.
