笔记101:OSQP求解器的底层算法 -- ADMM算法

前言1:这篇博客仅限于介绍拉格朗日乘子法,KKT条件,ALM算法,ADMM算法等最优化方法的使用以及简版代码实现,但不会涉及具体的数学推导;不过在下面我会给出具体数学推导的相关文章和截图,供学有余力的同志参考; 

前言2:从OSQP求解器的官方网站可以得知,OSQP求解器使用的最优化方法即ADMM算法;


a

a

a

a

a

a

1. 拉格朗日乘子法与KKT条件

注:这一章节不从数学推导上,而是从实际数学意义上,直观的讲解拉格朗日乘子法和KKT条件是如何推导得来的;


1.1 无约束的优化问题 -- 梯度下降法

首先从无约束的优化问题讲起,一般就是要使表达式 minf(x) 取到最小值;对于这类问题在高中就学过怎么做,只要对它的每一个变量求导,然后让偏导为零,解方程组就行了;

在极值点处一定满足 \frac{df(x)}{dx}=0 (只是必要条件),然后对它进行求解,再代入验证是否真的是极值点就行了;

但是有很多问题通过直接求导是解不出来的,或者很难解,所以就需要梯度下降法、牛顿法、坐标下降法之类的数值迭代算法;对于这些迭代算法就像下面这张图一样,我们希望找到其中的最小值,一个比较直观的想法是先找一个起点,然后不断向最低点靠近。就先把一个小球放到一个碗里一样;

一开始要找一个起始点,然后确定走的方向和距离,最后还要知道什么时候停止;这三步中最难的是确定走的方向,走的慢点还可以接受,要是方向错了就找不到最小值了;

  • 走的距离可以简单的设为一个比较小的值;
  • 起始点可以随机选一个 (x_0,y_0) ;
  • 关键是方向,可以选择 (x_0,y_0) 处的梯度的反方向,这是函数在这个点下降最快的方向;
  • 最终的停止条件是梯度的大小很接近于 0(在极值点处的梯度大小就是 0);

这种方法依靠梯度确定下降方向的方法叫做梯度下降法;

 对 f(x) 求极小值的流程:

  1. 随机选定起点 x_0
  2. 得到函数在 x_0 的梯度,然后从 x_0 向前走一步( x_0^{new}=x_0^{old}-\alpha \bigtriangledown f(x) )
  3. 重复第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 单个等式约束

我们可能要在满足一定约束条件的情况下最小化目标函数,比如有一个等式约束:

  • minf(x)
  • h(x)=0

如图,其中的圆圈是目标函数 𝑓(𝑥,𝑦) 投影在平面上的等值线,黑线是约束条件 ℎ(𝑥)=0 的函数图像;等值线与函数图像相交的点其实就是所有满足约束的点;

那么极值点只有可能在等值线与函数图像相切的地方取到,因为如果在相交的地方取到,那么沿着 ℎ(𝑥) 的图像往前走或者往后走,一定还有其它的等值线与它相交,也就是 𝑓(𝑥,𝑦) 的值还能变大和变小,所以交点不是极值点,只有相切的时候才有可能是极值点(不可能同时变大和变小,往前后两个方向走只能同时变大/变小,这才是极值点);

在相切的地方 ℎ(𝑥) 的梯度和 𝑓(𝑥,𝑦) 的梯度应该是在同一条直线上的(在切点处两个函数的梯度都与切平面垂直,所以在一条直线上);

所以满足条件的极值点一定满足:∇𝑓(𝑥,𝑦)=𝜆∇ℎ(𝑥,𝑦)  和 ℎ(𝑥,𝑦)=0

那么只要解出这个方程组,就可以得到问题的解析解了;当然也存在解不出来的情况,就需要用罚函数法之类的方法求数值解了;

为了方便和好记,就把原来的优化问题写成 𝑓(𝑥,𝑦)+𝜆ℎ(𝑥,𝑦) 的形式,然后分别对 𝜆,𝑥,𝑦 求偏导,并且令偏导为 0 就行了(对 x,y 的偏导为0可以得到式子 ∇𝑓(𝑥,𝑦)-𝜆∇ℎ(𝑥,𝑦)=0 )( 对 𝜆 的偏导为0可以得到式子 ℎ(𝑥,𝑦)=0 ),和之前得到的方程组是一样的;这种方法叫拉格朗日乘数法;

