单调递增函数的反函数，log10() 的 OverflowError答案

【问题标题】：Inverse function for monotonically increasing function, OverflowError for log10()单调递增函数的反函数，log10() 的 OverflowError
【发布时间】：2012-06-28 04:00:53
【问题描述】：

对于赋值，我们被要求创建一个返回反函数的函数。基本问题是从平方函数创建平方根函数。我想出了一个使用二分搜索的解决方案和另一个使用牛顿方法的解决方案。我的解决方案似乎适用于立方根和平方根，但不适用于 log10。以下是我的解决方案：

#Binary Search
def inverse1(f, delta=1e-8):
    """Given a function y = f(x) that is a monotonically increasing function on
    non-negative numbers, return the function x = f_1(y) that is an approximate
    inverse, picking the closest value to the inverse, within delta."""
    def f_1(y):
        low, high = 0, float(y)
        last, mid = 0, high/2
        while abs(mid-last) > delta:
            if f(mid) < y:
                low = mid
            else:
                high = mid
            last, mid = mid, (low + high)/2
        return mid
    return f_1

#Newton's Method
def inverse(f, delta=1e-5):
    """Given a function y = f(x) that is a monotonically increasing function on
    non-negative numbers, return the function x = f_1(y) that is an approximate
    inverse, picking the closest value to the inverse, within delta."""
    def derivative(func): return lambda y: (func(y+delta) - func(y)) / delta
    def root(y): return lambda x: f(x) - y
    def newton(y, iters=15):
        guess = float(y)/2
        rootfunc = root(y)
        derifunc = derivative(rootfunc)
        for _ in range(iters):
            guess = guess - (rootfunc(guess)/derifunc(guess))
        return guess
    return newton

无论使用哪种方法，当我在教授的测试函数中输入 n = 10000 for log10() 时，我都会收到此错误：（例外：使用我的牛顿方法函数时，log10() 很远，而这种二分搜索方法在达到输入阈值之前是相对准确的，但无论哪种方式，当 n = 10000 时，两种解决方案都会抛出此错误）

   2: sqrt =     1.4142136 (    1.4142136 actual); 0.0000 diff; ok
   2: log =     0.3010300 (    0.3010300 actual); 0.0000 diff; ok
   2: cbrt =     1.2599211 (    1.2599210 actual); 0.0000 diff; ok
   4: sqrt =     2.0000000 (    2.0000000 actual); 0.0000 diff; ok
   4: log =     0.6020600 (    0.6020600 actual); 0.0000 diff; ok
   4: cbrt =     1.5874011 (    1.5874011 actual); 0.0000 diff; ok
   6: sqrt =     2.4494897 (    2.4494897 actual); 0.0000 diff; ok
   6: log =     0.7781513 (    0.7781513 actual); 0.0000 diff; ok
   6: cbrt =     1.8171206 (    1.8171206 actual); 0.0000 diff; ok
   8: sqrt =     2.8284271 (    2.8284271 actual); 0.0000 diff; ok
   8: log =     0.9030900 (    0.9030900 actual); 0.0000 diff; ok
   8: cbrt =     2.0000000 (    2.0000000 actual); 0.0000 diff; ok
  10: sqrt =     3.1622777 (    3.1622777 actual); 0.0000 diff; ok
  10: log =     1.0000000 (    1.0000000 actual); 0.0000 diff; ok
  10: cbrt =     2.1544347 (    2.1544347 actual); 0.0000 diff; ok
  99: sqrt =     9.9498744 (    9.9498744 actual); 0.0000 diff; ok
  99: log =     1.9956352 (    1.9956352 actual); 0.0000 diff; ok
  99: cbrt =     4.6260650 (    4.6260650 actual); 0.0000 diff; ok
 100: sqrt =    10.0000000 (   10.0000000 actual); 0.0000 diff; ok
 100: log =     2.0000000 (    2.0000000 actual); 0.0000 diff; ok
 100: cbrt =     4.6415888 (    4.6415888 actual); 0.0000 diff; ok
 101: sqrt =    10.0498756 (   10.0498756 actual); 0.0000 diff; ok
 101: log =     2.0043214 (    2.0043214 actual); 0.0000 diff; ok
 101: cbrt =     4.6570095 (    4.6570095 actual); 0.0000 diff; ok
1000: sqrt =    31.6227766 (   31.6227766 actual); 0.0000 diff; ok
Traceback (most recent call last):
  File "/CS212/Unit3HW.py", line 296, in <module>
    print test()
  File "/CS212/Unit3HW.py", line 286, in test
    test1(n, 'log', log10(n), math.log10(n))
  File "/CS212/Unit3HW.py", line 237, in f_1
    if f(mid) < y:
  File "/CS212/Unit3HW.py", line 270, in power10
    def power10(x): return 10**x
OverflowError: (34, 'Result too large')

这里是测试函数：

