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

Question

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

以下算法不收敛：

而以下确实会收敛：

您可以假设函数 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)

注意：我知道浮点数是如何工作的（请不要 post 你最喜欢的 link 网站告诉我永远不要比较两个浮点数）。另外，我不是在寻找问题的解决方案，而是在寻找一种解释为什么其中一种算法收敛而另一种算法收敛。

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

Answer 1

"I'm aware of how floating point numbers work..."。也许浮点运算的工作原理比想象的要复杂。

这是使用牛顿法循环迭代的经典示例。将差值与 epsilon 进行比较是 "mathematical thinking" 并且在使用浮点数时可能会烧伤您。在您的示例中，您访问了 x 的几个浮点值，然后您陷入了两个数字之间的循环。 "floating-point thinking" 更好地表述如下（抱歉，我的首选语言是 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);
}

Answer 2

我认为尝试强制完全相等（而不是 err < small）总是会经常失败。在您的示例中，对于 1 到 10 之间的 100,000 个随机数（而不是您的 2.0），第一种方法大约有 1/3 的时间失败，第二种方法大约有 1/6 的时间失败。我敢打赌有一种方法可以预测！

这需要大约 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()

Answer 3

None的方法很完美：

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

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

Answer 4

I'm trying to figure out why one of the implementations of Newton's method is converging exactly within floating point precision, wheres another is not.

从技术上讲，它不会收敛到正确的值。尝试打印更多数字，或使用 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.

时产生的舍入误差决定

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

Different implementations of Newton's method in floating point arithmetic

python

math

floating-point

newtons-method