We can implement the above computation by exponentiating each number, then
summing them, then taking a logarithm as follows (written in
Julia).

logsumexp_naive(X)=log(sum(exp(X)))

When the above function returns a finite number then it is numerically
accurate. However, the above computation is not robust if one of the elements is very large (say, larger than 710 for double precision IEEE floating point).
Then \(\exp(x_i)\) returns a floating point infinity and the entire computation
returns a floating point infinity as well.

Standard Batch Solution

The standard solution to this problem is to use the mathematical identity

which holds for any \(\alpha \in \mathbb{R}\).
By selecting \(\alpha = \max_{i=1,\dots,n} x_i\) no argument to the
\(\exp\)-function will be larger than zero and the above naive computation can
be applied on the transformed numbers.
The code is as follows.

Code such as the above is used in almost all packages for performing
statistical computation and is described as the standard solution, see e.g.
here
and here.

However, there are the following problems:

It requires two scans over the data array, one to find the maximum, one to
compute the summation. For modern systems and large input arrays the above
computation is memory-bandwidth limited so two memory scans mean twice the
runtime.

It requires knowledge of the number of elements in the sum prior to
computation.

Streaming log-sum-exp Computation

The solution is to also compute the maximum element in a streaming manner and
to correct a running estimate whenever a new maximum is found.
I have not seen this solution elsewhere, but I hope you may find it useful.