Ctrl+K
scipy.optimize.

least_squares#

scipy.optimize.least_squares(fun, x0, jac='2-point', bounds=(-inf, inf), method='trf', ftol=1e-08, xtol=1e-08, gtol=1e-08, x_scale=1.0, loss='linear', f_scale=1.0, diff_step=None, tr_solver=None, tr_options=None, jac_sparsity=None, max_nfev=None, verbose=0, args=(), kwargs=None)[原始碼]#

解決具有變數邊界的非線性最小平方法問題。

給定殘差 f(x) (n 個實變數的 m 維實函數) 和損失函數 rho(s) (純量函數)，least_squares 尋找成本函數 F(x) 的局部最小值

minimize F(x) = 0.5 * sum(rho(f_i(x)**2), i = 0, ..., m - 1)
subject to lb <= x <= ub

損失函數 rho(s) 的目的是減少離群值對解的影響。

參數:
fun可呼叫物件

計算殘差向量的函數，具有簽名 fun(x, *args, **kwargs)，即，相對於其第一個參數進行最小化。傳遞給此函數的參數 x 是一個形狀為 (n,) 的 ndarray (永遠不是純量，即使對於 n=1 也是如此)。它必須分配並傳回形狀為 (m,) 的一維類陣列或純量。如果參數 x 是複數或函數 fun 傳回複數殘差，則必須將其包裝在實參數的實函數中，如範例章節末尾所示。

x0形狀為 (n,) 或 float 的類陣列

獨立變數的初始猜測值。如果為 float，則將其視為具有一個元素的一維陣列。當 method 為 ‘trf’ 時，初始猜測值可能會稍微調整，使其充分位於給定的 bounds 內。

jac{‘2-point’, ‘3-point’, ‘cs’, 可呼叫物件}, 選用

計算 Jacobian 矩陣的方法 (一個 m 乘 n 的矩陣，其中元素 (i, j) 是 f[i] 相對於 x[j] 的偏導數)。關鍵字選擇用於數值估計的有限差分方案。方案 ‘3-point’ 更準確，但需要的運算量是 ‘2-point’ (預設) 的兩倍。方案 ‘cs’ 使用複數步長，雖然可能最準確，但僅適用於當 fun 正確處理複數輸入並且可以解析地延拓到複數平面時。方法 ‘lm’ 始終使用 ‘2-point’ 方案。如果為可呼叫物件，則將其用作 jac(x, *args, **kwargs)，並且應傳回 Jacobian 的良好近似值 (或精確值)，形式為類陣列 (套用 np.atleast_2d)、稀疏矩陣 (csr_matrix 為效能首選) 或 scipy.sparse.linalg.LinearOperator。

bounds2 元組的類陣列或 Bounds, 選用

有兩種方法可以指定邊界

  1. Bounds 類別的實例

  2. 獨立變數的下限和上限。預設為無邊界。每個陣列都必須符合 x0 的大小，或者是純量，在後一種情況下，邊界對於所有變數都相同。使用具有適當符號的 np.inf 來停用所有或某些變數的邊界。

method{‘trf’, ‘dogbox’, ‘lm’}, 選用

執行最小化的演算法。

  • ‘trf’ : 信賴域反射演算法，特別適用於具有邊界的大型稀疏問題。通常是穩健的方法。

  • ‘dogbox’ : 具有矩形信賴域的 dogleg 演算法，典型用例是具有邊界的小型問題。不建議用於 Jacobian 秩不足的問題。

  • ‘lm’ : MINPACK 中實作的 Levenberg-Marquardt 演算法。不處理邊界和稀疏 Jacobian。通常是小型無約束問題的最有效方法。

預設值為 ‘trf’。請參閱 Notes 以取得更多資訊。

ftolfloat 或 None, 選用

成本函數變更的終止容忍度。預設值為 1e-8。當 dF < ftol * F 時，最佳化程序會停止，並且在最後一步中，局部二次模型與真實模型之間存在足夠的一致性。

如果為 None 且 ‘method’ 不是 ‘lm’，則會停用依此條件終止。如果 ‘method’ 為 ‘lm’，則此容忍度必須高於機器 epsilon。

