导购网站如何做免费推广,什么是企业网站建设,自建网站套现,辽宁网站建设前言1#xff1a;这篇博客仅限于介绍拉格朗日乘子法#xff0c;KKT条件#xff0c;ALM算法#xff0c;ADMM算法等最优化方法的使用以及简版代码实现#xff0c;但不会涉及具体的数学推导#xff1b;不过在下面我会给出具体数学推导的相关文章和截图#xff0c;供学有余力…前言1这篇博客仅限于介绍拉格朗日乘子法KKT条件ALM算法ADMM算法等最优化方法的使用以及简版代码实现但不会涉及具体的数学推导不过在下面我会给出具体数学推导的相关文章和截图供学有余力的同志参考
前言2从OSQP求解器的官方网站可以得知OSQP求解器使用的最优化方法即ADMM算法 a
a
a
a
a
a
1. 拉格朗日乘子法与KKT条件
注这一章节不从数学推导上而是从实际数学意义上直观的讲解拉格朗日乘子法和KKT条件是如何推导得来的 1.1 无约束的优化问题 -- 梯度下降法
首先从无约束的优化问题讲起一般就是要使表达式 取到最小值对于这类问题在高中就学过怎么做只要对它的每一个变量求导然后让偏导为零解方程组就行了
在极值点处一定满足 只是必要条件然后对它进行求解再代入验证是否真的是极值点就行了
但是有很多问题通过直接求导是解不出来的或者很难解所以就需要梯度下降法、牛顿法、坐标下降法之类的数值迭代算法对于这些迭代算法就像下面这张图一样我们希望找到其中的最小值一个比较直观的想法是先找一个起点然后不断向最低点靠近。就先把一个小球放到一个碗里一样 一开始要找一个起始点然后确定走的方向和距离最后还要知道什么时候停止这三步中最难的是确定走的方向走的慢点还可以接受要是方向错了就找不到最小值了
走的距离可以简单的设为一个比较小的值起始点可以随机选一个 关键是方向可以选择 处的梯度的反方向这是函数在这个点下降最快的方向最终的停止条件是梯度的大小很接近于 0在极值点处的梯度大小就是 0
这种方法依靠梯度确定下降方向的方法叫做梯度下降法 对 求极小值的流程
随机选定起点 得到函数在 的梯度然后从 向前走一步 重复第2步直到梯度接近于0小于一个事先设定的小值或到达指定的迭代次数上限 1.2 有约束优化问题 -- 拉格朗日乘子法和KKT条件
摘自文章
https://www.cnblogs.com/xinchen1111/p/8804858.html什么是仿射函数-CSDN博客https://www.cnblogs.com/fuleying/p/4481334.html
a
a
1.2.1 单个等式约束
我们可能要在满足一定约束条件的情况下最小化目标函数比如有一个等式约束 如图其中的圆圈是目标函数 () 投影在平面上的等值线黑线是约束条件 ℎ()0 的函数图像等值线与函数图像相交的点其实就是所有满足约束的点
那么极值点只有可能在等值线与函数图像相切的地方取到因为如果在相交的地方取到那么沿着 ℎ() 的图像往前走或者往后走一定还有其它的等值线与它相交也就是 (,) 的值还能变大和变小所以交点不是极值点只有相切的时候才有可能是极值点不可能同时变大和变小往前后两个方向走只能同时变大/变小这才是极值点
在相切的地方 ℎ() 的梯度和 (,) 的梯度应该是在同一条直线上的在切点处两个函数的梯度都与切平面垂直所以在一条直线上
所以满足条件的极值点一定满足∇(,)∇ℎ(,) 和 ℎ(,)0
那么只要解出这个方程组就可以得到问题的解析解了当然也存在解不出来的情况就需要用罚函数法之类的方法求数值解了
为了方便和好记就把原来的优化问题写成 (,)ℎ(,) 的形式然后分别对 ,, 求偏导并且令偏导为 0 就行了对 x,y 的偏导为0可以得到式子 ∇(,)-∇ℎ(,)0 对 的偏导为0可以得到式子 ℎ(,)0 和之前得到的方程组是一样的这种方法叫拉格朗日乘数法
a
a
1.2.2 多个等式约束
将拉格朗日乘子法扩展到含有多个等式约束的情况 这里的平面和球面分别代表了两个约束 ℎ1() 和 ℎ2()那么这个问题的可行域就是它们相交的那个圆这里蓝色箭头表示平面的梯度黑色箭头表示球面的梯度那么相交的圆的梯度就是它们的线性组合所以在极值点处目标函数的梯度和约束的梯度的线性组合在一条直线上所以满足 为了好记将原来的约束的问题写成 的形式然后对 、 求偏导让它们为 0结果像上面一样直接列方程组是一样的这可以看做是一种简记或者是对偶问题这个函数叫做拉格朗日函数
a
a
1.2.3 同时含有等式约束和不等式约束
如果问题中既有等式约束又有不等式约束怎么办呢对于 当然也同样约定不等式是 ≤如果是 ≥ 只要取反就行了
对于这个问题先不看等式约束对于不等式约束和目标函数的图 阴影部分就是可行域也就是说可行域从原来的一条线变成了一块区域那么能取到极值点的地方可能有两种情况
还是在 ℎ() 和 等值线相切的地方() 的极值点本身就在可行域里面
因为如果不是相切那么同样的对任意一个在可行域中的点如果在它附近往里走或者往外走() 一般都会变大或者变小所以绝大部分点都不会是极值点除非这个点刚好在交界处且和等值线相切或者这个点在可行域内部但是本身就是 f(x)() 的极值点 对于第一种情况不等式约束就变成等式约束了对 ()ℎ()() 用拉格朗日乘子法 对于第二种情况不等式约束就相当于没有对 ()ℎ() 用拉格朗日乘子法 把两种情况用同一组方程表示出来对比一下两个问题
第一种情况中有 ≥0 且 ()0第二种情况 0 且 ()≤0
综合两种情况可以写成 ()0 且 ≥0 且 ()≤0 这个就是 KKT 条件它的含义是这个优化问题的极值点一定满足这组方程组不是极值点也可能会满足但是不会存在某个极值点不满足的情况它也是原来的优化问题取得极值的必要条件解出来了极值点之后还是要代入验证的但是因为约束比较多情况比较复杂KKT 条件并不是对于任何情况都是满足的
要满足 KKT 条件需要有一些规范性条件Regularity conditions就是要求约束条件的质量不能太差常见的比如
LCQ如果 ℎ() 和 () 都是形如 的仿射函数那么极值一定满足 KKT 条件LICQ起作用的 () 函数即 () 相当于等式约束的情况和 ℎ() 函数在极值点处的梯度要线性无关那么极值一定满足 KKT 条件Slater 条件如果优化问题是个凸优化问题且至少存在一个点满足 ℎ()0 和 ()0极值一定满足 KKT 条件并且满足强对偶性质 a
a
a
a
a
a
2. ALM 算法
注1ALMAugmented Lagrange Multiplier算法即增广型拉格朗日乘子法
注2ALM 算法常用于线性规划问题 不能有高阶项/根号项也不能有两个变量之间的交乘积项
注3ALM 算法的推导过程 -- http://faculty.bicmr.pku.edu.cn/~wenzw/optbook/lect/16-lect-alm.pdf 2.1 ALM 算法介绍 只考虑含有等式约束的问题 a a 注含有不等式约束问题的 ALM 算法见文章 -- http://faculty.bicmr.pku.edu.cn/~wenzw/optbook/lect/16-lect-alm.pdf 构造增广拉格朗日函数 拉格朗日乘子矩阵惩罚项权重标量 a a 注意增广拉格朗日乘子法就是在拉格朗日乘子法获得的函数后面加上针对所有等式约束的惩罚项 a a 优点 原拉格朗日函数的收敛性不太好抖动很大不稳定加上了惩罚项可以增加原拉格朗日函数的凸性使得我们可以通过更简单的方法如梯度下降法去求解这个函数的最优解加强了等式约束的限制作用针对等式约束增加了惩罚项使得在迭代的过程中迭代点一直是在约束附近的区域进行梯度下降不会跑太远从而加快了求解速度 ALM 算法的求解过程梯度下降 先对变量 进行梯度下降再对拉格朗日乘子 进行梯度下降上述两步为一次迭代不断这样迭代下去直至收敛或者到达指定迭代次数 2.2 ALM 算法代码示例 举例只含等式约束 a
使用 Scipy 中的函数 linprog 求解线性规划问题将求解得到的结果作为参考答案
函数 linprogscipy.optimize.linprog函数参数最全详解_optimize的linprog-CSDN博客 solve the following problem with scipy
min f(x) -3x[0] - 5x[1]
s.t. x[0] x[2] 42x[1] x[3] 123x[0] 2x[1] x[4] 18x[0], x[1], x[2], x[3], x[4] 0optimal solution:
fun: -36.0
x: [ 2.000e00 6.000e00 2.000e00 0.000e00 0.000e00]
from scipy.optimize import linprog
import torchc torch.tensor([-3, -5, 0, 0, 0], dtypetorch.float32)
A torch.tensor([[1, 0, 1, 0, 0], [0, 2, 0, 1, 0], [3, 2, 0, 0, 1]], dtypetorch.float32)
b torch.tensor([4, 12, 18], dtypetorch.float32)res linprog(c, A_eqA, b_eqb, bounds(0, None))
print(res)
a
ALM 算法示例 solve the following problem with Augmented Lagrange Multiplier method
min f(x) -3x[0] - 5x[1]
s.t. x[0] x[2] 42x[1] x[3] 123x[0] 2x[1] x[4] 18x[0], x[1], x[2], x[3], x[4] 0问题形式
min f(x) c^T x
s.t. Ax b
x 0
import torchdef lagrangian_function(x, lambda_):return f(x) lambda_ (A x - b) alpha / 2 * ((A x - b)**2).sum()def f(x):return c xdef update_x(x, lambda_): update x with gradient descent lagrangian_function(x, lambda_).backward()new_x x - eta * x.gradx.data new_x.clamp(min0)x.grad.zero_()def update_lambda(lambda_):new_lambda lambda_ alpha * (A x - b)lambda_.data new_lambdadef pprint(i, x, lambda_):print(f\n{i1}th iter, L:{lagrangian_function(x, lambda_):.2f}, f: {f(x):.2f})print(fx: {x})print(flambda: {lambda_})print(constraints violation: )print(A x - b)def solve(x, lambda_):for i in range(500):pprint(i, x, lambda_) # 更新 xupdate_x(x, lambda_) # 更新拉格朗日乘子update_lambda(lambda_) # 打印信息if __name__ __main__:eta 0.03 # 学习率alpha 1 # 惩罚权重c torch.tensor([-3, -5, 0, 0, 0], dtypetorch.float32)A torch.tensor([[1, 0, 1, 0, 0], [0, 2, 0, 1, 0], [3, 2, 0, 0, 1]], dtypetorch.float32)b torch.tensor([4, 12, 18], dtypetorch.float32)lambda_ torch.tensor([0, 0, 0], dtypetorch.float32)x torch.tensor([2, 0, 2, 0, 0], dtypetorch.float32, requires_gradTrue)solve(x, lambda_) a
a
a
a
a
a
3. ADMM 算法
注1ADMMAlternating Direction Method of Multipliers算法即交替方向乘子法
注2关于详细的数学推导 -- 一文详解从拉格朗日乘子法、KKT条件、对偶上升法到罚函数与增广Lagrangian乘子法再到ADMM算法交替方向乘子法_拉格朗日乘子 二次惩罚系数-CSDN博客 3.1 ADMM 算法介绍 只考虑含有等式约束的问题 a a 和ALM算法的不同点在于将所有的变量分成几部分假设为两部分一部分是向量 一部分是向量 然后分别进行梯度下降而不是一次性将所有的变量全部都进行梯度下降 构造增广拉格朗日函数 拉格朗日乘子矩阵惩罚项权重标量 ADMM 算法的求解过程 先对变量 进行梯度下降再对变量 进行梯度下降再对变量 进行梯度下降上述三步为一次迭代不断这样迭代下去直至收敛或者到达指定迭代次数 3.2 ADMM 算法代码示例 solve the following problem with Alternating Direction Multiplier Method
min f(x) -3x[0] - 5x[1]
s.t. x[0] x[2] 42x[1] x[3] 123x[0] 2x[1] x[4] 18x[0], x[1], x[2], x[3], x[4] 0问题形式
min f(x) c^T x
s.t. Ax b
x 0
import torchdef lagrangian_function(x, lambda_):return f(x) lambda_ (A x - b) alpha / 2 * ((A x - b)**2).sum()def f(x):return c xdef update_x(x, lambda_): update x with gradient descent for i in range(len(x)):# not the best way to calculate gradient# 这里直接将所有的变量拆成了5份这并不是最好的求解方式lagrangian_function(x, lambda_).backward()x_i x[i] - eta * x.grad[i]x.data[i] max(0, x_i)x.grad.zero_()def update_lambda(lambda_):new_lambda lambda_ alpha * (A x - b)lambda_.data new_lambdadef pprint(i, x, lambda_):print(f\n{i1}th iter, L:{lagrangian_function(x, lambda_):.2f}, f: {f(x):.2f})print(fx: {x})print(flambda: {lambda_})print(constraints violation: )print(A x - b)def solve(x, lambda_):for i in range(500):pprint(i, x, lambda_)update_x(x, lambda_)update_lambda(lambda_)if __name__ __main__:eta 0.03 # 学习率alpha 1 # 惩罚权重c torch.tensor([-3, -5, 0, 0, 0], dtypetorch.float32)A torch.tensor([[1, 0, 1, 0, 0], [0, 2, 0, 1, 0], [3, 2, 0, 0, 1]], dtypetorch.float32)b torch.tensor([4, 12, 18], dtypetorch.float32)lambda_ torch.tensor([0, 0, 0], dtypetorch.float32)x torch.tensor([2, 0, 2, 0, 0], dtypetorch.float32, requires_gradTrue)solve(x, lambda_)
