本文主要是介绍最优化方法Python计算：二次规划的拉格朗日算法，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
目标函数为二次式，约束条件为线性式的最优化问题称为二次规划。其一般形式为
{ 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)。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！
这篇关于最优化方法Python计算：二次规划的拉格朗日算法的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