xtolfloat 或 None, 選用

獨立變數變更的終止容忍度。預設值為 1e-8。確切條件取決於使用的 method

  • 對於 ‘trf’ 和 ‘dogbox’ : norm(dx) < xtol * (xtol + norm(x))。

  • 對於 ‘lm’ : Delta < xtol * norm(xs)，其中 Delta 是信賴域半徑，而 xs 是根據 x_scale 參數 (請參閱下文) 縮放的 x 值。

如果為 None 且 ‘method’ 不是 ‘lm’，則會停用依此條件終止。如果 ‘method’ 為 ‘lm’，則此容忍度必須高於機器 epsilon。

gtolfloat 或 None, 選用

梯度範數的終止容忍度。預設值為 1e-8。確切條件取決於使用的 method

  • 對於 ‘trf’ : norm(g_scaled, ord=np.inf) < gtol，其中 g_scaled 是縮放的梯度值，以考量邊界的存在 [STIR]。

  • 對於 ‘dogbox’ : norm(g_free, ord=np.inf) < gtol，其中 g_free 是相對於不在邊界上的最佳狀態的變數的梯度。

  • 對於 ‘lm’ : Jacobian 的列與殘差向量之間角度的餘弦最大絕對值小於 gtol，或殘差向量為零。

如果為 None 且 ‘method’ 不是 ‘lm’，則會停用依此條件終止。如果 ‘method’ 為 ‘lm’，則此容忍度必須高於機器 epsilon。

x_scale類陣列或 ‘jac’, 選用

每個變數的特徵尺度。設定 x_scale 等效於以縮放變數 xs = x / x_scale 重新制定問題。另一種觀點是，沿著第 j 個維度的信賴域大小與 x_scale[j] 成正比。透過設定 x_scale，使得沿著任何縮放變數的給定大小的步長對成本函數具有相似的影響，可以實現更好的收斂性。如果設定為 ‘jac’，則使用 Jacobian 矩陣列的反範數迭代更新尺度 (如 [JJMore] 中所述)。

lossstr 或可呼叫物件, 選用

決定損失函數。允許以下關鍵字值

  • ‘linear’ (預設) : rho(z) = z。給出標準最小平方法問題。

  • ‘soft_l1’ : rho(z) = 2 * ((1 + z)**0.5 - 1)。l1 (絕對值) 損失的平滑近似值。通常是穩健最小平方法的良好選擇。

  • ‘huber’ : rho(z) = z if z <= 1 else 2*z**0.5 - 1。作用類似於 ‘soft_l1’。

  • ‘cauchy’ : rho(z) = ln(1 + z)。嚴重削弱離群值的影響，但可能會在最佳化過程中造成困難。

  • ‘arctan’ : rho(z) = arctan(z)。限制單個殘差的最大損失，具有與 ‘cauchy’ 類似的屬性。

如果為可呼叫物件，則它必須採用一維 ndarray z=f**2 並傳回形狀為 (3, m) 的類陣列，其中第 0 列包含函數值，第 1 列包含一階導數，第 2 列包含二階導數。方法 ‘lm’ 僅支援 ‘linear’ 損失。

f_scalefloat, 選用

離群值殘差和離群值殘差之間的軟邊界值，預設值為 1.0。損失函數的評估方式如下 rho_(f**2) = C**2 * rho(f**2 / C**2)，其中 C 是 f_scale，而 rho 由 loss 參數決定。此參數對 loss='linear' 沒有影響，但對於其他 loss 值，它至關重要。

max_nfevNone 或 int, 選用

終止前的最大函數評估次數。如果為 None (預設值)，則會自動選擇值

  • 對於 ‘trf’ 和 ‘dogbox’ : 100 * n。

  • 對於 ‘lm’ : 如果 jac 是可呼叫物件，則為 100 * n，否則為 100 * n * (n + 1) (因為 ‘lm’ 會計算 Jacobian 估計中的函數呼叫次數)。

diff_stepNone 或類陣列, 選用

