目标函数为二次式，约束条件为线性式的最优化问题称为二次规划。其一般形式为
{
minimize
1
2
x
⊤
H
x
+
c
⊤
x
s.t.
A
e
q
x
−
b
e
q
=
o
A
i
q
x
−
b
i
q
≥
o
.
\begin{cases} \text{minimize}\quad \frac{1}{2}\boldsymbol{x}^\top\boldsymbol{Hx}+\boldsymbol{c}^\top\boldsymbol{x}\\ \text{s.t.\ \ }\quad\quad\quad\boldsymbol{A}_{eq}\boldsymbol{x}-\boldsymbol{b}_{eq}=\boldsymbol{o}\\ \quad\quad\quad\quad\quad\boldsymbol{A}_{iq}\boldsymbol{x}-\boldsymbol{b}_{iq}\geq\boldsymbol{o} \end{cases}.
⎩
⎨
⎧minimize21x⊤Hx+c⊤xs.t. Aeqx−beq=oAiqx−biq≥o.
其中，
H
∈
R
n
×
n
\boldsymbol{H}\in\text{R}^{n\times n}
H∈Rn×n对称，
c
∈
R
n
\boldsymbol{c}\in\text{R}^n
c∈Rn，
A
e
q
∈
R
l
×
n
\boldsymbol{A}_{eq}\in\text{R}^{l\times n}
Aeq∈Rl×n，
b
e
q
∈
R
l
\boldsymbol{b}_{eq}\in\text{R}^l
beq∈Rl，
A
i
q
∈
R
m
×
n
\boldsymbol{A}_{iq}\in\text{R}^{m\times n}
Aiq∈Rm×n，
b
i
q
∈
R
m
\boldsymbol{b}_{iq}\in\text{R}^m
biq∈Rm。
仅含等式约束的二次规划形如
{
minimize
1
2
x
⊤
H
x
+
c
⊤
x
s.t.
A
x
−
b
=
o
.
(
1
)
\begin{cases} \text{minimize}\quad \frac{1}{2}\boldsymbol{x}^\top\boldsymbol{Hx}+\boldsymbol{c}^\top\boldsymbol{x}\\ \text{s.t.\ \ }\quad\quad\quad\boldsymbol{Ax}-\boldsymbol{b}=\boldsymbol{o} \end{cases}.\quad\quad(1)
{minimize21x⊤Hx+c⊤xs.t. Ax−b=o.(1)
假定
H
\boldsymbol{H}
H对称正定，
A
∈
R
l
×
n
\boldsymbol{A}\in\text{R}^{l\times n}
A∈Rl×n，rank
A
=
l
\boldsymbol{A}=l
A=l。正定二次式
1
2
x
⊤
H
x
+
c
⊤
x
\frac{1}{2}\boldsymbol{x}^\top\boldsymbol{Hx}+\boldsymbol{c}^\top\boldsymbol{x}
21x⊤Hx+c⊤x在凸集
Ω
=
{
x
∣
A
x
−
b
=
o
}
\Omega=\{\boldsymbol{x}|\boldsymbol{Ax}-\boldsymbol{b}=\boldsymbol{o}\}
Ω={x∣Ax−b=o}上有唯一满足必要条件的KKT点
(
x
0
λ
0
)
\begin{pmatrix}\boldsymbol{x}_0\\\boldsymbol{\lambda}_0\end{pmatrix}
(x0λ0)。为算得该KKT点，写出问题的拉格朗日函数
L
(
x
,
λ
)
=
1
2
x
⊤
H
x
+
c
⊤
x
−
λ
⊤
(
A
x
−
b
)
.
L(\boldsymbol{x},\boldsymbol{\lambda})=\frac{1}{2}\boldsymbol{x}^\top\boldsymbol{Hx}+\boldsymbol{c}^\top\boldsymbol{x}-\boldsymbol{\lambda}^\top(\boldsymbol{Ax}-\boldsymbol{b}).
L(x,λ)=21x⊤Hx+c⊤x−λ⊤(Ax−b).
关于
x
\boldsymbol{x}
x，
λ
\boldsymbol{\lambda}
λ的梯度为
∇
x
L
(
x
,
λ
)
=
H
x
+
c
−
A
⊤
λ
∇
λ
L
(
x
,
λ
)
=
−
A
x
+
b
\begin{array}{l} \nabla_{\boldsymbol{x}}L(\boldsymbol{x},\boldsymbol{\lambda})=\boldsymbol{Hx}+\boldsymbol{c}-\boldsymbol{A}^\top\boldsymbol{\lambda}\\ \nabla_{\boldsymbol{\lambda}}L(\boldsymbol{x},\boldsymbol{\lambda})=-\boldsymbol{Ax}+\boldsymbol{b} \end{array}
∇xL(x,λ)=Hx+c−A⊤λ∇λL(x,λ)=−Ax+b
令
∇
x
L
(
x
,
λ
)
=
o
\nabla_{\boldsymbol{x}}L(\boldsymbol{x},\boldsymbol{\lambda})=\boldsymbol{o}
∇xL(x,λ)=o且
∇
λ
L
(
x
,
λ
)
=
o
\nabla_{\boldsymbol{\lambda}}L(\boldsymbol{x},\boldsymbol{\lambda})=\boldsymbol{o}
∇λL(x,λ)=o，得线性方程组
{
H
x
+
c
−
A
⊤
λ
=
o
−
A
x
+
b
=
o
,
\begin{cases} \boldsymbol{Hx}+\boldsymbol{c}-\boldsymbol{A}^\top\boldsymbol{\lambda}=\boldsymbol{o}\\ -\boldsymbol{Ax}+\boldsymbol{b}=\boldsymbol{o} \end{cases},
{Hx+c−A⊤λ=o−Ax+b=o,
等价地表示为
(
H
−
A
⊤
−
A
O
)
(
x
λ
)
=
(
−
c
−
b
)
.
\begin{pmatrix}\boldsymbol{H}&-\boldsymbol{A}^\top\\-\boldsymbol{A}&\boldsymbol{O}\end{pmatrix}\begin{pmatrix}\boldsymbol{x}\\\boldsymbol{\lambda}\end{pmatrix} =\begin{pmatrix}-\boldsymbol{c}\\-\boldsymbol{b}\end{pmatrix}.
(H−A−A⊤O)(xλ)=(−c−b).
系数矩阵
(
H
−
A
⊤
−
A
O
)
\begin{pmatrix}\boldsymbol{H}&-\boldsymbol{A}^\top\\-\boldsymbol{A}&\boldsymbol{O}\end{pmatrix}
(H−A−A⊤O)称为拉格朗日矩阵。由
H
\boldsymbol{H}
H对称正定且rank
A
=
l
\boldsymbol{A}=l
A=l的假设，拉格朗日矩阵可逆，设
(
H
−
A
⊤
−
A
O
)
−
1
=
(
Q
−
R
⊤
−
R
S
)
,
\begin{pmatrix}\boldsymbol{H}&-\boldsymbol{A}^\top\\-\boldsymbol{A}&\boldsymbol{O}\end{pmatrix}^{-1}=\begin{pmatrix}\boldsymbol{Q}&-\boldsymbol{R}^\top\\-\boldsymbol{R}&\boldsymbol{S}\end{pmatrix},
(H−A−A⊤O)−1=(Q−R−R⊤S),
根据
(
H
−
A
⊤
−
A
O
)
(
Q
−
R
⊤
−
R
S
)
=
I
\begin{pmatrix}\boldsymbol{H}&-\boldsymbol{A}^\top\\-\boldsymbol{A}&\boldsymbol{O}\end{pmatrix}\begin{pmatrix}\boldsymbol{Q}&-\boldsymbol{R}^\top\\-\boldsymbol{R}&\boldsymbol{S}\end{pmatrix}=\boldsymbol{I}
(H−A−A⊤O)(Q−R−R⊤S)=I
算得
{
H
Q
+
A
⊤
R
=
I
−
H
R
⊤
−
A
⊤
S
=
O
−
A
Q
=
O
A
R
⊤
=
I
\begin{cases} \boldsymbol{HQ}+\boldsymbol{A}^\top\boldsymbol{R}=\boldsymbol{I}\\ -\boldsymbol{H}\boldsymbol{R}^\top-\boldsymbol{A}^\top\boldsymbol{S}=\boldsymbol{O}\\ -\boldsymbol{AQ}=\boldsymbol{O}\\ \boldsymbol{AR}^\top=\boldsymbol{I} \end{cases}
⎩
⎨
⎧HQ+A⊤R=I−HR⊤−A⊤S=O−AQ=OAR⊤=I
由于
A
\boldsymbol{A}
A行满秩，故
A
H
−
1
A
⊤
\boldsymbol{AH}^{-1}\boldsymbol{A}^\top
AH−1A⊤可逆。
(
A
H
−
1
A
⊤
)
−
1
A
H
−
1
(\boldsymbol{AH}^{-1}\boldsymbol{A}^\top)^{-1}\boldsymbol{AH}^{-1}
(AH−1A⊤)−1AH−1是
A
⊤
\boldsymbol{A}^\top
A⊤的伪逆。解上列连立式得
{
S
=
−
(
A
H
−
1
A
⊤
)
−
1
R
=
−
S
A
H
−
1
Q
=
H
−
1
−
H
−
1
A
⊤
R
\begin{cases}\boldsymbol{S}=-(\boldsymbol{AH}^{-1}\boldsymbol{A}^\top)^{-1}\\\boldsymbol{R}=-\boldsymbol{SAH}^{-1}\\\boldsymbol{Q}=\boldsymbol{H}^{-1}-\boldsymbol{H}^{-1}\boldsymbol{A}^\top\boldsymbol{R}\end{cases}
⎩
⎨
⎧S=−(AH−1A⊤)−1R=−SAH−1Q=H−1−H−1A⊤R
于是，二次规划(1)的KKT点
(
x
0
λ
0
)
=
(
Q
−
R
⊤
−
R
S
)
(
−
c
−
b
)
=
(
−
Q
c
+
R
⊤
b
R
c
−
S
b
)
.
\begin{pmatrix}\boldsymbol{x}_0\\\boldsymbol{\lambda}_0\end{pmatrix}=\begin{pmatrix}\boldsymbol{Q}&-\boldsymbol{R}^\top\\-\boldsymbol{R}&\boldsymbol{S}\end{pmatrix}\begin{pmatrix}-\boldsymbol{c}\\-\boldsymbol{b}\end{pmatrix}=\begin{pmatrix}-\boldsymbol{Qc}+\boldsymbol{R}^\top\boldsymbol{b}\\\boldsymbol{Rc}-\boldsymbol{Sb} \end{pmatrix}.
(x0λ0)=(Q−R−R⊤S)(−c−b)=(−Qc+R⊤bRc−Sb).
下列代码实现求解等式约束二次规划(1)的拉格朗日算法。
import numpy as np #导入numpy
def qlag(H, A, b, c):
H1 = np.linalg.inv(H) #H的逆阵
S = -np.linalg.inv(np.matmul(np.matmul(A, H1), A.T))
R = -np.matmul(np.matmul(S, A), H1)
Q = H1 - np.matmul(np.matmul(H1, A.T), R)
x0 = -np.matmul(Q, c) + np.matmul(R.T, b) #最优解
lamd0 = np.matmul(R, c) - np.matmul(S, b) #拉格朗日乘子
return x0, lamd0
程序的第2~9行定义的函数qlag实现拉格朗日算法。qlag的4个参数H，A，b和c分别表示二次规划(1)中的正定矩阵
H
\boldsymbol{H}
H，行满秩阵
A
\boldsymbol{A}
A，向量
b
\boldsymbol{b}
b和
c
\boldsymbol{c}
c。
函数体内的第3行调用numpy.linalg的inv函数计算
H
\boldsymbol{H}
H的逆阵
H
−
1
\boldsymbol{H}^{-1}
H−1，赋予H1。第4~6行分别计算
S
=
−
(
A
H
−
1
A
⊤
)
−
1
R
=
−
S
A
H
−
1
Q
=
H
−
1
−
H
−
1
A
⊤
R
\begin{array}{l} \boldsymbol{S}=-(\boldsymbol{AH}^{-1}\boldsymbol{A}^\top)^{-1}\\ \boldsymbol{R}=-\boldsymbol{SAH}^{-1}\\ \boldsymbol{Q}=\boldsymbol{H}^{-1}-\boldsymbol{H}^{-1}\boldsymbol{A}^\top\boldsymbol{R} \end{array}
S=−(AH−1A⊤)−1R=−SAH−1Q=H−1−H−1A⊤R
并赋予S，R和Q。第7、8行分别计算最优解和对应的拉格朗日乘子
x
0
=
−
Q
c
+
R
⊤
b
λ
0
=
R
c
−
S
b
\begin{array}{l} \boldsymbol{x}_0=-\boldsymbol{Qc}+\boldsymbol{R}^\top\boldsymbol{b}\\ \boldsymbol{\lambda}_0=\boldsymbol{Rc}-\boldsymbol{Sb} \end{array}
x0=−Qc+R⊤bλ0=Rc−Sb
并赋予x0和lamd0。
例1用qlag函数求解下列二次规划
{
minimize
x
1
2
+
2
x
2
2
+
x
3
2
−
2
x
1
x
2
+
x
3
s.t.
x
1
+
x
2
+
x
3
=
4
2
x
1
−
x
2
+
x
3
=
2
.
\begin{cases} \text{minimize}\quad x_1^2+2x_2^2+x_3^2-2x_1x_2+x_3\\ \text{s.t.\ \ }\quad\quad\quad x_1+x_2+x_3=4\\ \quad\quad\quad\quad\quad 2x_1-x_2+x_3=2 \end{cases}.
⎩
⎨
⎧minimizex12+2x22+x32−2x1x2+x3s.t. x1+x2+x3=42x1−x2+x3=2.
解:本问题中，
H
=
(
2
−
2
0
−
2
4
0
0
0
2
)
,
A
=
(
1
1
1
2
−
1
1
)
,
b
=
(
4
2
)
,
c
=
(
0
0
1
)
\boldsymbol{H}=\begin{pmatrix}2&-2&0\\-2&4&0\\0&0&2\end{pmatrix},\boldsymbol{A}=\begin{pmatrix}1&1&1\\2&-1&1\end{pmatrix},\boldsymbol{b}=\begin{pmatrix}4\\2\end{pmatrix},\boldsymbol{c}=\begin{pmatrix}0\\0\\1\end{pmatrix}
H=
2−20−240002
,A=(121−111),b=(42),c=
001
下列代码利用这些数据进行计算。
import numpy as np #导入numpy
from fractions import Fraction as F #设置输出格式
np.set_printoptions(formatter={'all':lambda x:
str(F(x).limit_denominator())})
H = np.array([[2, -2, 0], #矩阵H
[-2, 4, 0],
[0, 0, 2]])
A = np.array([[1, 1, 1], #矩阵A
[2, -1, 1]])
b = np.array([4, 2]) #向量b
c = np.array([0, 0, 1]) #向量c
print(qlag(H, A, b, c)) #计算最优解
程序的第2~4行设置数组输出格式为有理数。5~11设置表示本二次规划问题的矩阵H、A和向量b、c。第12行调用函数qlag，计算本二次规划最优解。运行程序，输出
(array([21/11, 43/22, 3/22]), array([29/11, -15/11]))
意味着最优解
x
0
=
(
21
11
43
22
3
22
)
\boldsymbol{x}_0=\begin{pmatrix}\frac{21}{11}\\\frac{43}{22}\\\frac{3}{22}\end{pmatrix}
x0=
11212243223
，对应的拉格朗日乘子
λ
0
=
(
29
11
−
15
11
)
\boldsymbol{\lambda}_0=\begin{pmatrix}\frac{29}{11}\\-\frac{15}{11}\end{pmatrix}
λ0=(1129−1115)。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！
扫一扫
40万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
