4
$\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?

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

1 Answer 1

Reset to default
4
$\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') ).$

Improve this answer
$\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$
    – Josh Albert
    Commented Sep 26, 2020 at 15:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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