Conceptual vs Numerical

Posted on by John

One of the things that makes numerical computation interesting is that it often reverses the usual conceptual order of things, using advanced math to compute things that are introduced earlier.

Here’s an example I recently ran across [1]. The hyperbolic functions are defined in terms of the exponential function:

\begin{align*} \sinh(x) &= \frac{\exp(x) - \exp(-x)}{2} \\ \cosh(x) &= \frac{\exp(x) + \exp(-x)}{2} \\ \tanh(x) &= \frac{\exp(x) - \exp(-x)}{\exp(x) + \exp(-x)} \\ \end{align*}

But it may be more efficient to compute the exponential function in terms of hyperbolic functions. Specifically,

\exp(x) = \sinh(x) + \sqrt{\sinh(x)^2 + 1}

Why would you want to compute the simple thing on the left in terms of the more complicated thing on the right? Because hyperbolic sine is an odd function. This means that half its power series coefficients are zero: an odd function only has odd powers in its power series.

Suppose you need to compute the exponential function in an environment where you only have a limited number of registers to store constants. You can get more accuracy out of the same number of registers by using them to store coefficients in the power series for hyperbolic sine than for exp.

If we store n coefficients for sinh, we can include powers of x up to 2n − 1. And since the coefficient of x2n is zero, the error in our Taylor approximation is O(x2n+1). See this post for more on this trick.

If we stored n coefficients for exp, we could include powers of x up to n − 1, and our error would be O(xn).

To make things concrete, suppose n = 3 and x = 0.01. The error in the exp function would be on the order of 10−6, but the error in the sinh function would be on the order of 10−14. That is, we could almost compute exp to single precision and sinh to almost double precision.

(I’m glossing over a detail here. Would you really need to store, for example, the 1 coefficient in front of x in either series? For simplicity in the argument above I’ve implicitly said yes. Whether you’d actually need to store it depends on the details of your implementation.)

The error estimate above uses big-oh arguments. Let’s do an actual calculation to see what we get.

    from math import exp, sinh, sqrt
    
    def exp_taylor(x):
        return 1 + x + 0.5*x**2
    
    def sinh_taylor(x):
        return x + x**3/6 + x**5/120
    
    def exp_via_sinh(x):
        s = sinh_taylor(x)
        return s + sqrt(s*s + 1)
    
    def print_error(approx, exact):
        print(f"Computed: {approx}")
        print(f"Exact:    {exact}")
        print(f"Error:    {approx - exact}")
    
    x = 0.01
    approx1 = exp_taylor(x)
    approx2 = exp_via_sinh(x)
    exact   = exp(x)
    
    print("Computing exp directly:\n")
    print_error(approx1, exact)
    
    print()
    
    print("Computing exp via sinh:\n")
    print_error(approx2, exact)

This produces

    Computing exp directly:

    Computed: 1.0100500000000001
    Exact:    1.010050167084168
    Error:    -1.6708416783473012e-07

    Computing exp via sinh:

    Computed: 1.0100501670841682
    Exact:    1.010050167084168
    Error:    2.220446049250313e-16

Our errors are roughly what we expected from our big-oh arguments.

More numerical analysis posts

[1] I saw this identity in Matters Computational: Ideas, Algorithms, Source Code by Jörg Arndt.

7 thoughts on “Conceptual vs Numerical”

  1. How well does this work for other functions? If you take a function, re-write it in terms of it’s odd (or even) part, and compute values, is there generally a speedup?

  2. John

    It’s generally more efficient to apply power series expansions to either even or odd functions. But if you simply break a function into an even part and an odd part, there’s no savings. For example, if you write exp(x) as sinh(x) + cosh(x), you don’t save anything; even though sinh(x) and cosh(x) have half as many coefficients as exp(x), they don’t overlap, and so you end up storing all the coefficients of exp(x).

  4. John

    I’ve implicitly assumed that sqrt is available. In any case, sqrt is not computed via power series. It’s usually computed via Newton’s method or some variation on that method.

  5. Andrew Hunter

    Of course a simple Taylor expansion is often not the best series for a numerical computation, right? From my memory modern processors/libraries are more likely to use rational approximations, and even without that I’d expect that the minimal-error polynomial for sinh() over some interval to not be the Taylor polynomial, at least not in general.

    It’s not immediately clear to me if the even/odd trick here still works? Though I’d hope you’d still get something out of it.

  7. Severin Pappadeux

    Well, the best case, I would venture, would be computing sech(x) and expressing exp(x) via sech(x). Sech(x) has the advantage of odd and alternating sign series, so you could quickly estimate error, how many terms to compute etc

Comments are closed.