a

a

1.2.2 多个等式约束

将拉格朗日乘子法扩展到含有多个等式约束的情况:

这里的平面和球面分别代表了两个约束 ℎ1(𝑥) 和 ℎ2(𝑥),那么这个问题的可行域就是它们相交的那个圆;这里蓝色箭头表示平面的梯度,黑色箭头表示球面的梯度,那么相交的圆的梯度就是它们的线性组合,所以在极值点处目标函数的梯度和约束的梯度的线性组合在一条直线上,所以满足:

为了好记,将原来的约束的问题写成 L(x,\lambda )=f(x)+\sum_{i=1}^{n}[\lambda _i\bigtriangledown h_i(x)] 的形式,然后对 𝑥、𝜆 求偏导,让它们为 0;结果像上面一样直接列方程组是一样的,这可以看做是一种简记,或者是对偶问题,这个函数叫做拉格朗日函数;

a

a

1.2.3 同时含有等式约束和不等式约束

如果问题中既有等式约束,又有不等式约束怎么办呢?对于:

当然也同样约定不等式是 ≤,如果是 ≥ 只要取反就行了;

对于这个问题先不看等式约束,对于不等式约束和目标函数的图:

阴影部分就是可行域,也就是说可行域从原来的一条线变成了一块区域;那么能取到极值点的地方可能有两种情况:

  1. 还是在 ℎ(𝑥) 和 等值线相切的地方
  2. 𝑓(𝑥) 的极值点本身就在可行域里面

因为如果不是相切,那么同样的,对任意一个在可行域中的点,如果在它附近往里走或者往外走,𝑓(𝑥) 一般都会变大或者变小,所以绝大部分点都不会是极值点;除非这个点刚好在交界处,且和等值线相切;或者这个点在可行域内部,但是本身就是 f(x)𝑓(𝑥) 的极值点;

 对于第一种情况,不等式约束就变成等式约束了,对 𝑓(𝑥)+𝜆ℎ(𝑥)+𝜇𝑔(𝑥) 用拉格朗日乘子法:

对于第二种情况,不等式约束就相当于没有,对 𝑓(𝑥)+𝜆ℎ(𝑥) 用拉格朗日乘子法:

把两种情况用同一组方程表示出来,对比一下两个问题:

  • 第一种情况中有 𝜇≥0 且 𝑔(𝑥)=0
  • 第二种情况 𝜇=0 且 𝑔(𝑥)≤0

综合两种情况,可以写成 𝜇𝑔(𝑥)=0 且 𝜇≥0 且 𝑔(𝑥)≤0:

这个就是 KKT 条件;它的含义是这个优化问题的极值点一定满足这组方程组(不是极值点也可能会满足,但是不会存在某个极值点不满足的情况);它也是原来的优化问题取得极值的必要条件,解出来了极值点之后还是要代入验证的,但是因为约束比较多,情况比较复杂,KKT 条件并不是对于任何情况都是满足的;

要满足 KKT 条件需要有一些规范性条件(Regularity conditions),就是要求约束条件的质量不能太差,常见的比如:

  1. LCQ:如果 ℎ(𝑥) 和 𝑔(𝑥) 都是形如 𝐴𝑥+𝑏 的仿射函数,那么极值一定满足 KKT 条件;
  2. LICQ:起作用的 𝑔(𝑥) 函数(即 𝑔(𝑥) 相当于等式约束的情况)和 ℎ(𝑥) 函数在极值点处的梯度要线性无关,那么极值一定满足 KKT 条件;
  3. Slater 条件:如果优化问题是个凸优化问题,且至少存在一个点满足 ℎ(𝑥)=0 和 𝑔(𝑥)=0,极值一定满足 KKT 条件,并且满足强对偶性质;

a

a

a

a

a

a

2. ALM 算法

注1:ALM(Augmented Lagrange Multiplier)算法,即增广型拉格朗日乘子法

注2:ALM 算法常用于线性规划问题( f(x) 不能有高阶项/根号项,也不能有两个变量之间的交乘积项);

注3:ALM 算法的推导过程 -- http://faculty.bicmr.pku.edu.cn/~wenzw/optbook/lect/16-lect-alm.pdf


2.1 ALM 算法介绍

只考虑含有等式约束的问题:

minf(x)

Ax-b=0

