Finding good approximation for $x^{1/2.4}$

Question

I would like to a good (8 bits accuracy) approximation for $x^{1/2.4}$ in the range $[0, 1]$. This transform is used for converting linear intensities to SRGB compressed values, so it's important that I make it run fast.

Plot of function:

Using a simple polynomial isn't practical because

1. the function has lots of high order derivatives
2. the function is roughly asymptotic to $x$, which is very different from the asymptotic behavior of high order polynomials

I've already have code that constructs an arbitrary degree polynomial for any function by minimizing the square error and even for a 10th degree polynomial, the accuracy is still only like 6 significant bits.

Then I learned about rational function approximations, which will have a much better asymptotic behavior. But the problem is I don't know how to find the optimal coefficients. There's the Pade formulation which creates an approximation around a single point, but since it doesn't use global information, it can be a very bad fit overall just like Taylor series.

I had Mathematica create an approximation of the form $(a_0 + a_1 x + a_2 x^2) / (1 + a_3 x)$ with PadeApproximant[x^(1/2.4), {x, 0.2, {3, 2}}], which is much better than a simple 3rd degree polynomial, but still not good enough, so I want to find a globally optimal solution, probably of the same form.

I tried finding a least squares solution like before, but it involves 4 huge, non-linear polynomial equations, that is taking Mathematica forever (I've waited 1/2 hour so far) to solve.

Can anyone suggest how to solve those non-linear equations, or another way to find a rational function approximation, or an entirely different approximation?

Thanks for any help

Answer 1

I think your reasoning is correct -- a rational approximation is probably the best solution. But, constructing rational approximations is difficult.

You have to decide how you will measure error. If you just look at maximum difference between the value of the function and it's approximation, then you are doing "uniform" or "minimax" approximation, which is usually more difficult than least-squares approximation. As you mentioned, Pade approximation is not very good for this, because it's a Taylor-series-like approach that's fixated on a single point.

I would recommend that you try the Mathematica RationalInterpolation function, if you haven't already.

Or, if you have access to Matlab, there is an add-on called Chebfun that does a very good job of constructing minimax polynomial and rational approximations. There are commands named ratinterp and remez, and a couple of others. The name "remez" comes from the Remez Exchange Algorithm, which is the standard way of computing minimax approximations. This section of the Chebfun documentation explains the techniques that are available. Section 4.8 covers rational approximations.

There is a nice book by Trefethen that provides a modern computationally-oriented account of approximation theory, with a focus on Chebfun. Insightful, and easy to read.

Edit

I tried the following Mathematica code

gx = GeneralMiniMaxApproximation[{t, 1 - t + t^(1/2.4)}, 
                             {t, {0, 1}, 3, 3}, x, 
                              MaxIterations -> 200, WorkingPrecision -> 50]

and then

hx = gx - 1 + x

This gave an approximation $hx$ of your function

$$ \frac{x^4 - 1.3730926529x^3 - 1.2231795079x^2 - 0.0120332034x - 3.892203211 \times 10^{-7}} {x^3 - 2.47604929086x^2 - 0.13996763348x - 0.00008072719}$$

The error function looks like this:

Not quite 8 bits, but getting close. You can fiddle with the degrees of the numerator and denominator to try to get something better.

Answer 2

One observes that $a^{\frac{1}{2.4}} = a^{\frac{5}{12}} = a \cdot a^{-\frac{1}{3}} \cdot a^{-\frac{1}{4}} = a^{\frac{1}{2}} \cdot a^{-\frac{1}{12}}$. There is a division-free iteration of the Householder type with cubic convergence for computing the reciprocal of the $n$-th root, $a^{-\frac{1}{n}}$, as follows: $$h_{i} \colon= 1-x_{i}^{n}a$$ $$x_{i+1} \colon= x_{i} + x_{i} h_{i}\left(\frac{1}{n}+\frac{1+n}{2n^{2}}h_{i}\right)$$ where $\frac{1}{n}$ and $\frac{1+n}{2n^2}$can obviously be pre-computed as constants. On platforms with a reasonably fast square root, the latter approximation is likely preferable. On the other hand, many high-end processors provide sufficient instruction-level parallelism to support concurrent computation of reciprocal cube root and reciprocal quad root.

One way of generating starting approximations $x_0$ is to use bit-manipulation of the IEEE-754 single-precision floating-point data, as shown in this question on Stackoverflow. Since the accuracy of the starting approximations is around $3.5 \cdot 10^{-2}$, one iteration for each of the reciprocal roots suffices to achieve the desired accuracy in the final result.

I performed an exhaustive test on the interval $[2^{-126}, 2^{126}]$, although sRGB conversion mentioned in the question only requires the computation of $a^{\frac{5}{12}}$ on the much narrower interval $[0.0031808, 1]$.

Implementation of the exemplary ISO-C99 code below using the SIMD instruction set of an x86 processor (this was indicated by asker in a comment) would appear to require at least SSE 4.1, as a 32-bit integer multiply is required.

#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <math.h>

#define USE_SQRT  (1)

float root_m3 (float a);  // compute a**(1/(-3))
float root_m4 (float a);  // compute a**(1/(-4))
float root_m12 (float a); // compute a**(1/(-12))

// Compute a**(5/12) on [2**(-126), 2**126]
// Maximum relative error = 1.4921e-3 (USE_SQRT==1); 4.2707e-4 (USE_SQRT==0)
float pow_5over12 (float a)
{
#if USE_SQRT
    return sqrtf (a) * root_m12 (a);
#else // USE_SQRT
    return a * root_m3 (a) * root_m4 (a);
#endif // USE_SQRT
}

// Compute a**3
float pow_3 (float a) 
{ 
    return a * a * a; 
}

// Compute a**4
float pow_4 (float a) 
{ 
    a = a * a; 
    return a * a; 
}

// Compute a**12
float pow_12 (float a) 
{ 
    a = a * a * a; 
    a = a * a; 
    return a * a; 
}

// Reinterpret bits of IEEE-754 'binary32' as 32-bit unsigned integer
uint32_t float_as_uint32 (float a) {
    uint32_t r;
    memcpy (&r, &a, sizeof r); 
    return r;
}

// Reinterpret bits of 32-bit unsigned integer as IEEE-754 'binary32'
float uint32_as_float (uint32_t a) {
    float r; 
    memcpy (&r, &a, sizeof r); 
    return r;
}

// Compute a**(1/(-3)). Maximum relative error = 2.04e-4
float root_m3 (float a)
{
    float r, h;
    // initial approximation
    r = uint32_as_float (0x54a1da3a - (float_as_uint32 (a) >> 16) * 0x5555);
    // Householder iteration with cubic convergence
    h = 1.0f - pow_3 (r) * a;
    r = (1.0f/3.0f + 4.0f/18.0f * h) * h * r + r;
    return r;
}

// Compute a**(1/(-4)). Maximum relative error = 2.25e-4
float root_m4 (float a)
{
    float r, h;
    // initial approximation
    r = uint32_as_float (0x4f583c27 - (float_as_uint32 (a) / 4));
    // Householder iteration with cubic convergence
    h = 1.0f - pow_4 (r) * a;
    r = (1.0f/4.0f + 5.0f/32.0f * h) * h * r + r;
    return r;
}

// Compute a**(1/(-12)). Maximum relative error = 1.43e-3
float root_m12 (float a)
{
    float r, h;
    // initial approximation
    r = uint32_as_float (0x44c39f16 - (float_as_uint32 (a) >> 16) * 0x1555);
    // Householder iteration with cubic convergence
    h = 1.0f - pow_12 (r) * a;
    r = (1.0f/12.0f + 13.0f/288.0f * h) * h * r + r;
    return r;
}

// Test approximation on interval [2**(-126), 2**126]
int main (void)
{
    float a, res, start, stop;
    double ref, relerr, maxrelerr;

    start = exp2f (-126.0f);
    stop = exp2f (126.0f);

    maxrelerr = 0;
    a = start;
    while (a <= stop) {
        res = pow_5over12 (a);
        ref = pow ((double)a, 5.0/12.0); 
        relerr = fabs ((res - ref) / ref);
        if (relerr > maxrelerr) maxrelerr = relerr;
        a = uint32_as_float (float_as_uint32 (a) + 1);
    }
    printf ("maxrelerr = %11.4e\n", maxrelerr);
    return EXIT_SUCCESS;
}
```

Answer 3

Consider $\frac{1}{2.4} \approx 0.41$, perhaps starting with a square root and approximating the ratio is easier. Or just use approximations for the logarithm and exponential. Those functions are directly implemented in hardware nowadays, and quite cheap. Or select judiciously some points and interpolate. Or approximate through splines. The magic fast inverse square root algorithm (really, its justification) might also give some ideas.

Just make sure this operation is really relevant performance wise before going down this path.

Answer 4

See also stackoverflow optimizations-for-pow-with-const-non-integer-exponent for other approaches to the same problem.

If you look at your graph of $x^{5/12}$ its obviously got a different shape near $x=0$, than near $x=1$. This is a good hint that simple polynomials wont fit so well over a domain including both $0$ and $1$.

So I would consider a strength reduction:

consider that in floating point representation $x = m \cdot 2^k$ where $0.5 \le m < 1$ so $$x^{5/12} = m^{5/12} \cdot 2^{5k/12}$$

MiniMaxApproximation[x^(5/12), {x, {1/2, 1}, 2, 2}]

Offers a rational O(2,2) polynomnial with a relative error less than $5.7\times10^{-7}$ for $m$ on $\left[\frac12,1\right]$.

For the second part of the strength reduction you can either use $Power\left[2,\frac {5k}{12}\right]$ or, since $k$ is an integer, optimise it with a 12 element lookup table and some modulo arithmetic.

Noting that, on most machines, division is very slow (and not good for numerical accuracy) I also considered simple polynomials and

MiniMaxApproximation[x^(5/12), {x, {1/2, 1}, 7, 0}]

Offers an O(7) polynomial with a relative error less than $2.4\times10^{-8}$ or better than half an ULP in single precision. Which should be about the same speed as the O(2,2). For completeness, here it is in Horner Form

0.246873 + 
   x (F954 + 
     x (-2.71714 + 
        x (3.86663 + 
          x (-3.85026 + x (2.47365 + (-0.919353 + 0.15006 x) x)))))

Answer 5

$$x^\frac{1}{2.4} =x^\frac{5}{12} = (x^\frac{5}{4})^\frac{1}{3} $$

Thus, assuming you have access to square root, you can calculate $x^\frac{5}{4}$ precisely, and thus you only need a good approximation to cube root in the range 0 to 1.

As noted, there are multiple ways one can do such an approximation. For example, use Newton-Raphson to quickly converge to a solution; one can get this to run faster by pre-calculating some values at intervals and caching them, so one can converge faster to the result by looking up the nearest starting point.

A simple approximation though can be achieved by taking $\sqrt[3]{x} \approx \frac{x^\frac{10}{32} + 2x^\frac{11}{32}}{3}$ which gives better than 10 bit accuracy for the cube root for numbers in this range, and only relies on having a square root function available

Answer 6

Continuing on from the simple approximation method I provided in my other answer, it appears that this method works quite well for doing such approximations.

Thus, for example, noting that $\frac{5}{12} = \frac{6.6666...}{16}$, gives us the following approximation $x^\frac{5}{12} \approx \frac{x^\frac{6}{16} + 2x^\frac{7}{16}}{3}$, which gives 9-bit precision over the desired range, and can be calculated quickly given square root.