cvxpy简易入门教程

2022-04-10

本文案例翻译整理自Welcome to CVXPY 1.2 — CVXPY 1.2 documentation

百度翻了翻，发现几乎没有cvxpy的中文入门教程，感觉这个库计算凸优化问题时挺有用的，所以参照文档试着写写入门教程，自己也顺便学习一下，有所纰漏见谅。

安装

我的mac是直接python3.9 -m pip install cvxpy就行了（我装了好几个python所以需要声明是3.9的pip，一般正常pip3 install *** 就行），安装很顺利，看网上貌似win比较麻烦，我没试过。

如果遇到报错，可能是因为：

缺少CMake导致无法构建qdldl的wheel，从而导致安装cvxpy失败。在报错信息中可以看到以下关键信息：

RuntimeError: CMake must be installed to build qdldl：缺少CMake导致无法构建qdldl。

FileNotFoundError: [Errno 2] No such file or directory: 'cmake'：找不到’CMake’文件或目录。

解决方法是安装CMake。你可以通过以下步骤解决这个问题：
安装CMake：在终端中运行以下命令安装CMake：
1
brew install cmake
如果你没有Homebrew，请先安装Homebrew，然后再运行上述命令。
安装cvxpy：再次尝试运行python3.9 -m pip install cvxpy安装cvxpy。
安装了CMake后，应该能够成功构建qdldl并安装cvxpy了。

CVXPY 1.2

CVXPY是一种用于凸优化问题的开源Python库。它允许您以自然的方式表达数学问题，而不是以编程求解要求的程序化标准形式表达问题，说白了就是，你可以直接说人话，不用说计算机鸟语。

例如，以下代码解决了简单界约束的最小二乘问题：

# Import packages.
import cvxpy as cp
import numpy as np

# Generate a random non-trivial【非平凡】 linear program.
m = 15
n = 10
np.random.seed(1)
s0 = np.random.randn(m)
lamb0 = np.maximum(-s0, 0)
s0 = np.maximum(s0, 0)
x0 = np.random.randn(n)
A = np.random.randn(m, n)
b = A @ x0 + s0
c = -A.T @ lamb0

# Define and solve the CVXPY problem.
x = cp.Variable(n)
prob = cp.Problem(cp.Minimize(c.T@x),
                 [A @ x <= b])
prob.solve()

# Print result.
print("\nThe optimal value is", prob.value)
print("A solution x is")
print(x.value)
print("A dual solution is")
print(prob.constraints[0].dual_value)

结果：

The optimal value is -15.220912605552863
A solution x is
[-1.10133381 -0.16360111 -0.89734939  0.03216603  0.6069123  -1.12687348
  1.12967856  0.88176638  0.49075229  0.8984822 ]
A dual solution is
[6.98805172e-10 6.11756416e-01 5.28171747e-01 1.07296862e+00
 3.93759300e-09 2.30153870e+00 4.25704434e-10 7.61206896e-01
 8.36906030e-09 2.49370377e-01 1.30187120e-09 2.06014070e+00
 3.22417207e-01 3.84054343e-01 1.59493839e-09]

先不用管这堆程序是怎么写出来的，后面会慢慢介绍。

最小二乘¶

在最小二乘或线性回归问题中，我们有 $A \in \mathcal{R}^{m \times n}$ 编辑和 $b \in \mathcal{R}^m$ 编辑，寻找一个向量 $x \in \mathcal{R}^{n}$ 编辑使得 $Ax$ 编辑接近 $b$ 编辑。接近度定义为平方差的和：

$\sum_{i=1}^m (a_i^Tx - b_i)^2,$ 编辑

也称为 $\ell_2$ 编辑规范平方， $\|Ax - b\|_2^2$ 编辑。

在以下代码中，我们用CVXPY解决了最小二乘问题

# Import packages.
import cvxpy as cp
import numpy as np

# Generate data.
m = 20
n = 15
np.random.seed(1)
A = np.random.randn(m, n)
b = np.random.randn(m)

# Define and solve the CVXPY problem.
x = cp.Variable(n)
cost = cp.sum_squares(A @ x - b)
prob = cp.Problem(cp.Minimize(cost))
prob.solve()

# Print result.
print("\nThe optimal value is", prob.value)
print("The optimal x is")
print(x.value)
print("The norm of the residual is ", cp.norm(A @ x - b, p=2).value)

结果如下：

The optimal value is 7.005909828287485
The optimal x is
[ 0.17492418 -0.38102551  0.34732251  0.0173098  -0.0845784  -0.08134019
  0.293119    0.27019762  0.17493179 -0.23953449  0.64097935 -0.41633637
  0.12799688  0.1063942  -0.32158411]
The norm of the residual is  2.6468679280023557

☺️：

这里A为20*15的矩阵，b为20维列向量。

x = cp.Variable(n) 定义了一个叫x的变量，它是一个15维列向量，具体数值这一步不确定。cp.sum_squares函数就是计算平方和的函数，prob=cp.Problem() 定义了一个“问题”，“问题”函数里填写凸优化的目标，目前的目标就是那个“平方和”cost最小，使用cp.Minimize函数表示。prob.solve() 求解，运行完这一步才能确定x的具体数值。

这里也能看出cvxpy编程书写步骤与自然语言接近，完全不是一般编程那样定义变量需要有确定的值，然后再循环调整之类的繁琐语言顺序。

prob.value储存的是minimize(cost)的值，就是优化后目标的值。查看变量x使用x.value

关于cp.norm计算向量范数：

线性规划

常见的标准形式如下：

$\begin{array}{ll} \mbox{minimize} & c^Tx \\ \mbox{subject to} & Ax \leq b. \end{array}$ 编辑

考虑到这个一般中国人初中就学过了，不解释了。直接给例子：

# Import packages.
import cvxpy as cp
import numpy as np

# Generate a random non-trivial linear program.
m = 15
n = 10
np.random.seed(1)
s0 = np.random.randn(m)
lamb0 = np.maximum(-s0, 0)
s0 = np.maximum(s0, 0)
x0 = np.random.randn(n)
A = np.random.randn(m, n)
b = A @ x0 + s0
c = -A.T @ lamb0

# Define and solve the CVXPY problem.
x = cp.Variable(n)
prob = cp.Problem(cp.Minimize(c.T@x),
                 [A @ x <= b])
prob.solve()

# Print result.
print("\nThe optimal value is", prob.value)
print("A solution x is")
print(x.value)
print("A dual solution is")
print(prob.constraints[0].dual_value)

结果：

The optimal value is -15.220912605552863
A solution x is
[-1.10133381 -0.16360111 -0.89734939  0.03216603  0.6069123  -1.12687348
  1.12967856  0.88176638  0.49075229  0.8984822 ]
A dual solution is
[6.98805172e-10 6.11756416e-01 5.28171747e-01 1.07296862e+00
 3.93759300e-09 2.30153870e+00 4.25704434e-10 7.61206896e-01
 8.36906030e-09 2.49370377e-01 1.30187120e-09 2.06014070e+00
 3.22417207e-01 3.84054343e-01 1.59493839e-09]

hhgw：

s0，lamb0，x0是构造线性规划方程组，虽然规划得有点绕，但是不重要，
prob = cp.Problem(cp.Minimize(c.T@x),[A @ x <= b])
看懂这一步就行了。

后面dual_value指线性规划的对偶规划的解。

二次规划

未完待续。