決定 Jacobian 有限差分近似的相對步長大小。實際步長計算為 x * diff_step。如果為 None (預設值)，則 diff_step 會被視為用於有限差分方案的傳統「最佳」機器 epsilon 冪 [NR]。

tr_solver{None, ‘exact’, ‘lsmr’}, 選用

用於解決信賴域子問題的方法，僅與 ‘trf’ 和 ‘dogbox’ 方法相關。

  • ‘exact’ 適用於 Jacobian 矩陣稠密且不太大的問題。每次迭代的計算複雜度與 Jacobian 矩陣的奇異值分解相當。

  • ‘lsmr’ 適用於 Jacobian 矩陣稀疏且大的問題。它使用迭代程序 scipy.sparse.linalg.lsmr 來尋找線性最小平方法問題的解，並且僅需要矩陣向量乘積評估。

如果為 None (預設值)，則會根據第一次迭代時傳回的 Jacobian 類型來選擇求解器。

tr_optionsdict, 選用

傳遞給信賴域求解器的關鍵字選項。

  • tr_solver='exact': tr_options 會被忽略。

  • tr_solver='lsmr': scipy.sparse.linalg.lsmr 的選項。此外，method='trf' 支援 ‘regularize’ 選項 (bool, 預設值為 True)，它會將正規化項新增至正規方程式，這會在 Jacobian 秩不足時改善收斂性 [Byrd] (方程式 3.4)。

jac_sparsity{None, 類陣列, 稀疏矩陣}, 選用

定義 Jacobian 矩陣的稀疏結構，用於有限差分估計，其形狀必須為 (m, n)。如果 Jacobian 在每列中只有少數非零元素，則提供稀疏結構將大大加快計算速度 [Curtis]。零項目表示 Jacobian 中的對應元素恆為零。如果提供，則強制使用 ‘lsmr’ 信賴域求解器。如果為 None (預設值)，則將使用稠密差分。對 ‘lm’ 方法沒有影響。

verbose{0, 1, 2}, 選用

演算法的詳細程度

  • 0 (預設值) : 靜默工作。

  • 1 : 顯示終止報告。

  • 2 : 在迭代期間顯示進度 (‘lm’ 方法不支援)。

args, kwargstuple 和 dict, 選用

傳遞給 fun 和 jac 的其他引數。預設情況下皆為空。呼叫簽名為 fun(x, *args, **kwargs)，jac 也是如此。

傳回值:
resultOptimizeResult

具有以下已定義欄位的 OptimizeResult

xndarray, 形狀 (n,)

找到的解。

costfloat

解處的成本函數值。

funndarray, 形狀 (m,)

解處的殘差向量。

jacndarray, 稀疏矩陣或 LinearOperator, 形狀 (m, n)

解處的修改 Jacobian 矩陣，從某種意義上來說，J^T J 是成本函數 Hessian 的 Gauss-Newton 近似值。類型與演算法使用的類型相同。

gradndarray, 形狀 (m,)

解處的成本函數梯度。

optimalityfloat

一階最佳化測量。在無約束問題中，它始終是梯度的均勻範數。在約束問題中，它是與迭代期間與 gtol 比較的量。

active_maskndarray of int, 形狀 (n,)

每個元件都顯示對應的約束是否處於活動狀態 (也就是說，變數是否在邊界上)

  • 0 : 約束未處於活動狀態。

  • -1 : 下限處於活動狀態。

  • 1 : 上限處於活動狀態。

對於 ‘trf’ 方法，可能會有些任意，因為它會產生嚴格可行的迭代序列，並且 active_mask 是在容忍度閾值內決定的。

nfevint

完成的函數評估次數。方法 ‘trf’ 和 ‘dogbox’ 不會計算數值 Jacobian 近似的函數呼叫次數，這與 ‘lm’ 方法相反。

njevint 或 None

完成的 Jacobian 評估次數。如果在 ‘lm’ 方法中使用數值 Jacobian 近似，則設定為 None。

statusint

