Python 中有效的求根演算法

作為數據科學家和計算機科學家，我們經常在日常生活中處理求根演算法，即使我們沒有意識到這一點。這些演算法被設計成定位特定值、局部/全局最大值或最小值的附近。

我們利用求根演算法來搜索特定值、局部/全局最大值或最小值的接近度。

在數學中，求根通常意味著我們試圖求解一個像 f(X) = 0 這樣的方程組。這將使求根演算法也成為非常有效的搜索演算法。我們要做的就是定義 g(X) = f(X) – Y，其中 Y 是搜索目標，而不是像 g(X) = f(X) – Y = 0 那樣求解 X。

這些方法主要分為兩個不同的家族:

括弧逼近(例如二等分演算法)
迭代方法(例如牛頓法、割線法、斯特芬森法等等)

在下面的教程中，我們將了解 Python 編程語言中一些演算法的實現，並將它們相互比較。這些演算法如下:

二分演算法
正則-法爾西演算法
伊利諾伊演算法
割線演算法
斯特芬森演算法

在我們開始之前，讓我們假設我們有一個連續的函數 f，並且想要搜索一個值 y。因此，我們正在求解方程:f(x) – y = 0。

理解等分演算法

二分演算法，也因其離散版本(二分搜索)或樹變體(二叉查找樹)而聞名，是一種在邊界內搜索目標值的有效演算法。因此，這種演算法也被稱為尋找演算法根的包圍方法。

關鍵強度:

二分演算法是一種穩健的演算法，保證了對目標值的合理收斂速度。

關鍵弱點:

該演算法需要知道根的估計面積。比如 3 ≤ π ≤ 4。
這個演算法很好用；估計區域中只有一個根。

假設我們知道 x 在 f(p)和 f(q)之間，這就形成了搜索括弧。演算法將檢查 x 是否大於或小於 f((p + q) / 2)，這是括弧的中點。

當搜索一個連續函數時，很可能，我們將永遠無法定位精確的值(例如，定位？).需要一定的誤差範圍來檢查支架的中點。當計算值接近目標值時，我們可以將誤差幅度視為提前停止。例如，如果誤差幅度為 0.001%，3.141624 是否足夠接近？，約 3.1415926…

