CODEWITHSUNDEEP
pythonnumerical-stability
Emil Eriksson
Emil Eriksson

Reputation: 2120

How can I make this function more numerically stable?

The following function is supposed to work similarly to pow(x, 1/k) but to be symmetric around the line y = 1 - x as well as not having a 0 or 1 slope at either end of [0, 1]:

def sym_gamma(x, k):
  if k == 1.0:
    return x
  a = 1.0 / k - 1.0
  b = 1.0 / a
  c = k + 1.0 / k - 2.0;
  return 1.0 / (a - c * x) - b

As can be seen, it is not defined when k = 1 so when that is the case, I simply return x. However, this special case handling is not enough since the function also behaves poorly when x is not equal to but very close to 1.0. For example sym_gamma(0.5, 1.00000001) yields 0.0 while it's supposed to return something very close to 0.5.

How can achieve the same thing without the poor stability? I know that I can introduce a tolerance with respect to k equaling 1.0 but it feels like a hack and I would also want to make sure that the function is perfectly smooth with regards to k.

Upvotes: 3

Views: 414

Answers (2)

ftorre
ftorre

Reputation: 606

Simplifying your expression seems to help with the precision. Numerical errors tends to accumulate in each operation. Thus, reducing the number of operation will reduce the chance of numerical errors.

We can notice that:

a = (1 - k) / k
b = k / (1 - k)
c = (1 - k) ** 2 / k
a - c * x = (1 - k) * (1 + x*k - x) / k
1.0 / (a - c * x) - b = x*k / (1 - x * (1 - k))

Then you can simply rewrite your method:

def sym_gamma(x, k):
  return x*k / (1 - x * (1 - k))

Instead of performing several division, only one division is computed. This method returns 0.5000000025 for sym_gamma(0.5, 1.00000001).

Upvotes: 5

Dev Bhuyan
Dev Bhuyan

Reputation: 620

As already answered by ftorre. Simplifying the equations into one would reduce the chance of numerical errors. This problem is, however, caused by the limited precision in the floating-point representation of the variables. Floating-point representation says about the way numbers are stored in actual memory. The default precision for float type variables in Python is double precision (meaning 8+8=16 values after the decimal). This can vary with the specific Python implementation, and when multiple operations are cascaded, the final result often comes out unexpected as above.

A more robust solution to your problem would be to use the decimal module, which lets you choose the precision. Here's an example of your code using the module:

from decimal import Decimal, getcontext

def sym_gamma(x, k):
    if k == 1.0:
        return x

    precision = 30
    getcontext().prec = precision

    a = Decimal(1.0) / Decimal(k) - Decimal(1.0)
    b = Decimal(1.0) / a
    c = Decimal(k) + Decimal(1.0) / Decimal(k) - Decimal(2.0)

    x = Decimal(x)

    return float(Decimal(1.0) / (a - c * x) - b)

This code also preserves your values of a, b, and c. You can set the precision to your desired value (I have preset it to 30)

Upvotes: 1

Related Questions