優化是機器學習領域最有趣的主題之一。我們日常生活中遇到的大多數問題都是通過數值優化方法解決的。在這篇文章中，讓我們研究一些基本的數值優化演算法，以找到任意給定函式(這對於凸函式最有效)的區域性最優解。讓我們從簡單的凸函式開始，其中區域性和全域性最小值是相同的，然後轉向具有多個區域性和全域性最小值的高度非線性函式。
整個優化圍繞線性代數和微積分的基本概念展開。最近的深度學習更新引起了數值和隨機優化演算法領域的巨大興趣，為深度學習網路所展示的驚人定性結果提供了理論支援。在這些型別的學習演算法中，沒有任何明確已知的優化函式，但我們只能訪問0階和1階的Oracles。Oracles是在任何給定點返回函式值（0階），梯度（1階）或Hessian（2階）的黑盒子。本本提供了對無約束和約束優化函式的基本理論和數值理解，還包括它們的python實現。
一個點成為區域性最小值的必要和充分條件：
設f(.)是一個連續的二階可微函式。對於任一點為極小值，它應滿足以下條件:
- 一階必要條件:
- 二階必要條件：
- 二階充足條件：
梯度下降：
梯度下降是學習演算法(機器學習、深度學習或深度強化學習)領域所有進展的支柱。在這一節中，我們將看到為了更快更好的收斂，梯度下降的各種修改。讓我們考慮線性迴歸的情況，我們估計方程的係數:
假設該函式是所有特徵的線性組合。通過最小化損失函式來確定最佳係數集：
這是線性迴歸任務的最大似然估計。最小二乘法(Ordinary Least Square)涉及到求特徵矩陣的逆，其表示式為:
對於實際問題，資料的維數很大，很容易造成計算量的激增。例如,讓我們考慮問題的影像特徵分析:一般影像的大小1024 x1024,這意味著特徵的數量級為10⁶。由於具有大量的特徵，這類優化問題只能通過迭代的方式來解決，這就導致了我們所熟知的梯度下降法和牛頓-拉弗森法。
梯度下降演算法：
梯度下降演算法利用先前指定的學習率（eta）在負梯度的方向（最陡的下降方向）上更新迭代（X）。學習率用於在任何給定的迭代中防止區域性最小值的overshoot。
下圖顯示了函式f（x）=x²的梯度下降演算法的收斂性，其中eta = 0.25
找出一個最優eta是目前的任務，這需要事先了解函式的理解和操作域。
import matplotlib.pyplot as plt import numpy as np # assuming function to be x**2 def get_gradient(x): return 2*x def get_value(x): return np.sum(x**2) # python implementation of vanilla gradient descent update ruledef def gradient_descent_update (x, eta): """ get_gradient is 1st order oracle """ return x - eta*get_gradient(x)
Armijo Goldstein條件梯度下降:
為了減少手動設定的工作，應用Armijo Goldstein (AG)條件來查詢下一個迭代的(eta)。AG條件的形式推導需要線性逼近、Lipchitz條件和基本微積分的知識。
我們定義兩個函式f1(x)和f2(x)作為兩個不同係數和的f(x)的線性逼近，其具有兩個不同的係數α和β，由下式給出：
在AG條件的每次迭代中，滿足以下關係的特定eta：
找到並且相應地更新當前迭代。
下圖顯示了下一次迭代的範圍，對於函式f（x）=x²的收斂，alpha = 0.25，beta = 0.5：
上圖中紅色、藍色和綠色的線對應於綠色顯示的下一個可能迭代的範圍。
# python implementation of gradient descent with AG condition update rule def gradient_descent_update_AG(x, alpha =0.5, beta =0.25): eta =0.5 max_eta =np.inf min_eta =0. value = get_value(x) grad = get_gradient(x) while True : x_cand = x - (eta)*grad f = get_value(x_cand) f1 = value - eta *a lpha*np.sum(np.abs(grad)* *2 ) f2 = value - eta *be ta *np.sum(np.abs(grad)* *2 ) if f<=f2 and f>=f1: return x_cand if f <= f1: if eta == max_eta: eta = np.min([2.*eta, eta + max_eta/2.]) else: eta = np.min([2.*eta, (eta + max_eta)/2.]) if f>=f2: max_eta = eta eta = (eta+min_eta)/2.0
完全鬆弛條件的梯度下降：
在完全鬆弛條件的情況下，新的函式g（eta）被最小化以獲得隨後用於尋找下一次迭代的eta。
此方法涉及解決查詢每個下一個迭代的優化問題，其執行方式如下：
# The python implementation for FR is shown below: def gradient_descent_update_FR (x): eta = 0.5 thresh = 1e-6 max_eta = np.inf min_eta = 0. while True : x_cand = x - eta*get_gradient(x) g_diff = -1. *get_gradient(x)*get_gradient(x_cand) if np.sum(np.abs(g_diff)) < thresh and eta > 0 : return x_cand if g_diff > 0: if eta == max_eta: eta = np.min([2.*eta, eta + max_eta/2.]) else: eta = np.min([2.*eta, (eta + max_eta)/2.]) else: max_eta = eta eta = (min_eta + eta)/2.0
隨機梯度下降
我們知道，在實際設定中，資料的維數會非常大，這使得對所有特徵進行進一步的梯度計算非常昂貴。在這種情況下，隨機選擇一批點(特徵)並計算期望的梯度。整個過程的收斂只是在預期的意義上。
在數學上它意味著：隨機選擇一個點（p）來估計梯度。
在上面的迭代中，wt可以被視為噪聲。只有當E（wt）趨於0時，該迭代才會導致區域性最優。
同樣，可以看出wt的方差也是有限的。通過以上證明，保證了SGD的收斂性。
# SGD implementation in python
def SGD (self, X, Y, batch_size, thresh= 1 ):
loss = 100
step = 0
if self.plot: losses = []
while loss >= thresh:
# mini_batch formation
index = np.random.randint( 0 , len(X), size = batch_size)
trainX, trainY = np.array(X)[index], np.array(Y)[index]
self.forward(trainX)
loss = self.loss(trainY)
self.gradients(trainY)
# update
self.w0 -= np.squeeze(self.alpha*self.grad_w0)
self.weights -= np.squeeze(self.alpha*self.grad_w)
if self.plot: losses.append(loss)
if step % 1000 == 999 :
print "Batch number: {}" .format(step)+ " current loss: {}" .format(loss)
step += 1
if self.plot : self.visualization(X, Y, losses)
pass
AdaGrad
Adagrad是一個優化器，可幫助自動調整優化問題中涉及的每個特徵的學習率。這是通過跟蹤所有梯度的歷史來實現的。該方法也僅在期望意義上收斂。
def AdaGrad (self, X, Y, batch_size,thresh=0.5 ,epsilon=1e-6 ):
loss = 100
step = 0
if self.plot: losses = []
G = np.zeros((X.shape[ 1 ], X.shape[ 1 ]))
G0 = 0
while loss >= thresh:
# mini_batch formation
index = np.random.randint( 0 , len(X), size = batch_size)
trainX, trainY = np.array(X)[index], np.array(Y)[index]
self.forward(trainX)
loss = self.loss(trainY)
self.gradients(trainY)
G += self.grad_w.T*self.grad_w
G0 += self.grad_w0** 2
den = np.sqrt(np.diag(G)+epsilon)
delta_w = self.alpha*self.grad_w / den
delta_w0 = self.alpha*self.grad_w0 / np.sqrt(G0 + epsilon)
# update parameters
self.w0 -= np.squeeze(delta_w0)
self.weights -= np.squeeze(delta_w)
if self.plot: losses.append(loss)
if step % 500 == 0 :
print "Batch number: {}".format (step)+ " current loss: {}".format(loss)
step += 1
if self.plot : self.visualization(X, Y, losses)
pass
讓我們轉向基於Hessian的方法，牛頓和擬牛頓方法：
基於hessian的方法是基於梯度的二階優化方法，幾何上涉及到梯度和曲率資訊來更新權重，因此收斂速度比僅基於梯度的方法快得多，牛頓方法的更新規則定義為:
該演算法的收斂速度遠快於基於梯度的方法。數學上,梯度下降法的收斂速度正比於O (1 / t),而對於牛頓法它正比於O (1 / t²)。但是對於高維資料，為每個迭代估計二階oracle的計算成本很高，這導致使用一階oracle模擬二階oracle。這就給出了擬牛頓演算法。這類擬牛頓法中最常用的演算法是BFGS和LMFGS演算法。在這一節中，我們只討論BFGS演算法，它涉及到對函式Hessian的rank one矩陣更新。該方法的總體思想是隨機初始化Hessian，並使用rank one更新規則在每次迭代中不斷更新Hessian。數學上可以表示為:
# python implementation for BFGS def BFGS_update (H, s, y): smooth = 1e-3 s = np.expand_dims(s, axis= -1 ) y = np.expand_dims(y, axis= -1 ) rho = 1. /(s.T.dot(y) + smooth) p = (np.eye(H.shape[ 0 ]) - rho*s.dot(y.T)) return p.dot(H.dot(p.T)) + rho*s.dot(s.T) def quasi_Newton (x0, H0, num_iter= 100 , eta= 1 ): xo, H = x0, H0 for _ in range(num_iter): xn = xo - eta*H.dot(get_gradient(xo)) y = get_gradient(xn) - get_gradient(xo) s = xn - xo H = BFGS_update(H, s, y) xo = xn return xo
約束優化
現在，是時候討論一些圍繞約束優化(包括問題的制定和解決策略)的關鍵概念了。本節還將討論一種稱為SVM(支援向量機)的演算法的理論和Python實現。當我們在現實生活中遇到問題時，提出一個理想的優化函式是相當困難的，有時是不可行的，所以我們通過對問題施加額外的約束來放鬆優化函式，優化這個約束設定將提供一個近似，我們將盡可能接近問題的實際解決方案，但也是可行的。求解約束優化問題的方法有拉格朗日公式法、懲罰法、投影梯度下降法、內點法等。在這一節中，我們將學習拉格朗日公式和投影梯度下降法。本節還將詳細介紹用於優化CVXOPT的開源工具箱，並介紹使用此工具箱的SVM實現。
約束優化問題的一般形式：
其中f(x)為目標函式，g(x)、h(x)分別為不等式約束和等式約束。如果f(x)是凸的約束條件形成一個凸集，(即g(x)為凸函式h(x)為仿射函式），該優化演算法保證收斂於全域性最小值。對於其他問題，它收斂於區域性極小值。
投影梯度下降
求解約束優化設定的第一步(也是最明顯的一步)是對約束集進行迭代投影。這是求解約束優化問題中最強大的演算法。這包括兩個步驟(1)在最小化(下降)方向上找到下一個可能的迭代，(2)在約束集上找到迭代的投影。
# projected gradient descent implementation def projection_oracle_l2 (w, l2_norm): return (l2_norm/np.linalg.norm(w))*w def projection_oracle_l1 (w, l1_norm): # first remember signs and store them. Modify w signs = np.sign(w) w = w*signs # project this modified w onto the simplex in first orthant. d=len(w) if np.sum(w)<=l1_norm: return w*signs for i in range(d): w_next = w+ 0 w_next[w> 1e-7 ] = w[w> 1e-7 ] - np.min(w[w> 1e-7 ]) if np.sum(w_next)<=l1_norm: w = ((l1_norm - np.sum(w_next))*w + (np.sum(w) - l1_norm)*w_next)/(np.sum(w)-np.sum(w_next)) return w*signs else : w=w_next def main (): eta= 1.0 /smoothness_coeff for l1_norm in np.arange( 0 , 10 , 0.5 ): w=np.random.randn(data_dim) for i in range( 1000 ): w = w - eta*get_gradient(w) w = projection_oracle_l1(w, l1_norm) pass
瞭解拉格朗日公式
在大多數優化問題中，找到約束集上迭代的投影是一個困難的問題(特別是在複雜約束集的情況下)。它類似於在每次迭代中解決優化問題，在大多數情況下，優化問題是非凸的。在現實中，人們試圖通過解決對偶問題而不是原始問題來擺脫約束。
在深入拉格朗日對偶和原始公式之前，讓我們先了解一下KKT條件及其意義
- 對於任意點為具有等式約束的區域性/全域性最小值:
- 類似地，不等式約束:
這兩個條件可以很容易地通過將單位圓看作一個約束集來觀察。在第一部分中，我們只考慮一個(mu)符號不重要的邊界，這是等式約束的結果。在第二種情況下，考慮一個單位圓的內部集合，其中-ve符號表示(lambda)，表示可行解區域。
KKT (Karush-Kuhn-Tucker)條件被認為是一階必要條件，當一個點被認為是一個平穩點(區域性極小點、區域性極大點、鞍點)時，該條件必須滿足。x ^ * 是區域性最小值：
LICQ條件:所有活動約束
它們應該是線性無關的。
拉格朗日函式
對於任何函式f (x)與不等式約束g_i (x)≤0和等式約束h_i (x) = 0,拉格朗日L (x λμ)
優化函式：
以上優化相當於：
上面的公式被遊戲理論的主張稱為原始問題（p ^ *）（即第二個玩家將總是有更好的優化機會），可以很容易地看出：
這個公式被稱為對偶問題(d ^ *)。
當且僅當目標函式為凸且約束集為凸時，原始公式和對偶公式的最優值相同。這項要求的證明理由如下:
函式是凸的，這意味著
SVM（支援向量機）
SVM屬於用於分類和迴歸問題的監督機器學習演算法類。SVM易於擴充套件，可以解決線性和非線性問題（通過使用基於核的技巧）。在大多數問題中，我們無法對兩類不同的資料進行單獨的區分，因此我們需要在決策邊界的構建中提供一點餘地，這很容易用SVM來表示。支援向量機的思想是在兩組不同的資料點之間建立分離的超平面。一旦獲得了分離超平面，對其中一個類中的資料點(測試用例)進行分類就變得相對容易了。支援向量機甚至可以很好地用於高維資料。與其他機器學習（ML）模型相比，svm的優點是記憶體效率高、準確、快速。我們來看看支援向量機的數學
SVM原始問題
使用拉格朗日演算法的SVM對偶問題
Derivation of Dual
CVXOPT
在本節中，我們將討論使用CVXOPT庫在python中實現上述SVM對偶演算法。
CVXOPT通用格式的問題
# SVM using CVXOPT
import numpy as np
from cvxopt import matrix,solvers
def solve_SVM_dual_CVXOPT (x_train, y_train, x_test):
"""
Solves the SVM training optimisation problem (the Arguments:
x_train: A numpy array with shape (n,d), denoting R^d.
y_train: numpy array with shape (n,) Each element
x_train: A numpy array with shape (m,d), denoting dual) using cvxopt.
n training samples in takes +1 or -1
m test samples in R^d.
Limits:
n<200, d<100, m<1000
Returns:
y_pred_test : A numpy array with shape (m,). This is the result of running the learned model on the
test instances x_test. Each element is +1 or -1.
"""
n, d = x_train.shape
c = 10 # let max lambda value be 10
y_train = np.expand_dims(y_train, axis= 1 )* 1.
P = (y_train * x_train).dot((y_train * x_train).T)
q = -1. *np.ones((n, 1 ))
G = np.vstack((np.eye(n)* -1 ,np.eye(n)))
h = np.hstack((np.zeros(n), np.ones(n) * c))
A = y_train.reshape( 1 , -1 )
b = np.array([ 0.0 ])
P = matrix(P); q = matrix(q)
G = matrix(G); h = matrix(h)
A = matrix(A); b = matrix(b)
sol = solvers.qp(P, q, G, h, A, b)
lambdas = np.array(sol[ 'x' ])
w = ((y_train * lambdas).T.dot(x_train)).reshape( -1 , 1 )
b = y_train - np.dot(x_train, w)
prediction = lambda x, w, b: np.sign(np.sum(w.T.dot(x) + b))
y_test = np.array([prediction(x_, w, b) for x_ in x_test])
return y_test
if __name__ == "__main__" :
# Example format of input to the functions
n= 100
m= 100
d= 10
x_train = np.random.randn(n,d)
x_test = np.random.randn(m,d)
w = np.random.randn(d)
y_train = np.sign(np.dot(x_train, w))
y_test = np.sign(np.dot(x_test, w))
y_pred_test = solve_SVM_dual_CVXOPT(x_train, y_train, x_test)
check1 = np.sum(y_pred_test == y_test)
print ( "Score: {}/{}" .format(check1, len(y_Test)))