x\geq 0

a

a

注:含有不等式约束问题的 ALM 算法,见文章 -- http://faculty.bicmr.pku.edu.cn/~wenzw/optbook/lect/16-lect-alm.pdf

构造增广拉格朗日函数:

L(x,\lambda )=f(x)+\lambda (Ax-b)+\frac{\alpha }{2}\left \| Ax-b \right \|_2^2

  • \lambda:拉格朗日乘子(矩阵)
  • \alpha:惩罚项权重(标量)

a

a

注意:增广拉格朗日乘子法就是在拉格朗日乘子法获得的函数后面,加上针对所有等式约束的惩罚项;

a

a

优点:

  • 原拉格朗日函数的收敛性不太好,抖动很大不稳定;
  • 加上了惩罚项可以增加原拉格朗日函数的凸性,使得我们可以通过更简单的方法,如梯度下降法去求解这个函数的最优解;
  • 加强了等式约束的限制作用(针对等式约束增加了惩罚项),使得在迭代的过程中迭代点一直是在约束附近的区域进行梯度下降,不会跑太远,从而加快了求解速度;

ALM 算法的求解过程:梯度下降

  • 先对变量 x 进行梯度下降:x_{k+1}=x_{k}-\eta \bigtriangledown _{x}L
  • 再对拉格朗日乘子 \lambda 进行梯度下降:\lambda _{k+1}=\lambda _{k}+\alpha (Ax_{k+1}-b)
  • 上述两步为一次迭代,不断这样迭代下去,直至收敛或者到达指定迭代次数;

2.2 ALM 算法代码示例

举例:只含等式约束

minf(x)=C^Tx

Ax-b=0

x\geq 0

a

使用 Scipy 中的函数 linprog 求解线性规划问题,将求解得到的结果作为参考答案:

函数 linprog:scipy.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.000e+00  6.000e+00  2.000e+00  0.000e+00  0.000e+00]
"""from scipy.optimize import linprog
import torchc = torch.tensor([-3, -5, 0, 0, 0], dtype=torch.float32)
A = torch.tensor([[1, 0, 1, 0, 0], [0, 2, 0, 1, 0], [3, 2, 0, 0, 1]], dtype=torch.float32)
b = torch.tensor([4, 12, 18], dtype=torch.float32)res = linprog(c, A_eq=A, b_eq=b, 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(min=0)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{i+1}th iter, L:{lagrangian_function(x, lambda_):.2f}, f: {f(x):.2f}')print(f'x: {x}')print(f'lambda: {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], dtype=torch.float32)A = torch.tensor([[1, 0, 1, 0, 0], [0, 2, 0, 1, 0], [3, 2, 0, 0, 1]], dtype=torch.float32)b = torch.tensor([4, 12, 18], dtype=torch.float32)lambda_ = torch.tensor([0, 0, 0], dtype=torch.float32)x = torch.tensor([2, 0, 2, 0, 0], dtype=torch.float32, requires_grad=True)solve(x, lambda_)

a

a

a

a

a

a

3. ADMM 算法

注1:ADMM(Alternating Direction Method of Multipliers)算法,即交替方向乘子法

注2:关于详细的数学推导 -- 一文详解从拉格朗日乘子法、KKT条件、对偶上升法到罚函数与增广Lagrangian乘子法再到ADMM算法(交替方向乘子法)_拉格朗日乘子 二次惩罚系数-CSDN博客


3.1 ADMM 算法介绍

 只考虑含有等式约束的问题:

minf(x,y)

Ax+By-b=0

x,y\geq 0

a

a

和ALM算法的不同点在于:将所有的变量分成几部分(假设为两部分,一部分是向量 x,一部分是向量 y),然后分别进行梯度下降,而不是一次性将所有的变量全部都进行梯度下降;

构造增广拉格朗日函数:

L(x,y,\lambda )=f(x,y)+\lambda (Ax+By-b)+\frac{\alpha }{2}\left \| Ax+By-b \right \|_2^2

  • \lambda:拉格朗日乘子(矩阵)
  • \alpha:惩罚项权重(标量)

ADMM 算法的求解过程:

  • 先对变量 x 进行梯度下降:x_{k+1}=x_{k}-\eta \bigtriangledown _{x}L(x_k,y_k,\lambda _k)
  • 再对变量 y 进行梯度下降:y_{k+1}=y_{k}-\eta \bigtriangledown _{y}L(x_{k+1},y_k,\lambda _k)
  • 再对变量 \lambda 进行梯度下降:\lambda_{k+1}=\lambda_{k}+\alpha (Ax_{k+1}+By_{k+1}-b)
  • 上述三步为一次迭代,不断这样迭代下去,直至收敛或者到达指定迭代次数;

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{i+1}th iter, L:{lagrangian_function(x, lambda_):.2f}, f: {f(x):.2f}')print(f'x: {x}')print(f'lambda: {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], dtype=torch.float32)A = torch.tensor([[1, 0, 1, 0, 0], [0, 2, 0, 1, 0], [3, 2, 0, 0, 1]], dtype=torch.float32)b = torch.tensor([4, 12, 18], dtype=torch.float32)lambda_ = torch.tensor([0, 0, 0], dtype=torch.float32)x = torch.tensor([2, 0, 2, 0, 0], dtype=torch.float32, requires_grad=True)solve(x, lambda_)

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/860903.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

图书馆书籍管理系统

项目名称与项目简介 图书馆书籍管理系统 本项目是一个计算机管理系统,也就是将传统手工的管理方式转变为智能化、标准化、规范化的管理管理模式,对图书馆中所有的图书、文献资料、音像资料、报刊、期刊等各种类型的资料实现采编、收集图书信息、检索、归…

探讨数字化背景下VSM(价值流程图)的挑战和机遇

在信息化、数字化飞速发展的今天,各行各业都面临着前所未有的挑战与机遇。作为源自丰田生产模式的VSM(价值流程图),这一曾经引领制造业革命的工具,在数字化背景下又将如何乘风破浪,应对新的市场格局和技术变…

Java中的方法重载

方法重载(Method Overloading)是Java中面向对象编程的重要特性之一。它允许一个类中有多个方法名称相同但参数列表不同的方法,从而实现不同的功能。本文将深入探讨方法重载的概念、规则及其应用场景,并通过示例代码展示其实际用法…

屏幕翻译下载哪个软件好?好用的屏幕翻译推荐

想象一下,当我们在阅读外文文档或是观看外语电影时,如果能有一款翻译工具同步提供译文,那将是多么令人愉悦的体验! 如果这种翻译服务能够在不影响其他应用的情况下进行,那就是double快乐了。 其实,现在要…

JAVA之Guava

Guava 是由 Google 开发的一个开源 Java 库,包含了大量被认为是 Java 编程中常用的实用工具类和方法。Guava 提供了丰富的集合处理、缓存、字符串处理、并发工具和 I/O 操作等功能。 1. 集合工具 Guava 提供了强大的集合框架扩展,简化了集合的创建、操…

【LeetCode】每日一题:二叉树的层次遍历

给你二叉树的根节点 root ,返回其节点值的 层序遍历 。 (即逐层地,从左到右访问所有节点)。 解题思路 水题 AC代码 # Definition for a binary tree node. # class TreeNode: # def __init__(self, val0, leftNone, rightN…

【ajax07基础】回调函数地狱

一:什么是回调函数地狱 在一个回调函数中嵌套另一个回调函数(甚至一直嵌套下去),形成回调函数地狱 回调函数地狱存在问题: 可读性差异常捕获严重耦合性严重 // 1. 获取默认第一个省份的名字axios({url: http://hmaj…

Postman使用简要步骤

目录 下载与安装 Postman 启动 Postman 创建请求 设置请求信息 发送请求 编写测试脚本 使用环境变量和全局变量 使用 Pre-request Script 生成 API 文档 模拟和测试 API 持续集成与自动化测试 下载与安装 Postman 访问 Postman 官网(https://www.postman…

在 C 语言中使用 UT_hash_handle 简化实现哈希表

在 C 语言中使用 UT_hash_handle 简化实现哈希表的步骤如下: 1、包含头文件:首先需要包含 uthash 库的头文件 uthash.h。 #include "uthash.h"2、定义包含哈希表支持的结构体:在定义结构体时,需要在结构体中添加 UT_h…

浏览器自动填充登录用户名和密码,如何清除

文章目录 刷新网页的时候浏览器会自动填充用户名和密码刷新之后效果图解决方案完整的login.vue代码核心代码原理(添加 readonly 和监听 focus 事件) 刷新网页的时候浏览器会自动填充用户名和密码 刷新之后效果图 解决方案 完整的login.vue代码 <template><div class…

【线上】如何解决积压消费?

我是小米,一个喜欢分享技术的29岁程序员。如果你喜欢我的文章,欢迎关注我的微信公众号“软件求生”,获取更多技术干货!​​​​​​​ Hello, 各位亲爱的读者朋友们!我是你们的小米,一个积极活泼的技术分享达人,今天我们要聊聊一个大家在分布式系统中经常遇到的棘手问题…

【Linux】使用ntp同步时间

ntp介绍 NTP&#xff08;Network Time Protocol&#xff0c;网络时间协议&#xff09;是一种用于同步计算机时间的协议&#xff0c;工作在UDP的123端口上。它是一种客户端-服务器协议&#xff0c;用于同步计算机的时钟。通过连接到网络上的时间服务器&#xff0c;计算机可以获…

PFA 反应罐盖特氟龙 润滑绝缘行业加工 匠心工艺

PFA反应罐别名也叫反应瓶&#xff0c;储样罐&#xff0c;清洗罐等。可作为样品前处理实验中消解样品和中低压溶样的反应容器&#xff0c;广泛应用于半导体分析、新材料、新能源、同位素分析等。 PFA反应罐规格参考&#xff1a;250ml、300ml、350ml、500ml、1L等。 产品特点&am…

项目管理:如何解决项目延期的那些问题?

在项目管理中&#xff0c;项目延期是一种普遍现象&#xff0c;也管理者最为头疼的一个问题。为了有效地解决项目延期问题&#xff0c;我们需要从多个方面入手&#xff1a; 1、快速识别原因&#xff1a;当项目出现延期迹象时&#xff0c;首要任务是迅速识别并定位导致延期的根…

【语言模型】探索AI模型、AI大模型、大模型、大语言模型与大数据模型的关系与协同

一、引言 随着人工智能&#xff08;AI&#xff09;技术的飞速发展&#xff0c;各种AI模型如雨后春笋般涌现&#xff0c;其中AI模型、AI大模型、大模型、大语言模型以及大数据模型等概念在学术界和工业界引起了广泛关注。这些模型不仅各自具有独特的特点和应用场景&#xff0c;…

【深度学习】基于因果表示学习的CITRIS模型原理和实验

1.引言 1.1.本文的主要内容 理解动态系统中的潜在因果因素&#xff0c;对于智能代理在复杂环境中进行有效推理至关重要。本文将深入介绍CITRIS&#xff0c;这是一种基于变分自编码器&#xff08;VAE&#xff09;的框架&#xff0c;它能够从时间序列图像中提取并学习因果表示&…

QT QSlider控件-主介绍 触发函数常用函数

QSlider控件是Qt库中用于提供一个可拖动滑块以选择数值或范围的界面元素。它广泛应用于需要用户进行数值调节的场景&#xff0c;如音量控制、亮度调整等。 一、QAbstractSlider的6个信号量触发函数&#xff1a; 1、void actionTriggered (int action): 当滑块上的某个可定义动…

量化系统----开源简化版qmt实盘交易系统,提供源代码

量化系统----开源简化版qmt实盘交易系统&#xff0c;提供源代码 https://mp.weixin.qq.com/s/qeqH8XtUeoDjIJIXMe5D-w 最近有读者反应开源的qmt_trader内容太多了不知道怎么样使用&#xff0c;我独立了一个简单板块的easy_qmt_tarder方面大家的使用 qmt_tarder开源下载 量化系…

企业网络安全策略转型

企业网络安全策略的转型是伴随着数字化转型的深入推进而产生的&#xff0c;它涉及到从传统的安全防护向智能化、主动防御的转变。在这一过程中&#xff0c;企业不仅要面对日益复杂的网络威胁&#xff0c;还需适应新的技术环境&#xff0c;如云计算、大数据等新兴技术带来的挑战…

【网络通信】计算机网络安全技术总结

一、概述 在数字时代的浪潮下&#xff0c;计算机网络安全技术已成为保护数据完整性和安全性的基石。这项技术不仅是计算机科学的重要组成部分&#xff0c;也是应对各种网络威胁和挑战的关键手段。 二、核心技术和应用 2.1 加密技术 作为网络安全技术的核心&#xff0c;加密技…