演算法終止的原因

  • -1 : 從 MINPACK 傳回的不正確輸入參數狀態。

  • 0 : 超過最大函數評估次數。

  • 1 : 符合 gtol 終止條件。

  • 2 : 符合 ftol 終止條件。

  • 3 : 符合 xtol 終止條件。

  • 4 : 符合 ftol 和 xtol 終止條件。

messagestr

終止原因的文字描述。

successbool

如果符合其中一個收斂準則 ( status > 0)，則為 True。

另請參閱

leastsq

Levenberg-Marquadt 演算法 MINPACK 實作的舊版包裝器。

curve_fit

應用於曲線擬合問題的最小平方法最小化。

註解

方法 ‘lm’ (Levenberg-Marquardt) 呼叫 MINPACK (lmder, lmdif) 中實作的最小平方法演算法的包裝器。它執行 Levenberg-Marquardt 演算法，該演算法制定為信賴域類型演算法。此實作基於論文 [JJMore]，它非常穩健且有效，具有許多巧妙的技巧。它應該是您無約束問題的首選。請注意，它不支援邊界。此外，當 m < n 時，它無法運作。

方法 ‘trf’ (信賴域反射) 的動機是解決方程式系統的過程，該方程式系統構成了 [STIR] 中制定的邊界約束最小化問題的一階最佳化條件。該演算法迭代地解決信賴域子問題，該子問題由特殊的對角二次項增強，並且信賴域形狀由與邊界的距離和梯度的方向決定。這些增強功能有助於避免直接步入邊界並有效探索整個變數空間。為了進一步改善收斂性，該演算法會考慮從邊界反射的搜尋方向。為了遵守理論要求，該演算法使迭代嚴格可行。對於稠密 Jacobian，信賴域子問題由一種非常類似於 [JJMore] (並在 MINPACK 中實作) 中描述的精確方法來解決。與 MINPACK 實作的不同之處在於，Jacobian 矩陣的奇異值分解在每次迭代中完成一次，而不是 QR 分解和 Givens 旋轉消除序列。對於大型稀疏 Jacobian，使用 2 維子空間方法來解決信賴域子問題 [STIR], [Byrd]。子空間由縮放梯度和由 scipy.sparse.linalg.lsmr 提供的近似 Gauss-Newton 解所跨越。當未施加約束時，該演算法與 MINPACK 非常相似，並且通常具有可比較的效能。該演算法在無邊界和有邊界問題中都非常穩健地運作，因此被選為預設演算法。

方法 ‘dogbox’ 在信賴域框架中運作，但考慮矩形信賴域，而不是傳統的橢球體 [Voglis]。目前信賴域與初始邊界的交集再次為矩形，因此在每次迭代中，都會透過 Powell 的 dogleg 方法 [NumOpt] 近似地解決受邊界約束的二次最小化問題。對於稠密 Jacobian，可以精確計算所需的 Gauss-Newton 步長，或者對於大型稀疏 Jacobian，可以透過 scipy.sparse.linalg.lsmr 近似計算。當 Jacobian 的秩小於變數數時，該演算法可能會表現出收斂速度緩慢。在具有少量變數的有邊界問題中，該演算法通常優於 ‘trf’。

穩健損失函數的實作如 [BA] 中所述。其想法是在每次迭代時修改殘差向量和 Jacobian 矩陣，使得計算出的梯度和 Gauss-Newton Hessian 近似值與成本函數的真實梯度和 Hessian 近似值相符。然後，演算法以正常方式繼續進行，也就是說，穩健損失函數實作為標準最小平方法演算法的簡單包裝器。

在版本 0.17.0 中新增。

參考文獻

[STIR] (1,2,3)

M. A. Branch, T. F. Coleman, and Y. Li, “A Subspace, Interior, and Conjugate Gradient Method for Large-Scale Bound-Constrained Minimization Problems,” SIAM Journal on Scientific Computing, Vol. 21, Number 1, pp 1-23, 1999.

[NR]

William H. Press et. al., “Numerical Recipes. The Art of Scientific Computing. 3rd edition”, Sec. 5.7.

[Byrd] (1,2)

