Improving Function Approximations using Newton–Raphson Iterations

Suppose we want to evaluate a function $y = f(x)$, which is either slow or not directly computable. If we also have (and can efficiently evaluate) all of:

• $f^{-1}(y)$, the function’s inverse
• $\left[f^{-1}\right]'(y)$, the derivative of the function’s inverse
• $f^\ast(x)$, an approximant to $f(x)$

then we can use the Newton–Raphson method to iteratively compute more accurate approximations.

Derivation

The Newton–Raphson method finds successively better approximations to a real function’s roots. By constructing a new function $g(y)$ with a root at the solution,

we can produce successive approximations $y^\ast_n$ to $f(x)$ by evaluating

Note that $x$ is treated as a constant when differentiating and evaluating.

Example: Square Root

// Evaluates a Newton-Raphson approximation to f(x) = sqrt(x)
float sqrtApproxNR(const float& x) {
    if (x < 0.0f) return NAN;
    
    // Avoid a division-by-zero later on
    if (x == 0.0f) return 0.0f;

    // Initial rough approximation to sqrt(x) using bit-level manipulation
    float y = std::bit_cast<float>((std::bit_cast<uint32_t>(x) + 0x3F800000) >> 1);
    
    // Compute NR step
    // f^-1(y) = y^2; [f^-1]'(y) = 2y
    // y = y - (y * y - x)/(2.0f * y);
    
    // After simplifying,
    y = 0.5f * (y + x/y);
    
    // Repeat as many times as desired; each iteration increases accuracy at
    //   the expense of execution time.
    y = 0.5f * (y + x/y);
    
    return y;
}

Function Transforms

Sometimes, the performance or stability can be improved by applying an additional transform $h(y)$ to both terms in $g(y)$:

The choice of $h(y)$ is not readily obvious, but, in general, should seek to identify functions that cancel terms in or expand the domain of $g(y)$.

Example: Wright Omega

Seeking to approximate the Wright omega function $f(x) = \omega(x)$ leads to

which can be problematic as $y \rightarrow 0$ due to limited floating point precision. By applying the transform $h(y) = \exp{y}$, we may instead calculate

which both removes the singularity at $y = 0$ and reduces the number of floating point operations performed each iteration.

Sources

• “Newton’s method,” Wikipedia, October 2024.
• “Wright omega function,” Wikipedia, October 2024.
• S. D’Angelo, L. Gabrielli, and L. Turchet, “Fast Approximation of the Lambert W Function for Virtual Analog Modelling,” Proceedings of the 22nd International Conference on Digital Audio Effects, Birmingham, UK, September 2019.