def test():
    import math
    nums = [2,4,6,8,10,99,100,101,1000,10000, 20000, 40000, 100000000]
    for n in nums:
        test1(n, 'sqrt', sqrt(n), math.sqrt(n))
        test1(n, 'log', log10(n), math.log10(n))
        test1(n, 'cbrt', cbrt(n), n**(1./3.))


def test1(n, name, value, expected):
    diff = abs(value - expected)
    print '%6g: %s = %13.7f (%13.7f actual); %.4f diff; %s' %(
        n, name, value, expected, diff,
        ('ok' if diff < .002 else '**** BAD ****'))

这些是测试的设置方式：

#Using inverse() or inverse1() depending on desired method
def power10(x): return 10**x
def square(x): return x*x
log10 = inverse(power10)
def cube(x): return x*x*x
sqrt = inverse(square)
cbrt = inverse(cube)
print test()

发布的其他解决方案似乎在运行全套测试输入时没有问题（我试图不查看发布的解决方案）。对此错误有任何见解吗？

似乎共识是数字的大小，但是，我教授的代码似乎适用于所有情况：

#Prof's code:
def inverse2(f, delta=1/1024.):
    def f_1(y):
        lo, hi = find_bounds(f, y)
        return binary_search(f, y, lo, hi, delta)
    return f_1

def find_bounds(f, y):
    x = 1
    while f(x) < y:
        x = x * 2
    lo = 0 if (x ==1) else x/2
    return lo, x

def binary_search(f, y, lo, hi, delta):
    while lo <= hi:
        x = (lo + hi) / 2
        if f(x) < y:
            lo = x + delta
        elif f(x) > y:
            hi = x - delta
        else:
            return x;
    return hi if (f(hi) - y < y - f(lo)) else lo

log10 = inverse2(power10)
sqrt = inverse2(square)
cbrt = inverse2(cube)

print test()

结果：

     2: sqrt =     1.4134903 (    1.4142136 actual); 0.0007 diff; ok
     2: log =     0.3000984 (    0.3010300 actual); 0.0009 diff; ok
     2: cbrt =     1.2590427 (    1.2599210 actual); 0.0009 diff; ok
     4: sqrt =     2.0009756 (    2.0000000 actual); 0.0010 diff; ok
     4: log =     0.6011734 (    0.6020600 actual); 0.0009 diff; ok
     4: cbrt =     1.5865107 (    1.5874011 actual); 0.0009 diff; ok
     6: sqrt =     2.4486818 (    2.4494897 actual); 0.0008 diff; ok
     6: log =     0.7790794 (    0.7781513 actual); 0.0009 diff; ok
     6: cbrt =     1.8162270 (    1.8171206 actual); 0.0009 diff; ok
     8: sqrt =     2.8289337 (    2.8284271 actual); 0.0005 diff; ok
     8: log =     0.9022484 (    0.9030900 actual); 0.0008 diff; ok
     8: cbrt =     2.0009756 (    2.0000000 actual); 0.0010 diff; ok
    10: sqrt =     3.1632442 (    3.1622777 actual); 0.0010 diff; ok
    10: log =     1.0009756 (    1.0000000 actual); 0.0010 diff; ok
    10: cbrt =     2.1534719 (    2.1544347 actual); 0.0010 diff; ok
    99: sqrt =     9.9506714 (    9.9498744 actual); 0.0008 diff; ok
    99: log =     1.9951124 (    1.9956352 actual); 0.0005 diff; ok
    99: cbrt =     4.6253061 (    4.6260650 actual); 0.0008 diff; ok
   100: sqrt =    10.0004883 (   10.0000000 actual); 0.0005 diff; ok
   100: log =     2.0009756 (    2.0000000 actual); 0.0010 diff; ok
   100: cbrt =     4.6409388 (    4.6415888 actual); 0.0007 diff; ok
   101: sqrt =    10.0493288 (   10.0498756 actual); 0.0005 diff; ok
   101: log =     2.0048876 (    2.0043214 actual); 0.0006 diff; ok
   101: cbrt =     4.6575475 (    4.6570095 actual); 0.0005 diff; ok
  1000: sqrt =    31.6220242 (   31.6227766 actual); 0.0008 diff; ok
  1000: log =     3.0000000 (    3.0000000 actual); 0.0000 diff; ok
  1000: cbrt =    10.0004883 (   10.0000000 actual); 0.0005 diff; ok
 10000: sqrt =    99.9991455 (  100.0000000 actual); 0.0009 diff; ok
 10000: log =     4.0009756 (    4.0000000 actual); 0.0010 diff; ok
 10000: cbrt =    21.5436456 (   21.5443469 actual); 0.0007 diff; ok
 20000: sqrt =   141.4220798 (  141.4213562 actual); 0.0007 diff; ok
 20000: log =     4.3019052 (    4.3010300 actual); 0.0009 diff; ok
 20000: cbrt =    27.1449150 (   27.1441762 actual); 0.0007 diff; ok
 40000: sqrt =   199.9991455 (  200.0000000 actual); 0.0009 diff; ok
 40000: log =     4.6028333 (    4.6020600 actual); 0.0008 diff; ok
 40000: cbrt =    34.2003296 (   34.1995189 actual); 0.0008 diff; ok
 1e+08: sqrt =  9999.9994545 (10000.0000000 actual); 0.0005 diff; ok
 1e+08: log =     8.0009761 (    8.0000000 actual); 0.0010 diff; ok
 1e+08: cbrt =   464.1597912 (  464.1588834 actual); 0.0009 diff; ok