R. H. Byrd, R. B. Schnabel and G. A. Shultz, “Approximate solution of the trust region problem by minimization over two-dimensional subspaces”, Math. Programming, 40, pp. 247-263, 1988.

[Curtis]

A. Curtis, M. J. D. Powell, and J. Reid, “On the estimation of sparse Jacobian matrices”, Journal of the Institute of Mathematics and its Applications, 13, pp. 117-120, 1974.

[JJMore] (1,2,3)

J. J. More, “The Levenberg-Marquardt Algorithm: Implementation and Theory,” Numerical Analysis, ed. G. A. Watson, Lecture Notes in Mathematics 630, Springer Verlag, pp. 105-116, 1977.

[Voglis]

C. Voglis and I. E. Lagaris, “A Rectangular Trust Region Dogleg Approach for Unconstrained and Bound Constrained Nonlinear Optimization”, WSEAS International Conference on Applied Mathematics, Corfu, Greece, 2004.

[NumOpt]

J. Nocedal and S. J. Wright, “Numerical optimization, 2nd edition”, Chapter 4.

[BA]

B. Triggs et. al., “Bundle Adjustment - A Modern Synthesis”, Proceedings of the International Workshop on Vision Algorithms: Theory and Practice, pp. 298-372, 1999.

範例

在此範例中，我們找到 Rosenbrock 函數在沒有獨立變數邊界時的最小值。

>>> import numpy as np
>>> def fun_rosenbrock(x):
...     return np.array([10 * (x[1] - x[0]**2), (1 - x[0])])

請注意，我們僅提供殘差向量。演算法會將成本函數建構為殘差平方和，這會產生 Rosenbrock 函數。精確最小值位於 x = [1.0, 1.0]。

>>> from scipy.optimize import least_squares
>>> x0_rosenbrock = np.array([2, 2])
>>> res_1 = least_squares(fun_rosenbrock, x0_rosenbrock)
>>> res_1.x
array([ 1.,  1.])
>>> res_1.cost
9.8669242910846867e-30
>>> res_1.optimality
8.8928864934219529e-14

我們現在約束變數，以使先前的解決方案變得不可行。具體來說，我們要求 x[1] >= 1.5，而 x[0] 保持不受約束。為此，我們將 bounds 參數指定為 least_squares，形式為 bounds=([-np.inf, 1.5], np.inf)。

我們也提供了分析 Jacobian 矩陣

>>> def jac_rosenbrock(x):
...     return np.array([
...         [-20 * x[0], 10],
...         [-1, 0]])

將這些整合在一起，我們可以看到新的解位於邊界上

>>> res_2 = least_squares(fun_rosenbrock, x0_rosenbrock, jac_rosenbrock,
...                       bounds=([-np.inf, 1.5], np.inf))
>>> res_2.x
array([ 1.22437075,  1.5       ])
>>> res_2.cost
0.025213093946805685
>>> res_2.optimality
1.5885401433157753e-07

現在，我們針對具有 100000 個變數的 Broyden 三對角向量值函數，求解方程組（即，成本函數在最小值處應為零）

>>> def fun_broyden(x):
...     f = (3 - x) * x + 1
...     f[1:] -= x[:-1]
...     f[:-1] -= 2 * x[1:]
...     return f

對應的 Jacobian 矩陣是稀疏的。我們告知演算法通過有限差分來估計它，並提供 Jacobian 的稀疏結構，以顯著加速此過程。

>>> from scipy.sparse import lil_matrix
>>> def sparsity_broyden(n):
...     sparsity = lil_matrix((n, n), dtype=int)
...     i = np.arange(n)
...     sparsity[i, i] = 1
...     i = np.arange(1, n)
...     sparsity[i, i - 1] = 1
...     i = np.arange(n - 1)
...     sparsity[i, i + 1] = 1
...     return sparsity
...
>>> n = 100000
>>> x0_broyden = -np.ones(n)
...
>>> res_3 = least_squares(fun_broyden, x0_broyden,
...                       jac_sparsity=sparsity_broyden(n))
>>> res_3.cost
4.5687069299604613e-23
>>> res_3.optimality
1.1650454296851518e-11

