最优化方法Python计算：二次规划的拉格朗日算法
目标函数为二次式，约束条件为线性式的最优化问题称为二次规划。其一般形式为
{minimize12x⊤Hx+c⊤xs.t. Aeqx−beq=oAiqx−biq≥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∈Rn×n\boldsymbol{H}\in\text{R}^{n\times n}H∈Rn×n对称，c∈Rn\boldsymbol{c}\in\text{R}^nc∈Rn，Aeq∈Rl×n\boldsymbol{A}_{eq}\in\text{R}^{l\times n}Aeq∈Rl×n，beq∈Rl\boldsymbol{b}_{eq}\in\text{R}^lbeq∈Rl，Aiq∈Rm×n\boldsymbol{A}_{iq}\in\text{R}^{m\times n}Aiq∈Rm×n，biq∈Rm\boldsymbol{b}_{iq}\in\text{R}^mbiq∈Rm。
仅含等式约束的二次规划形如
{minimize12x⊤Hx+c⊤xs.t. Ax−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∈Rl×n\boldsymbol{A}\in\text{R}^{l\times n}A∈Rl×n，rankA=l\boldsymbol{A}=lA=l。正定二次式12x⊤Hx+c⊤x\frac{1}{2}\boldsymbol{x}^\top\boldsymbol{Hx}+\boldsymbol{c}^\top\boldsymbol{x}21x⊤Hx+c⊤x在凸集Ω={x∣Ax−b=o}\Omega=\{\boldsymbol{x}|\boldsymbol{Ax}-\boldsymbol{b}=\boldsymbol{o}\}Ω={x∣Ax−b=o}上有唯一满足必要条件的KKT点(x0λ0)\begin{pmatrix}\boldsymbol{x}_0\\\boldsymbol{\lambda}_0\end{pmatrix}(x0λ0)。为算得该KKT点，写出问题的拉格朗日函数
L(x,λ)=12x⊤Hx+c⊤x−λ⊤(Ax−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}λ的梯度为
∇xL(x,λ)=Hx+c−A⊤λ∇λL(x,λ)=−Ax+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
令∇xL(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，得线性方程组
{Hx+c−A⊤λ=o−Ax+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⊤−AO)(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⊤−AO)\begin{pmatrix}\boldsymbol{H}&-\boldsymbol{A}^\top\\-\boldsymbol{A}&\boldsymbol{O}\end{pmatrix}(H−A−A⊤O)称为拉格朗日矩阵。由H\boldsymbol{H}H对称正定且rankA=l\boldsymbol{A}=lA=l的假设，拉格朗日矩阵可逆，设
(H−A⊤−AO)−1=(Q−R⊤−RS),\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⊤−AO)(Q−R⊤−RS)=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
算得
{HQ+A⊤R=I−HR⊤−A⊤S=O−AQ=OAR⊤=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行满秩，故AH−1A⊤\boldsymbol{AH}^{-1}\boldsymbol{A}^\topAH−1A⊤可逆。(AH−1A⊤)−1AH−1(\boldsymbol{AH}^{-1}\boldsymbol{A}^\top)^{-1}\boldsymbol{AH}^{-1}(AH−1A⊤)−1AH−1是A⊤\boldsymbol{A}^\topA⊤的伪逆。解上列连立式得
{S=−(AH−1A⊤)−1R=−SAH−1Q=H−1−H−1A⊤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点
(x0λ0)=(Q−R⊤−RS)(−c−b)=(−Qc+R⊤bRc−Sb).\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=−(AH−1A⊤)−1R=−SAH−1Q=H−1−H−1A⊤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行分别计算最优解和对应的拉格朗日乘子
x0=−Qc+R⊤bλ0=Rc−Sb\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函数求解下列二次规划
{minimizex12+2x22+x32−2x1x2+x3s.t. x1+x2+x3=42x1−x2+x3=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−20−240002),A=(1112−11),b=(42),c=(001)\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]))
意味着最优解x0=(21114322322)\boldsymbol{x}_0=\begin{pmatrix}\frac{21}{11}\\\frac{43}{22}\\\frac{3}{22}\end{pmatrix}x0=11212243223，对应的拉格朗日乘子λ0=(2911−1511)\boldsymbol{\lambda}_0=\begin{pmatrix}\frac{29}{11}\\-\frac{15}{11}\end{pmatrix}λ0=(1129−1115)。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！
北京人形旗下天工造物具身智能开源社区，聚焦具身天工与慧思开物两大平台
更多推荐
5
0- 0
已为社区贡献2条内容
所有评论(0)