注：該文是上了開智學堂數據科學基礎班的課后做的筆記，主講人是肖凱老師。
最優化
為什么要做最優化呢？因為在生活中，人們總是希望幸福值或其它達到一個極值，比如做生意時希望成本最小，收入最大，所以在很多商業情境中，都會遇到求極值的情況。
函數求根
這里「函數的根」也稱「方程的根」，或「函數的零點」。
先把我們需要的包加載進來。
import numpy as np
import scipy as sp
import scipy.optimize as opt
import matplotlib.pyplot as plt
%matplotlib inline
函數求根和最優化的關系？什么時候函數是最小值或最大值？
兩個問題一起回答：最優化就是求函數的最小值或最大值，同時也是極值，在求一個函數最小值或最大值時，它所在的位置肯定是導數為 0 的位置，所以要求一個函數的極值，必然要先求導，使其為 0，所以函數求根就是為了得到最大值最小值。
scipy.optimize 有什么方法可以求根？
可以用 scipy.optimize 中的 bisect 或 brentq 求根。
f = lambda x: np.cos(x) - x # 定義一個匿名函數
x = np.linspace(-5, 5, 1000) # 先生成 1000 個 x
y = f(x) # 對應生成 1000 個 f(x)
plt.plot(x, y); # 看一下這個函數長什么樣子
plt.axhline(0, color='k'); # 畫一根橫線，位置在 y=0
opt.bisect(f, -5, 5) # 求取函數的根
0.7390851332155535
plt.plot(x, y)
plt.axhline(0, color='k')
plt.scatter([_], [0], c='r', s=100); # 這里的 [_] 表示上一個 Cell 中的結果，這里是 x 軸上的位置，0 是 y 上的位置
求根有兩種方法，除了上面介紹的 bisect，還有 brentq，后者比前者快很多。
%timeit opt.bisect(f, -5, 5)
%timeit opt.brentq(f, -5, 5)
10000 loops, best of 3: 157 µs per loop
The slowest run took 11.65 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 35.9 µs per loop
函數求最小化
求最小值就是一個最優化問題。求最大值時只需對函數做一個轉換，比如加一個負號，或者取倒數，就可轉成求最小值問題。所以兩者是同一問題。
初始值對最優化的影響是什么？
舉例來說，先定義個函數。
f = lambda x: 1-np.sin(x)/x
x = np.linspace(-20., 20., 1000)
y = f(x)
當初始值為 3 值，使用 minimize 函數找到最小值。minimize 函數是在新版的 scipy 里，取代了以前的很多最優化函數，是個通用的接口，背后是很多方法在支撐。
x0 = 3
xmin = opt.minimize(f, x0).x # x0 是起始點，起始點最好離真正的最小值點不要太遠
plt.plot(x, y)
plt.scatter(x0, f(x0), marker='o', s=300); # 起始點畫出來，用圓圈表示
plt.scatter(xmin, f(xmin), marker='v', s=300); # 最小值點畫出來，用三角表示
plt.xlim(-20, 20);
初始值為 3 時，成功找到最小值。
現在來看看初始值為 10 時，找到的最小值點。
x0 = 10
xmin = opt.minimize(f, x0).x
plt.plot(x, y)
plt.scatter(x0, f(x0), marker='o', s=300)
plt.scatter(xmin, f(xmin), marker='v', s=300)
plt.xlim(-20, 20);
由上圖可見，當初始值為 10 時，函數找到的是局部最小值點，可見 minimize 的默認算法對起始點的依賴性。
那么怎么才能不管初始值在哪個位置，都能找到全局最小值點呢？
如何找到全局最優點？
可以使用 basinhopping 函數找到全局最優點，相關背后算法，可以看幫助文件，有提供論文的索引和出處。
我們設初始值為 10 看是否能找到全局最小值點。
x0 = 10
from scipy.optimize import basinhopping
xmin = basinhopping(f,x0,stepsize = 5).x
plt.plot(x, y);
plt.scatter(x0, f(x0), marker='o', s=300);
plt.scatter(xmin, f(xmin), marker='v', s=300);
plt.xlim(-20, 20);
當起始點在比較遠的位置，依然成功找到了全局最小值點。
如何求多元函數最小值？
以二元函數為例，使用 minimize 求對應的最小值。
def g(X):
x,y = X
return (x-1)**4 + 5 * (y-1)**2 - 2*x*y
X_opt = opt.minimize(g, (8, 3)).x # (8,3) 是起始點
print X_opt
[ 1.88292611 1.37658521]
fig, ax = plt.subplots(figsize=(6, 4)) # 定義畫布和圖形
x_ = y_ = np.linspace(-1, 4, 100)
X, Y = np.meshgrid(x_, y_)
c = ax.contour(X, Y, g((X, Y)), 50) # 等高線圖
ax.plot(X_opt[0], X_opt[1], 'r*', markersize=15) # 最小點的位置是個元組
ax.set_xlabel(r"$x_1$", fontsize=18)
ax.set_ylabel(r"$x_2$", fontsize=18)
plt.colorbar(c, ax=ax) # colorbar 表示顏色越深，高度越高
fig.tight_layout()
畫 3D 圖。
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
fig = plt.figure()
ax = fig.gca(projection='3d')
x_ = y_ = np.linspace(-1, 4, 100)
X, Y = np.meshgrid(x_, y_)
surf = ax.plot_surface(X, Y, g((X,Y)), rstride=1, cstride=1, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
cset = ax.contour(X, Y, g((X,Y)), zdir='z',offset=-5, cmap=cm.coolwarm)
fig.colorbar(surf, shrink=0.5, aspect=5);
曲線擬合
曲線擬合和最優化有什么關系？
曲線擬合的問題是，給定一組數據，它可能是沿着一條線散布的，這時要找到一條最優的曲線來擬合這些數據，也就是要找到最好的線來代表這些點，這里的最優是指這些點和線之間的距離是最小的，這就是為什么要用最優化問題來解決曲線擬合問題。
舉例說明，給一些點，找到一條線，來擬合這些點。
先給定一些點：
N = 50 # 點的個數
m_true = 2 # 斜率
b_true = -1 # 截距
dy = 2.0 # 誤差
np.random.seed(0)
xdata = 10 * np.random.random(N) # 50 個 x，服從均勻分布
ydata = np.random.normal(b_true + m_true * xdata, dy) # dy 是標准差
plt.errorbar(xdata, ydata, dy,
fmt='.k', ecolor='lightgray');
上面的點整體上呈現一個線性關系，要找到一條斜線來代表這些點，這就是經典的一元線性回歸。目標就是找到最好的線，使點和線的距離最短。要優化的函數是點和線之間的距離，使其最小。點是確定的，而線是可變的，線是由參數值，斜率和截距決定的，這里就是要通過優化距離找到最優的斜率和截距。
點和線的距離定義如下：
def chi2(theta, x, y):
return np.sum(((y - theta[0] - theta[1] * x)) ** 2)
上式就是誤差平方和。
誤差平方和是什么？有什么作用？
誤差平方和公式為：
誤差平方和大，表示真實的點和預測的線之間距離太遠，說明擬合得不好，最好的線，應該是使誤差平方和最小，即最優的擬合線，這里是條直線。
誤差平方和就是要最小化的目標函數。
找到最優的函數，即斜率和截距。
theta_guess = [0, 1] # 初始值
theta_best = opt.minimize(chi2, theta_guess, args=(xdata, ydata)).x
print(theta_best)
[-1.01442005 1.93854656]
上面兩個輸出即是預測的直線斜率和截距，我們是根據點來反推直線的斜率和截距，那么真實的斜率和截距是多少呢？-1 和 2，很接近了，差的一點是因為有噪音的引入。
xfit = np.linspace(0, 10)
yfit = theta_best[0] + theta_best[1] * xfit
plt.errorbar(xdata, ydata, dy,
fmt='.k', ecolor='lightgray');
plt.plot(xfit, yfit, '-k');
最小二乘（Least Square）是什么？
上面用的是 minimize 方法，這個問題的目標函數是誤差平方和，這就又有一個特定的解法，即最小二乘。
最小二乘的思想就是要使得觀測點和估計點的距離的平方和達到最小，這里的“二乘”指的是用平方來度量觀測點與估計點的遠近（在古漢語中“平方”稱為“二乘”），“最小”指的是參數的估計值要保證各個觀測點與估計點的距離的平方和達到最小。
關於最小二乘估計的計算，涉及更多的數學知識，這里不想詳述，其一般的過程是用目標函數對各參數求偏導數，並令其等於 0，得到一個線性方程組。具體推導過程可參考斯坦福機器學習講義 第 7 頁。
def deviations(theta, x, y):
return (y - theta[0] - theta[1] * x)
theta_best, ier = opt.leastsq(deviations, theta_guess, args=(xdata, ydata))
print(theta_best)
[-1.01442016 1.93854659]
最小二乘 leastsq 的結果跟 minimize 結果一樣。注意 leastsq 的第一個參數不再是誤差平方和 chi2，而是誤差本身 deviations，即沒有平方，也沒有和。
yfit = theta_best[0] + theta_best[1] * xfit
plt.errorbar(xdata, ydata, dy,
fmt='.k', ecolor='lightgray');
plt.plot(xfit, yfit, '-k');
非線性最小二乘
上面是給一些點，擬合一條直線，擬合一條曲線也是一樣的。
def f(x, beta0, beta1, beta2): # 首先定義一個非線性函數，有 3 個參數
return beta0 + beta1 * np.exp(-beta2 * x**2)
beta = (0.25, 0.75, 0.5) # 先猜 3 個 beta
xdata = np.linspace(0, 5, 50)
y = f(xdata, *beta)
ydata = y + 0.05 * np.random.randn(len(xdata)) # 給 y 加噪音
def g(beta):
return ydata - f(xdata, *beta) # 真實 y 和 預測值的差，求最優曲線時要用到
beta_start = (1, 1, 1)
beta_opt, beta_cov = opt.leastsq(g, beta_start)
print beta_opt # 求到的 3 個最優的 beta 值
[ 0.25525709 0.74270226 0.54966466]
拿估計的 beta_opt 值跟真實的 beta = (0.25, 0.75, 0.5) 值比較，差不多。
fig, ax = plt.subplots()
ax.scatter(xdata, ydata) # 畫點
ax.plot(xdata, y, 'r', lw=2) # 真實值的線
ax.plot(xdata, f(xdata, *beta_opt), 'b', lw=2) # 擬合的線
ax.set_xlim(0, 5)
ax.set_xlabel(r"$x$", fontsize=18)
ax.set_ylabel(r"$f(x, \beta)$", fontsize=18)
fig.tight_layout()
除了使用最小二乘，還可以使用曲線擬合的方法，得到的結果是一樣的。
beta_opt, beta_cov = opt.curve_fit(f, xdata, ydata)
print beta_opt
[ 0.25525709 0.74270226 0.54966466]
有約束的最小化
有約束的最小化是指，要求函數最小化之外，還要滿足約束條件，舉例說明。
邊界約束
def f(X):
x, y = X
return (x-1)**2 + (y-1)**2 # 這是一個碗狀的函數
x_opt = opt.minimize(f, (0, 0), method='BFGS').x # 無約束最優化
假設有約束條件，x 和 y 要在一定的范圍內，如 x 在 2 到 3 之間，y 在 0 和 2 之間。
bnd_x1, bnd_x2 = (2, 3), (0, 2) # 對自變量的約束
x_cons_opt = opt.minimize(f, np.array([0, 0]), method='L-BFGS-B', bounds=[bnd_x1, bnd_x2]).x # bounds 矩形約束
fig, ax = plt.subplots(figsize=(6, 4))
x_ = y_ = np.linspace(-1, 3, 100)
X, Y = np.meshgrid(x_, y_)
c = ax.contour(X, Y, f((X,Y)), 50)
ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15) # 沒有約束下的最小值，藍色五角星
ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15) # 有約束下的最小值，紅色星星
bound_rect = plt.Rectangle((bnd_x1[0], bnd_x2[0]),
bnd_x1[1] - bnd_x1[0], bnd_x2[1] - bnd_x2[0],
facecolor="grey")
ax.add_patch(bound_rect)
ax.set_xlabel(r"$x_1$", fontsize=18)
ax.set_ylabel(r"$x_2$", fontsize=18)
plt.colorbar(c, ax=ax)
fig.tight_layout()
不等式約束
介紹下相關理論，先來看下存在等式約束的極值問題求法，比如下面的優化問題。
目標函數是 f(w)，下面是等式約束，通常解法是引入拉格朗日算子，這里使用 \(\beta\) 來表示算子，得到拉格朗日公式為
l 是等式約束的個數。
然后分別對 w 和 \(\beta\) 求偏導，使得偏導數等於 0，然后解出 w 和 \(\beta_{i}\)，至於為什么引入拉格朗日算子可以求出極值，原因是 f(w) 的 dw 變化方向受其他不等式的約束，dw的變化方向與f(w)的梯度垂直時才能獲得極值，而且在極值處，f(w) 的梯度與其他等式梯度的線性組合平行，因此他們之間存在線性關系。（參考《最優化與KKT條件》）
對於不等式約束的極值問題
常常利用拉格朗日對偶性將原始問題轉換為對偶問題，通過解對偶問題而得到原始問題的解。該方法應用在許多統計學習方法中。有興趣的可以參閱相關資料，這里不再贅述。
def f(X):
return (X[0] - 1)**2 + (X[1] - 1)**2
def g(X):
return X[1] - 1.75 - (X[0] - 0.75)**4
x_opt = opt.minimize(f, (0, 0), method='BFGS').x
constraints = [dict(type='ineq', fun=g)] # 約束采用字典定義，約束方式為不等式約束，邊界用 g 表示
x_cons_opt = opt.minimize(f, (0, 0), method='SLSQP', constraints=constraints).x
fig, ax = plt.subplots(figsize=(6, 4))
x_ = y_ = np.linspace(-1, 3, 100)
X, Y = np.meshgrid(x_, y_)
c = ax.contour(X, Y, f((X, Y)), 50)
ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15) # 藍色星星，沒有約束下的最小值
ax.plot(x_, 1.75 + (x_-0.75)**4, 'k-', markersize=15)
ax.fill_between(x_, 1.75 + (x_-0.75)**4, 3, color="grey")
ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15) # 在區域約束下的最小值
ax.set_ylim(-1, 3)
ax.set_xlabel(r"$x_0$", fontsize=18)
ax.set_ylabel(r"$x_1$", fontsize=18)
plt.colorbar(c, ax=ax)
fig.tight_layout()
scipy.optimize.minimize 中包括了多種最優化算法，每種算法使用范圍不同，詳細參考官方文檔。