牛顿法在浮点运算中的不同实现答案

【问题标题】：Different implementations of Newton's method in floating point arithmetic牛顿法在浮点运算中的不同实现
【发布时间】：2015-06-11 11:38:20
【问题描述】：

我正在用牛顿法求解一维非线性方程。我试图弄清楚为什么牛顿方法的一种实现精确地收敛在浮点精度内，而另一种则没有。

以下算法不收敛：

而以下确实收敛：

您可以假设函数 f 和 f' 是平滑且表现良好的。我能想到的最好的解释是，这在某种程度上与所谓的迭代改进有关（Golub 和 Van Loan，1989 年）。任何进一步的见解将不胜感激！

这是一个简单的 Python 示例来说明问题

# Python
def f(x):
    return x*x-2.

def fp(x):
    return 2.*x

xprev = 0.

# converges
x = 1. # guess
while x != xprev:
    xprev = x
    x = (x*fp(x)-f(x))/fp(x)
print(x)

# does not converge
x = 1. # guess
while x != xprev:
    xprev = x
    dx = -f(x)/fp(x)
    x = x + dx
print(x)

注意：我知道浮点数的工作原理（请不要将您最喜欢的链接发布到告诉我永远不要比较两个浮点数的网站）。另外，我不是在寻找问题的解决方案，而是想解释为什么其中一种算法收敛而另一种不收敛。

更新： 正如@uhoh 指出的那样，在很多情况下第二种方法不会收敛。但是，我仍然不知道为什么第二种方法在我的现实世界场景中比第一种更容易收敛。所有的测试用例都有非常简单的函数f，而现实世界的f 有几百行代码（这就是我不想发布的原因）。所以也许f 的复杂性很重要。如果您对此有任何其他见解，请告诉我！

【问题讨论】：

你什么时候在第一个算法中分配Xprev？如果是在评估dx 之前，那么在代数上两者是相同的。也许非常小的f'(x) 会溢出？
错字。现已更正。
f'(x) 通常为 ~1，x 为 ~0.1。所以我认为不会在任何地方发生溢出。
从代数上讲，您的两个公式是相同的，因此要么实现中存在错误，要么您永远不会完全获得相同的浮点数。通常，您通过将两个浮点数的差异与一个非常小数进行比较来比较它们。
是的，通常我会将它与一个非常小的数字进行比较，这很好。但是，我仍然对为什么第二种算法完全收敛感到困惑（对于许多具有不同初始 x 和 f 的不同运行，即它不是一组特定数字取消的巧合）。

标签： python math floating-point newtons-method

【解决方案1】：

“我知道浮点数是如何工作的……”。也许浮点运算的工作原理比想象的要复杂。

这是使用牛顿法循环迭代的经典示例。将差异与 epsilon 进行比较是“数学思维”，在使用浮点时可能会烧毁你。在您的示例中，您访问了 x 的几个浮点值，然后您被困在两个数字之间的循环中。 “浮点思维”最好表述如下（抱歉，我的首选语言是 C++）

std::set<double> visited;
xprev = 0.0;
x = 1.0;
while (x != prev)
{
    xprev = x;
    dx = -F(x)/DF(x);
    x = x + dx;
    if (visited.find(x) != visited.end())
    {
        break;  // found a cycle
    }
    visited.insert(x);
}

【讨论】：

我明白这一点，但为什么其他实现永远不会发生这种情况？
(x*F'(x)-F(x))/F'(x) 和 x-F(x)/F'(x) 是使用实数算术的等价表达式。它们不是使用浮点运算的等效表达式。使用调试器单步执行并分析这两个表达式，看看它们可能不同。您的第一个实现的浮点舍入行为对您来说似乎“很好”。第二个实现的循环显然是一个浮点数。函数 F(x) 是凸的并且你接近一个根，所以理论上牛顿应该收敛。但浮点算术不是真正的算术......

【解决方案2】：

我认为试图强制一个完全相等（而不是 err

运行大约需要 30 秒，结果很可爱！：

def f(x, a):
    return x*x - a

def fp(x):
    return 2.*x

def A(a):
    xprev = 0.
    x = 1.
    n = 0
    while x != xprev:
        xprev = x
        x = (x * fp(x) - f(x,a)) / fp(x)
        n += 1
        if n >100:
            return n, x
    return n, x



def B(a):
    xprev = 0.
    x = 1.
    n = 0
    while x != xprev:
        xprev = x
        dx = - f(x,a) / fp(x)
        x = x + dx
        n += 1
        if n >100:
            return n, x
    return n, x

import numpy as np
import matplotlib.pyplot as plt


n = 100000

aa = 1. + 9. * np.random.random(n)

data_A = np.zeros((2, n))
data_B = np.zeros((2, n))

for i, a in enumerate(aa):
    data_A[:,i] = A(a)
    data_B[:,i] = B(a)

bins = np.linspace(0, 110, 12)

hist_A = np.histogram(data_A, bins=bins)
hist_B = np.histogram(data_B, bins=bins)
print "A: n<10: ", hist_A[0][0], " n>=100: ", hist_A[0][-1]
print "B: n<10: ", hist_B[0][0], " n>=100: ", hist_B[0][-1]

plt.figure()
plt.subplot(1,2,1)
plt.scatter(aa, data_A[0])
plt.subplot(1,2,2)
plt.scatter(aa, data_B[0])
plt.show()

【讨论】：

这就是我想了解的！在我感兴趣的实际案例中（它的发布方式很长），第一种方法似乎在数百万个测试用例中每次都收敛。也许这取决于函数f的复杂程度？
可能是。我刚刚在程序中添加了一个情节，有趣的是，第二种方法往往不会在 4 到 8 之间失败 - 看图！我希望我可以发布图形（没有足够高的分数）其他人可以吗？

【解决方案3】：

没有一种方法是完美的：

两种方法都会失败的一种情况是，如果根正好位于两个连续浮点数 f1 和 f2 的中间。然后，两种方法都到达 f1，将尝试计算该中间值，并很有可能找到 f2，反之亦然。

/f(x) / / / / f1 / --+----------+------> x /f2 / / /

【讨论】：

是不是几乎总是这样（包括上面的例子）？
@hanno 不，根也可以接近浮点值，这使得x 与xprev 计算相同变得更容易。但是话又说回来，如果 f 的计算中有重要的近似值，正如 Simon 指出的那样，需要了解计算根的函数是什么。

【解决方案4】：

我试图弄清楚为什么牛顿方法的一种实现精确地收敛在浮点精度范围内，而另一种则没有。

从技术上讲，它不会收敛到正确的值。尝试打印更多数字，或使用float.hex。

第一个给出

>>> print "%.16f" % x
1.4142135623730949
>>> float.hex(x)
'0x1.6a09e667f3bccp+0'

而正确舍入的值是下一个浮点值：

>>> print "%.16f" % math.sqrt(2)
1.4142135623730951
>>> float.hex(math.sqrt(2))
'0x1.6a09e667f3bcdp+0'

第二种算法实际上是在两个值之间交替，所以不会收敛。

问题是由于f(x)中的灾难性取消：因为x*x将非常接近2，当你减去2时，结果将由计算中产生的舍入误差决定x*x.

【讨论】：