3
$\begingroup$

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?

$\endgroup$
5
  • $\begingroup$ Tiny update: I guess the branch issue is irrelevant since I’m storing complex numbers in Cartesian coordinates. $\endgroup$ Jan 23 '20 at 0:32
  • $\begingroup$ Have you thought about the Kahan summation algorithm? $\endgroup$ Jan 23 '20 at 1:47
  • $\begingroup$ Unfortunately Kahan summation won't be enough -- my computation involves factorials so it needs to be on log scale. $\endgroup$ Jan 23 '20 at 3:34
  • 1
    $\begingroup$ 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. $\endgroup$ Jan 24 '20 at 16:31
  • $\begingroup$ Yes, probably a good idea for log-sum-exp generically as well! $\endgroup$ Jan 27 '20 at 16:36
3
$\begingroup$

After thinking about this some more, I can answer this one myself!

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') ).$

$\endgroup$
1
  • $\begingroup$ 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)$. $\endgroup$ Sep 26 '20 at 15:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.