log-sum-exp trick for signed/complex numbers

I need to evaluate a sum of values that are on many different orders of magnitude in scale but might be signed. I’ve had great luck with the “log-sum-exp” trick for an unsigned version of my problem, so I’m hoping to apply this to the signed values.

Suppose $$z,z’\in\mathbb C$$ are stored using separate real and imaginary parts in double precision. We wish to compute $$w:=\log(e^z + e^{z’})$$. Here I’m not concerned about which branch of log we get back—at the end of the day only $$e^w$$ matters, the logs are just for numerical stability.

For any $$a\in\mathbb C$$, by an identical argument to the log-sum-exp trick we know $$W=a+\log(e^{z-a}+ e^{z’-a}),$$ again up to branch of log.

Is there a good “policy” for choosing $$a$$ to promote numerical stability, similar to the max suggested for the real valued log-sum-exp trick?

• Tiny update: I guess the branch issue is irrelevant since I’m storing complex numbers in Cartesian coordinates. Commented Jan 23, 2020 at 0:32
• Have you thought about the Kahan summation algorithm? Commented Jan 23, 2020 at 1:47
• Unfortunately Kahan summation won't be enough -- my computation involves factorials so it needs to be on log scale. Commented Jan 23, 2020 at 3:34
• Just passing by and I am not an expert, but I would imagine that one needs to use log1p here to get some stability for inputs of wildly different size. Commented Jan 24, 2020 at 16:31
• Yes, probably a good idea for log-sum-exp generically as well! Commented Jan 27, 2020 at 16:36

I don't think the complex plane makes the log-sum-exp trick appreciably different, at least in Cartesian coordinates. In particular, if $$z=u+iv$$ then $$e^z=e^{u+iv}=e^u (\cos v + i\sin v).$$ Notice the $$v$$ part has magnitude 1 by construction, so overflow or underflow is principally caused by $$u$$.

Hence, a reasonable choice is to take $$a = \max( \mathrm{real}(z), \mathrm{real}(z') ).$$

• For numerical reasons it's better is to use the logp1 function. To that end choose $a = \mathrm{argmax}_{u\in\{z,z'\}} \mathrm{Re}[u]$. In this case, $w = a + \log \left(e^{z-a} + e^{z'-a}\right) = a + \log\left(1 + e^{z-z'}\right)$ if $\mathrm{Re}[z']>\mathrm{Re}[z]$ else $w=a + \log\left(1 + e^{z'-z}\right)$. Commented Sep 26, 2020 at 15:18