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.