如果計算值足夠接近目標值，則搜索完成，否則，如果 x 小於 f((p + q) / 2，則搜索下半部分的值，反之亦然。

現在讓我們考慮下面的代碼片段，它在 Python 中演示了同樣的內容。

示例:


def bisectionAlgorithm(f, p, q, y, margin = .00_001):
    ''' Bracketed approach of Root-finding with bisection method
    Arguments
    ----------
    f: callable, continuous function
    p: float, lower bound to be searched
    q: float, upper bound to be searched
    y: float, target value
    margin: float, error margin in absolute term

    Return Values
    ----------
    A float r, where f(r) is within the margin of y
    '''

    if (lower := f(p)) > (upper := f(q)):
        p, q = q, p
        lower, upper = upper, lower

    assert y >= lower, f"y is smaller than the lower bound. {y} < {lower}"
    assert y <= upper, f"y is larger than the upper bound. {y} > {upper}"

    while 1:
        r = (p + q) / 2
        if abs((y_r := f(r)) - y) < margin:
            # found!
            return r
        elif y < y_r:
            p, upper = r, y_r
        else:
            q, lower = r, y_r

說明:

在上面的代碼片段中，我們將一個函數定義為等分演算法，該演算法接受一些參數，如 f、p、q、y 和 margin ，其中 f 是一個可調用的連續函數， p 是一個浮點值和要搜索的下限， q 也是一個浮點值和要搜索的上限， y 又是一個浮點值和目標值然後我們使用 if 條件語句來檢查分配給 f(p) 的下界是否大於分配給 f(q) 的上界，對於這種情況， p 和 q 的值等於 q 和 p 並且上下等於上下。然後我們使用斷言關鍵字來調試 y 的值。然後我們使用了而循環，其中我們定義了 r 的值等於 p 和 q 的平均值。在這個循環中，我們還使用了 if-elif-else 條件語句來檢查分配給 f(r) – y 的 y_r 是否小於邊距，並返回相同的 r 。

理解虛假位置演算法

類似於二等分演算法，偽位置演算法，也稱為雷古拉法西，利用包圍方法。但是與二等分演算法不同，它並不使用每次迭代都將問題空間分成兩半的強力方法。相反，該演算法迭代地繪製一條從 f(p)到 f(q)的直線，並將截距與目標值進行比較。然而，不能保證搜索效率總是得到提高。例如，下圖描述了如何只有下限以降低的速率增加，而上限保持不變。

圖 1: 停滯邊界減緩收斂

關鍵強度:

這種演算法通常比二等分演算法收斂得更快。
Regula Falsi 的好處是，隨著括弧變小，連續函數將收斂到一條直線。

關鍵弱點:

當演算法達到停滯界限時，Regula Falsi 也顯示出較慢的收斂速度。
該演算法需要知道根的近似面積。比如 3 ≤ π ≤ 4。

Regula Falsi 和平分演算法在實現上的顯著區別是 r 不再是 p 和 q 的中點；然而，據估計:

讓我們考慮下面的代碼片段來證明這一點:

示例:


def regulaFalsiAlgorithm(f, p, q, y, margin = .00_001):
    ''' Bracketed approach of Root-finding with regula-falsi method
    Arguments
    ----------
    f: callable, continuous function
    p: float, lower bound to be searched
    q: float, upper bound to be searched
    y: float, target value
    margin: float, error margin in absolute term

    Return Values
    ----------
    A float r, where f(r) is within the margin of y
    '''

    assert y >= (lower := f(p)), f"y is smaller than the lower bound. {y} < {lower}"
    assert y <= (upper := f(q)), f"y is larger than the upper bound. {y} > {upper}"

    while 1:
        r = ((p * (upper - y)) - (q * (lower - y))) / (upper - lower)
        if abs((y_r := f(r)) - y) < margin:
            # found!
            return r
        elif y < y_r:
            q, upper = r, y_r
        else:
            p, lower = r, y_r

說明:

在上面的代碼片段中，我們定義了一個函數為 regulaFalsiAlgorithm ，該函數接受一些參數，如 f、p、q、y 和 margin ，其中 f 是一個可調用的連續函數， p 是一個浮點值和要搜索的下限， q 也是一個浮點值和要搜索的上限， y 又是一個浮點值然後我們使用斷言關鍵字來調試 y 的值。然後我們使用了和循環，在該循環中我們定義了 r 的值。在這個循環中，我們還使用了 if-elif-else 條件語句來檢查分配給 f(r) – y 的 y_r 是否小於邊距，並為其返回 r 。

理解伊利諾伊演算法(修正的雷古拉·法西)

為了越過停滯邊界，我們可以插入一個條件規則，當一個邊界保持停滯兩輪。就拿前面的例子來說，由於 q 已經兩輪沒有移動了，並且 r 還沒有靠近根 x，所以在下一輪的時候，線會畫到 f(q) / 2，而不是 f(q)。如果下界是停滯界，那麼下界也是如此。

圖 2: 伊利諾伊演算法避免了長時間停滯，收斂速度更快。

關鍵強度:

伊利諾伊演算法比二等分演算法和雷古拉法西演算法收斂更快。
我們可以通過將停滯界限與目標值的距離減半來避免停滯界限。

關鍵弱點:

當演算法達到停滯界限時，該演算法也表現出較慢的收斂。
該演算法需要知道根的估計面積。比如 3 ≤ π ≤ 4。

示例:


def illinoisAlgorithm(f, p, q, y, margin = .00_001):
    ''' Bracketed approach of Root-finding with illinois method
    Arguments
    ----------
    f: callable, continuous function
    p: float, lower bound to be searched
    q: float, upper bound to be searched
    y: float, target value
    margin: float, error margin in absolute term

    Return Values
    ----------
    A float r, where f(r) is within the margin of y
    '''

    assert y >= (lower := f(p)), f"y is smaller than the lower bound. {y} < {lower}"
    assert y <= (upper := f(q)), f"y is larger than the upper bound. {y} > {upper}"

    stagnant = 0

    while 1:
        r = ((p * (upper - y)) - (q * (lower - y))) / (upper - lower)
        if abs((y_r := f(r)) - y) < margin:
            # found!
            return r
        elif y < y_r:
            q, upper = r, y_r
            if stagnant == -1:
                # Lower bound is stagnant!
                lower += (y - lower) / 2
            stagnant = -1
        else:
            p, lower = r, y_r
            if stagnant == 1:
                # Upper bound is stagnant!
                upper -= (upper - y) / 2
            stagnant = 1

說明:

在上面的代碼片段中，我們定義了一個函數為 illinoisAlgorithm ，它接受一些參數，如 f、p、q、y 和 margin ，其中 f 是一個可調用的連續函數， p 是一個浮點值和要搜索的下限， q 也是一個浮點值和要搜索的上限， y 又是一個浮點值和目標值，【T12 然後我們使用斷言關鍵字來調試 y 的值。然後，我們將變數定義為等於零的停滯。然後我們使用了和循環，在該循環中我們定義了 r 的值。在這個循環中，我們還使用了 if-elif-else 條件語句來檢查分配給 f(r) – y 的 y_r 是否小於邊距，並為其返回 r 。

理解割線法(擬牛頓法)

到目前為止，我們一直在實施括弧方法。如果我們不知道括弧的位置呢？在這種情況下，割線法可能會有所幫助。割線法是一種迭代演算法，從兩個值開始，並試圖向目標值收斂。雖然我們可以在演算法收斂的同時獲得更好的性能，並且不需要知道粗略的根位置，但是如果兩個初始值離真正的根太遠，我們可能會遇到發散的風險。

關鍵強度:

割線方法不需要由根組成的括弧。
這種方法不需要知道根的估計面積。

關鍵弱點:

1.與所有早期的方法不同，割線法不能保證收斂性。

割線方法從檢查兩個用戶定義的種子開始。假設我們想求 x 2 – math.pi = 0 的根，從 x_0 = 4 和 x_1 = 5 開始；我們的種子將分別是 4 粒和 5 粒。(注:這個過程類似於像 x 一樣搜索 x 2 = math.pi)

圖 3: 割線基於 x1 和 x2 定位 x3 的方法

然後，我們通過在 f(x_0)和 f(x_1)之間畫一條直線來定位目標值為 x_2 的截距，就像我們在 Regula Falsi 演算法中所做的那樣。如果 f(x_2)沒有足夠接近目標值，我們必須重複該步驟並定位 x_3。通常，我們可以使用以下公式計算下一個 x:

讓我們考慮下面的代碼片段來證明這一點:

示例:


def secantAlgorithm(f, x0, x1, y, iterations, margin = .00_001):
    ''' Iterative approach of Root-finding with secant method
    Arguments
    ----------
    f: callable, continuous function
    x0: float, initial seed
    x1: float, initial seed
    y: float, target value
    iterations: int, maximum number of iterations to avoid indefinite divergence
    margin: float, margin of error in absolute term
    Return Values
    ----------
    A float x2, where f(x2) is within the margin of y
    '''

    assert x0 != x1, "Two different initial seeds are required."

    if abs((y0 := f(x0) - y)) < margin:
        # found!
        return x0
    if abs((y1 := f(x1) - y)) < margin:
        # found!
        return x1

    for i in range(iterations):
        x2 = x1 - y1 * (x1 - x0) / (y1 - y0)
        if abs((y2 := f(x2) - y)) < margin:
            # found!
            return x2
        x0, x1 = x1, x2
        y0, y1 = y1, y2
    return x2

說明:

在上面的代碼片段中，我們將一個函數定義為 secantAlgorithm ，該函數接受一些參數，如 f、x0、x1、y、迭代和 margin ，其中 f 是一個可調用的連續函數， x0 和 x1 是浮點值和初始種子， y 又是一個浮點值，目標值，迭代為整數值，存儲最大迭代次數的值，避免不定發散，餘量為絕對項誤差餘量，也是浮點值。然後我們使用斷言關鍵字來檢查 x0 的值是否不等於 x1 的值。然後，我們使用 if 條件語句檢查分配給 f(x0) – y 的 y0 是否小於邊距變數，並返回相同的 x0 變數。我們再次使用 if 條件語句檢查分配給 f(x1) – y 的 y1 是否小於邊距變數，並返回相同的 x1 變數。最後，我們使用【T44 for】循環來迭代存儲在迭代變數中的值，並定義公式來尋找根。在循環中，我們再次使用了 if 條件語句並返回 x2 。

理解斯特芬森方法

割線法通過消除由根組成的括弧的要求，進一步改進了 Regula Falsi 演算法。讓我們回想一下，直線只是兩個 x 值的切線(即導數)的一個天真值(或者在 Regula Falsi 和 Illinois 演算法的情況下是上下界)。隨著搜索的收斂，這個值將是完美的。在斯特芬森演算法中，我們將嘗試用更好的導數值來代替直線，以進一步改進割線法。

關鍵強度:

Steffensen 的方法不需要包含根的括弧。
該方法也不需要根的估計面積的知識。
如果可能的話，這種方法比割線法收斂得更快。

關鍵弱點:

如果初始種子離真正的根太遠，Steffensen 的方法不能保證收斂。
連續函數的求值將是割線法的兩倍，以便更好地計算導數。

我們可以在斯特芬森演算法的幫助下，通過基於用戶定義的初始種子 x0 計算以下內容來更好地估計導數。

這相當於下面的式，其中 h = f(x):

取 h 到 0 的極限，我們將得到？(?).

然後，我們將使用廣義評估斜率函數，按照與割線法相同的步驟定位下一步:

讓我們考慮下面的代碼片段來證明這一點:

示例:


def steffensenAlgorithm(f, x, y, iterations, margin = .00_001):
    ''' Iterative approach of Root-finding with steffensen's method
    Arguments
    ----------
    f: callable, continuous function
    x: float, initial seed
    y: float, target value
    iterations: int, maximum number of iterations to avoid indefinite divergence
    margin: float, error margin in absolute term
    Return Values
    ----------
    A float x2, where f(x2) is within the margin of y
    '''

    assert x != 0, "Initial seed can't be zero."

    if abs((yx := f(x) - y)) < margin:
        # found!
        return x

    for n in range(iterations):
        g = (f(x + yx) - y) / yx - 1
        if g * x == 0:
            # Division by zero, search stops
            return x
        x -= (f(x) - y) / (g * x)
        if abs((yx := f(x) - y)) < margin:
            # found!
            return x
    return x

說明:

在上面的代碼片段中，我們定義了一個函數為 secantAlgorithm ，它接受一些參數，如 f，x0，x1，y，迭代，和 margin ，其中 f 是一個可調用的連續函數， x0 和 x1 是浮點值和初始種子， y 又是一個浮點值，目標值，迭代為整數值，存儲最大迭代次數的值，避免不定發散，餘量為絕對項誤差餘量，也是浮點值。然後我們使用斷言關鍵字來檢查初始種子是否不等於零。然後我們使用 if 條件語句檢查分配給 f(x) – y 的 yx 是否小於邊距變數，並返回相同的 x 變數。最後，我們使用【T30 for-循環來迭代存儲在迭代變數中的值，並定義公式來尋找根。在循環中，我們再次使用了 if 條件語句並返回 x 。

結論

在上面的教程中，我們已經了解了用於搜索根的以下五種演算法的優點、缺點和實現。

二分演算法
雷古拉法爾西演算法
伊利諾伊演算法
割線演算法
斯特芬森演算法

現在讓我們考慮下表，它顯示了我們已經實現的演算法的比較。

| | 二分演算法 | 雷古拉法爾西演算法 | 伊利諾伊演算法 | 割線演算法 | 斯特芬森演算法 |
| 接近 | 支架 | 支架 | 支架 | 重複的 | 重複的 |
| 收斂 | 擔保 | 擔保 | 擔保 | 不保證 | 不保證 |
| 根部大致位置的知識 | 需要 | 需要 | 需要 | 不需要 | 不需要 |
| 初始種子數量 | Two | Two | Two | Two | one |
| 緩慢收斂的風險 | 無法使用 | 停滯界限 | 無法使用 | 初始種子離根不夠近 | 初始種子離根不夠近 |
| 減少問題空間的方法 | 暴力減半 | 有限差分近似導數 | 有限差分近似導數 | 有限差分近似導數 | 有限差分近似導數 |
| 收斂速度 | 線性的 | 線性的 | 超線性的 | 超線性的 | 二次的 |

一旦我們熟悉了這些演算法，本教程就沒有涉及到許多其他要探索的尋根演算法。其中一些是牛頓-拉夫森法、逆二次插值、布倫特法等等。繼續探索，將上述演算法添加到工具庫中。

原創文章，作者：小藍，如若轉載，請註明出處：https://www.506064.com/zh-tw/n/257991.html

Python 中有效的求根演算法

理解等分演算法

理解虛假位置演算法

理解伊利諾伊演算法(修正的雷古拉·法西)

理解割線法(擬牛頓法)

理解斯特芬森方法

結論

相關推薦

發表回復