None

【问题讨论】：

是的，10**1000.0 太大了。
@KennyTM 其他学生的解决方案似乎通过了测试。我觉得我的解决方案缺少一些东西。
也许你不应该从 y=1000 开始。在所有情况下尝试 y=1。
这里是为什么你的教授代码运行得更好的提示。看看他的“find bounds”方法的目的，看看你为什么会溢出
对于像这样的逐次逼近算法，一个足够好的初始猜测通常是成功的一半。二进制搜索就像伊索乌龟：如果 lo 和 hi 将解决方案括起来，则缓慢而可靠。牛顿的方法更像是一只野兔。它可能会很快跳到解（良好的初始猜测），但它可能会跳向无穷大（错误的初始猜测）。

标签： python binary-search newtons-method

【解决方案1】：

这其实是你对数学而不是程序的理解有问题。该算法很好，但提供的初始条件不是。

你这样定义inverse(f, delta)：

def inverse(f, delta=1e-5):
    ...
    def newton(y, iters=15):
        guess = float(y)/2
        ...
    return newton

所以你猜 1000 = 10^x 的结果是 500.0，但肯定 10⁵⁰⁰ 太大了。初始猜测应该选择对 f 有效，而不是选择对 f 的逆。

我建议您使用 1 的猜测值进行初始化，即将该行替换为

guess = 1

它应该可以正常工作。

顺便说一句，您的二进制搜索的初始条件也是错误的，因为您假设解决方案在 0 和 y 之间：

low, high = 0, float(y)

这对您的测试用例来说是正确的，但构建反例很容易，例如log₁₀ 0.1 (= -1), √0.36 (= 0.6) 等（你教授的find_bounds 方法确实解决了√0.36 的问题，但仍然无法处理日志₁₀ 0.1 个问题。）

【讨论】：

他教授的 find bounds 方法试图解决您的第二点。
我想我自动假设，这是任何 mon 递增函数的一般反函数，将范围从中间拆分是一个很好的起点……应该多考虑一下。所以一般来说盯着 1 会是一个更好的起点，还是只针对 log10？
另外，任何关于为什么我的牛顿法对所有 log10 测试都差强人意的见解（上面的例子是我的 bin-search 方法）？还是这与牛顿的初始猜测问题有关？
@Verbal_Kint 我想最好在math.stackexchange.com 上询问有关这些数字求根方法的详细信息。

【解决方案2】：

我跟踪了您的错误，但基本上归结为 10**10000000 导致 python 溢出。使用数学库快速检查

math.pow(10,10000000)

Traceback (most recent call last):
  File "<pyshell#3>", line 1, in <module>
    math.pow(10,10000000)
OverflowError: math range error

我为你做了一些研究，发现了这个

Handling big numbers in code

您需要重新评估为什么需要计算如此大的数字（并相应地更改您的代码 :: 建议）或开始寻找一些更大的数字处理解决方案。

您可以编辑您的反函数以检查某些输入是否会导致它失败（try 语句），如果函数不是单调递增，这也可以解决零除法的一些问题，并避免这些区域或

您可以在关于 y=x 的“有趣”区域中镜像多个点，并通过这些点使用插值方案来创建“逆”函数（厄米特、泰勒级数等）。

【讨论】：

我不知道大小限制，但似乎不会影响班级论坛中发布的其他解决方案，这些解决方案通过了所有测试。此外，任何关于为什么我的牛顿方法根本无法处理 log10 的见解（答案很糟糕）？谢谢
math 模块本质上提供了对 C 的快速 <math.h> 的访问，因此对于内置 pow() 或 ** 运算符不会的值将引发 OverflowError 异常。也就是说，您的重点是正确的，因为即使没有引发异常，计算如此庞大的数字也会非常缓慢。
这里是为什么你的教授代码运行得更好的提示。看看他的“find bounds”方法的目的，看看你为什么会溢出。你可能连一次计算都没有完成

【解决方案3】：

如果您使用的是 Python 2.x，int 和 long 是不同的类型，OverflowError 只能导致 ints（q.v.，Built-in Exceptions）。尝试显式使用 longs 代替（通过在整数值上使用 long() 内置函数，或将 L 附加到数字文字）。

编辑：显然，正如 Paul Seeb 和 KennyTM 在他们出色的答案中指出的那样，这不是算法缺陷的补救措施。

【讨论】：