在人工智能中,需要对于数据进行学习。而学习过程的步骤之一就是找到使得设定好的损失函数(loss function)达到最小值的参数。优化本质上就是找到一个函数最小值的过程。而为了简化问题,凸优化指的是研究定义于凸集中的凸函数最小化的问题。
我们现在想要找到函数$f_0(x)$的最小值。我们使用一种迭代的算法来找。
我们先随便指定一个$x^k$,我们只要找到一个$x^{k+1}$使得$f_0(x^k)>f_0(x^{k+1})$就可以使用$x^{k+1}$来代替$x^k$进入下一轮迭代。
为了让我们更快的计算,我们希望每次迭代效率更高,也就是使得函数$f_0$减少更多。 $$ x^{k+1} = x^k + \alpha^k \cdot d^k $$ 其中,$x$是需要被优化的变量。$\alpha$为每一步走的长度,因此为标量。d是一个单位矢量,指每一步走的方向。
对于一个函数$f_0(x)$,我们先找一个$x_0$。我们使用下面的方法对每一个$x_k$进行迭代。 $$ x^{k+1} = x^k - c \cdot \nabla f(x^k) $$ 其中$c$是一个常数。
当$x^k$与$x^{k+1}$接近到一定程度时停止迭代。
一个梯度下降法的python实现如下:
import numpy as np
def gradient_descent(f, grad, init_x, lr=0.01, tol=1e-6, max_iter=1000):
"""
梯度下降法求解函数的最小值
:param f: 目标函数
:param grad: 目标函数的梯度
:param init_x: 初始点
:param lr: 学习率
:param tol: 容差
:param max_iter: 最大迭代次数
:return: 最小值点和最小值
"""
x = init_x
for i in range(max_iter):
g = grad(x)
x_new = x - lr * g
if np.linalg.norm(x_new - x) < tol:
break
x = x_new
return x, f(x)
# 例子:求解函数 f(x) = (x-1)^2 + y^2 的最小值
def f(x):
return (x[0]-1) ** 2 + x[1] ** 2
def grad(x):
return np.array([2 * (x[0] - 1), 2 * x[1]])
init_x = np.array([1.0, 1.0])
x_min, f_min = gradient_descent(f, grad, init_x)
print("最小值点为:", x_min)
print("最小值为:", f_min)
也就是说,在上面的例子当中,每一步走的方向是$\nabla f(x^k)$,每一步的步长为向量$c \cdot \nabla f(x^k)$的长度。
从上面的例子可以看出,不同的优化算法本质上不同之处体现在步长选择和方向选择的不同。
但是上述的方法也有一定的缺点。比如每一步的步长可能相对较短。
Amijo Rule仍然选择了$\nabla f(x^k)$作为方向,但是改变了步长的选择方法。
若
则
否则停止
其中$\gamma$决定了下降的显著程度,一般在$[0,0.5]$之间,$\beta$一般在$[0,1]$ 之间。
这种方法也可以称为 Back-tracing Algorithm,因为一般会选一个较大的步长,随后逐渐减小。
python实现如下:
import numpy as np
def gradient_descent(f, grad, init_x, lr=1, c=0.1, rho=0.5, max_iter=1000):
"""
使用 Amijo rule 的梯度下降法求解函数的最小值
:param f: 目标函数
:param grad: 目标函数的梯度
:param init_x: 初始点
:param lr: 初始步长
:param c: Amijo rule 中的参数
:param rho: Amijo rule 中的参数
:param max_iter: 最大迭代次数
:return: 最小值点和最小值
"""
x = init_x
for i in range(max_iter):
g = grad(x)
alpha = lr
while f(x - alpha * g) > f(x) - c * alpha * np.dot(g, g):
alpha *= rho
x_new = x - alpha * g
if np.linalg.norm(x_new - x) < 1e-6:
break
x = x_new
return x, f(x)
# 例子:求解函数 f(x) = x^2 + y^2 的最小值
def f(x):
return (x[0] - 1) ** 2 + x[1] ** 2
def grad(x):
return np.array([2 * (x[0] - 1), 2 * x[1]])
init_x = np.array([1.0, 1.0])
x_min, f_min = gradient_descent(f, grad, init_x)
print("最小值点为:", x_min)
print("最小值为:", f_min)
但不仅我们可以改变选择步长的方法,我们也可以改变方向的办法。
坐标轮换本质上的目的在于,因为每次计算多元函数梯度的成本较高,因此每次只对一个坐标进行优化。
例如,对于函数$f_0(x_0,x_1,x_2,...,x_n)$, 可以在固定一部分坐标的情况下,对一个坐标进行优化,然后轮换。
python实现如下:
import numpy as np
def coordinate_descent(f, grad, x0, max_iter=1000, tol=1e-6):
"""
坐标轮换法求解函数的最小值
:param f: 目标函数
:param grad: 目标函数的梯度
:param x0: 初始点
:param max_iter: 最大迭代次数
:param tol: 收敛阈值
:return: 最小值点和最小值
"""
x = x0.copy()
n = x0.size
for iter in range(max_iter):
for i in range(n):
ei = np.zeros(n)
ei[i] = 1
gi = grad(x).dot(ei)
#这个地方采取了梯度点成方向向量的方法计算一个方向上的梯度,其实没有体现出可以只计算一个方向上偏导的有点。。。
if abs(gi) > 1e-10:
alpha = -f(x) / gi
x[i] += alpha
if np.linalg.norm(grad(x)) < tol:
break
return x, f(x)
# 例子:求解函数 f(x) = (x1-1)^2 + x2^2 + x3^2 的最小值
def f(x):
return (x[0] - 1) ** 2 + x[1] ** 2 + x[2] ** 2
def grad(x):
return np.array([2 * (x[0] - 1), 2 * x[1], 2 * x[2]])
x0 = np.array([1.0, 2.0, 3.0])
x_min, f_min = coordinate_descent(f, grad, x0)
print("最小值点为:", x_min)
print("最小值为:", f_min)
有时候也被叫做Zig-zag
这个方法也可以进行拓展。对于多元函数,不一定需要轮换坐标,而是把多个坐标拆分成若干子空间进行循环优化。
牛顿法是一种对于梯度下降的改进。梯度下降的缺点在于它只关注了一个函数的一阶导数。为了更好的知道这个函数的性质,牛顿法采用的不只是线性模拟,而包括了二次项。
为了让上面的式子取到最小值,
因此牛顿法的流程是: $$ d^k = -(\nabla^2f(x^k) )^{-1}\nabla f(x^k)\ \alpha^k = \arg \min f(x^k+\alpha^k\cdot d^k)\ x^{k+1} = x^k + \alpha^k \cdot d^k $$ 一个python实现的牛顿法求解如下
import numpy as np
def newton_method(f, x0, tol=1e-8, max_iter=1000):
# f是目标函数,x0是初始点,tol是容差,max_iter是最大迭代次数
x = x0
for i in range(max_iter):
# 计算函数值、一阶导数和二阶导数
fx = f(x)
grad_f = np.array([f(x + h) - fx for h in np.eye(len(x))]) / 1e-6
hess_f = np.array([(grad_f(x + h) - grad_f(x)) / 1e-6 for h in np.eye(len(x))])
# 计算搜索方向和步长
dx = np.linalg.solve(hess_f, -grad_f)
alpha = 1.0
while f(x + alpha * dx) > fx:
alpha *= 0.5
# 更新x
x = x + alpha * dx
# 判断是否达到容差
if np.linalg.norm(dx) < tol:
break
return x
# 定义目标函数
def f(x):
return (x[0]-1)**2 + x[1]**2
# 调用牛顿法
x0 = np.array([0.0, 0.0])
x_opt = newton_method(f, x0)
print("优化结果:", x_opt)
print("目标函数最小值:", f(x_opt))
牛顿法的问题主要在于计算Hessian矩阵计算量很大。基于这一点,后来也有不少拟牛顿法。这些方法通过找到一个方便计算的近似Hessian矩阵来实现目的。
一般的优化问题会被一定条件束缚。
因此一般的优化问题有以下的形式:
$$
min \space f_0(x)\ s.t. \space f_i(x) \le 0, \space i = 1,2,...,m\ h_i(x) = 0,i=1,2,...,n
$$
对于这种问题,最简单得方法是采用拉格朗日乘子法
我们先讨论只有等式约束的情况,即 $$ \min \space f_0(x)\ s.t. h_i(x) = 0,i=1,2,...,n $$
对于凸问题,我们直接可以说极值点在一阶导为零处,也就是说 $$ \frac{\partial}{\partial x}(f_0(x)+ \lambda_i \cdot h_i(x)) = 0 $$ 我们发现$h_i(x) = 0$也可以被改成类似的形式 $$ \frac{\partial}{\partial \lambda_i}(f_0(x)+ \lambda_i \cdot h_i(x)) = 0 $$ 因此整个问题可以转化为 $$ min \space g(x,\lambda)\ g(x,\lambda) = f_0(x)+ \lambda_i \cdot h_i(x) $$
一个拉格朗日乘子法求最小值的python实现
import numpy as np
from scipy.optimize import minimize
# 定义目标函数
def f(x):
return x[0]**2 + x[1]**2
# 定义等式约束
def g(x):
return x[0] + x[1] - 1
# 定义拉格朗日函数
def lagrangian(x, lmbda):
return f(x) - lmbda * g(x)
# 定义求解函数
def solve_lagrangian(x0):
# 定义约束条件
cons = {'type': 'eq', 'fun': g}
# 初始拉格朗日乘子为0
lmbda0 = 0.0
# 定义求解函数
def lagrangian_min(x):
return lagrangian(x, lmbda0)
# 调用优化函数求解
res = minimize(lagrangian_min, x0, constraints=cons)
# 提取优化结果和乘子
x_opt = res.x
lmbda_opt = res.fun / g(x_opt)
return x_opt, lmbda_opt
# 求解优化问题
x0 = [0.5, 0.5]
x_opt, lmbda_opt = solve_lagrangian(x0)
# 输出结果
print("优化结果:", x_opt)
print("拉格朗日乘子:", lmbda_opt)
print("目标函数最小值:", f(x_opt))
print("约束条件值:", g(x_opt))
我们继续讨论含有不等式约束的情况,以及继续使用拉格朗日乘子
$$
min \space f_0(x)\ s.t. \space f_i(x) \le 0, \space i = 1,2,...,m\ h_i(x) = 0,i=1,2,...,n
$$
定义
令$g(\lambda,\mu) = inf_x \space L(x,\lambda,\mu) $
那么一定有
因此
以下我们不加证明的给出下面的定理(Slater's Condition):
对于凸问题$min \space f_0(x) \space s.t. \space f_i(x) \le 0 \space (i = 1,2,...,n), \space Ax=b$, 满足
为什么会出现这种情况?下面不加数学的给出一个经济学解释。
所以我们可以知道
从上面的各种情况,我们可以总结出一下达到最小值的条件(Karush-Kuhn-Tucker Conditions)
对于能使得$f_0(x)$在符合约束条件下达到最小值的$x^$,一定符合KKT条件 $$ f_i(x^)\le0\ h_i(x^*)=0\ \lambda_i \ge 0\ \lambda_i \cdot f_i(x) = 0\ \nabla f_0(x) + \Sigma\lambda_i \cdot \nabla f_i(x) + \mu_i \cdot \nabla h_i(x) = 0 $$ (最后一个条件是求导等于零)
Water-filling Problem: $$ min \space \Sigma^{n}_{i=1} \log (\alpha_i + x_i) \ s.t. \space x \ge 0, \space 1^T x = 0 $$ 该问题的KKT条件: $$ x^* \ge 0\ 1^Tx^* = 1\ \lambda_i^* \cdot x_i^=0,\space i=1,2,...,n\ -\frac{1}{\alpha_i + x_i^}-\lambda_i^+v^ = 0 $$ 消掉$\lambda^$,得到 $$ x^\cdot(v^-\frac{1}{\alpha_i+x_i^}) = 0,i=1,2,...,n\ v^* \ge \frac{1}{\alpha_i+x_i^},i=1,2,...,n $$ 也就是 $$ v^ \ge \frac{1}{\alpha_i} \Rightarrow x^_i = 0\ v^ \lt \frac{1}{\alpha_i} \Rightarrow x^_i = \frac{1}{v_i^}-\alpha_i $$
需要解决的问题是: $$ \min \space f_0(x)\ s.t. h_i(x) = 0,i=1,2,...,n $$ 惩罚函数法:
不在纠结于一定要满足等式条件。把等式约束下的优化问题改成一个无约束优化问题 $$ \min \space f_0(x)+\frac{c}{2}\cdot \Sigma^n_{i=1} h_i^2(x) $$ 当参数c趋于无穷时,得到的最优解逐渐趋近于原本的最优解。
所以贪婪的人类必定会把两个好东西放在一起。
相比于传统的拉格朗日函数,我们引入一些很新的东西: $$ L_c(x,v) = f(x) + v^T h(x) + \frac{c}{2} |h(x)|^2 $$ 这个问题和原来的拉格朗日函数是等价的
我们考虑约束函数
原来的拉格朗日法迭代的函数是这样的: $$ x^{k+1}=x^k-\alpha^k \cdot (\nabla f(x^k)+A^Tv^k)\ v^{k+1}=v^k+\alpha^k(Ax^k-b) $$ 现在我们的迭代过程变成了这样: $$ x^{k+1}=\arg \min_x L_c(x,v_k)\ v^{k+1}=v^k+\alpha^k(Ax^k-b) $$ (这里的两个方法只是找到了下降的方向,至于具体的步长需要继续步长搜索,例如Amijo Rule)
import numpy as np
def augmented_lagrangian(f, g, x0, lambd0, mu0, maxiter, tol, rho=1.0, gamma=2.0):
x = x0
lambd = lambd0
mu = mu0
it = 0
converged = False
while not converged and it < maxiter:
# 求解Lagrange乘子
def L(x):
return f(x) + np.dot(lambd, g(x)) + (mu/2.0) * np.linalg.norm(g(x))**2
# 梯度下降
x = gradient_descent(lambda x: L(x), x, tol=tol, rho=rho, gamma=gamma)
# 更新Lagrange乘子
lambd = lambd - mu * g(x)
# 更新惩罚参数
mu = gamma * mu
# 检查是否收敛
if np.linalg.norm(g(x)) < tol:
converged = True
it += 1
return x
def gradient_descent(f, x0, tol, rho, gamma, maxiter=1000):
x = x0
it = 0
converged = False
while not converged and it < maxiter:
gradient = compute_gradient(f, x)
x_new = x - rho * gradient
if np.linalg.norm(x_new - x) < tol:
converged = True
x = x_new
rho *= gamma
it += 1
return x
def compute_gradient(f, x, eps=1e-8):
gradient = np.zeros_like(x)
for i in range(len(x)):
x_plus = x.copy()
x_plus[i] += eps
x_minus = x.copy()
x_minus[i] -= eps
gradient[i] = (f(x_plus) - f(x_minus)) / (2.0 * eps)
return gradient
对于大型人工智能或者计算程序,分布式计算有利于提高效率。同样,优化问题也有办法被拆分进行。例如下面的问题:
$$
\min (f(x) + g(x))
$$
可以转化为
$$
\min f(x) + g(x)\s.t.x-z=0
$$
我们定义增广拉格朗日函数
分成两步解决: $$ 1.[x^{k+1},z^{k+1}]=\arg \min_{x,z} (f(x)+g(z)+v(x-z)+\frac{c}{2}\cdot ||x-z||^2_2)\ 2.v^{k+1}=v^k+c \cdot(x^{k+1}-z^{k+1}) $$ 而对于第一步拆分成以下的步骤。考虑对于大循环第k步的情况 $$ 1A. x^{k|t+1}=\arg \min_x f(x)+\frac{c}{2}\cdot ||x^{k|t}-z^{k|t}+\frac{v^k}{2}||^2_2)\ 1B. z^{k|t+1}=\arg \min_z g(z)+\frac{c}{2}\cdot ||x^{k|t}-z^{k|t}+\frac{v^k}{2}||^2_2) $$ 我们可以把1A和1B拆分到两个小计算机上进行,把步骤2在大计算机上运行,实现优化的分布。
-
Cnvex Optimization Stephen Boyd
-
感谢chatGPT写的python