我們也來解決一個曲線擬合問題，使用穩健損失函數來處理資料中的離群值。將模型函數定義為 y = a + b * exp(c * t)，其中 t 是預測變數，y 是觀測值，而 a、b、c 是要估計的參數。

首先，定義產生帶有雜訊和離群值資料的函數，定義模型參數，並產生資料

>>> from numpy.random import default_rng
>>> rng = default_rng()
>>> def gen_data(t, a, b, c, noise=0., n_outliers=0, seed=None):
...     rng = default_rng(seed)
...
...     y = a + b * np.exp(t * c)
...
...     error = noise * rng.standard_normal(t.size)
...     outliers = rng.integers(0, t.size, n_outliers)
...     error[outliers] *= 10
...
...     return y + error
...
>>> a = 0.5
>>> b = 2.0
>>> c = -1
>>> t_min = 0
>>> t_max = 10
>>> n_points = 15
...
>>> t_train = np.linspace(t_min, t_max, n_points)
>>> y_train = gen_data(t_train, a, b, c, noise=0.1, n_outliers=3)

定義用於計算殘差和參數初始估計的函數。

>>> def fun(x, t, y):
...     return x[0] + x[1] * np.exp(x[2] * t) - y
...
>>> x0 = np.array([1.0, 1.0, 0.0])

計算標準最小平方法解

>>> res_lsq = least_squares(fun, x0, args=(t_train, y_train))

現在，使用兩個不同的穩健損失函數計算兩個解。參數 f_scale 設定為 0.1，表示內群殘差不應顯著超過 0.1（使用的雜訊水平）。

>>> res_soft_l1 = least_squares(fun, x0, loss='soft_l1', f_scale=0.1,
...                             args=(t_train, y_train))
>>> res_log = least_squares(fun, x0, loss='cauchy', f_scale=0.1,
...                         args=(t_train, y_train))

最後，繪製所有曲線。我們看到，通過選擇適當的 loss，即使在存在強烈離群值的情況下，我們也可以獲得接近最佳的估計值。但請記住，一般建議首先嘗試 ‘soft_l1’ 或 ‘huber’ 損失（如果必要的話），因為其他兩個選項可能會在最佳化過程中引起困難。

>>> t_test = np.linspace(t_min, t_max, n_points * 10)
>>> y_true = gen_data(t_test, a, b, c)
>>> y_lsq = gen_data(t_test, *res_lsq.x)
>>> y_soft_l1 = gen_data(t_test, *res_soft_l1.x)
>>> y_log = gen_data(t_test, *res_log.x)
...
>>> import matplotlib.pyplot as plt
>>> plt.plot(t_train, y_train, 'o')
>>> plt.plot(t_test, y_true, 'k', linewidth=2, label='true')
>>> plt.plot(t_test, y_lsq, label='linear loss')
>>> plt.plot(t_test, y_soft_l1, label='soft_l1 loss')
>>> plt.plot(t_test, y_log, label='cauchy loss')
>>> plt.xlabel("t")
>>> plt.ylabel("y")
>>> plt.legend()
>>> plt.show()
../../_images/scipy-optimize-least_squares-1_00_00.png

在下一個範例中，我們展示如何使用 least_squares() 來最佳化複數變數的複數值殘差函數。考慮以下函數

>>> def f(z):
...     return z - (0.5 + 0.5j)

我們將其封裝到一個返回實數殘差的實數變數函數中，方法是簡單地將實部和虛部作為獨立變數處理

>>> def f_wrap(x):
...     fx = f(x[0] + 1j*x[1])
...     return np.array([fx.real, fx.imag])

因此，我們最佳化的是 2n 個實數變數的 2m 維實數函數，而不是原始的 n 個複數變數的 m 維複數函數

>>> from scipy.optimize import least_squares
>>> res_wrapped = least_squares(f_wrap, (0.1, 0.1), bounds=([0, 0], [1, 1]))
>>> z = res_wrapped.x[0] + res_wrapped.x[1]*1j
>>> z
(0.49999999999925893+0.49999999999925893j)
在此頁面上