最优化方法Python计算：二次规划的拉格朗日算法

目标函数为二次式，约束条件为线性式的最优化问题称为二次规划。其一般形式为
$$\{\begin{array}{l}\text{minimize}\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}+{\mathbf{c}}^{\top }\mathbf{x}\\ \text{s.t. }{\mathbf{A}}_{eq}\mathbf{x}-{\mathbf{b}}_{eq}=\mathbf{o}\\ {\mathbf{A}}_{iq}\mathbf{x}-{\mathbf{b}}_{iq}\ge \mathbf{o}\end{array}.$$
其中，$\mathbf{H}\in {\text{R}}^{n\times n}$对称，$\mathbf{c}\in {\text{R}}^n$，${\mathbf{A}}_{eq}\in {\text{R}}^{l\times n}$，${\mathbf{b}}_{eq}\in {\text{R}}^l$，${\mathbf{A}}_{iq}\in {\text{R}}^{m\times n}$，${\mathbf{b}}_{iq}\in {\text{R}}^m$。
仅含等式约束的二次规划形如
$$\{\begin{array}{l}\text{minimize}\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}+{\mathbf{c}}^{\top }\mathbf{x}\\ \text{s.t. }\mathbf{Ax}-\mathbf{b}=\mathbf{o}\end{array}.(1)$$
假定$\mathbf{H}$对称正定，$\mathbf{A}\in {\text{R}}^{l\times n}$，rank$\mathbf{A}=l$。正定二次式$\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}+{\mathbf{c}}^{\top }\mathbf{x}$在凸集 Ω = { x ∣ A x − b = o } \Omega=\{\boldsymbol{x}|\boldsymbol{Ax}-\boldsymbol{b}=\boldsymbol{o}\} Ω={x∣Ax−b=o}上有唯一满足必要条件的KKT点$\left(\begin{array}{c}{\mathbf{x}}_0\\ {\mathbf{\lambda }}_0\end{array}\right)$。为算得该KKT点，写出问题的拉格朗日函数
$$L(\mathbf{x},\mathbf{\lambda })=\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}+{\mathbf{c}}^{\top }\mathbf{x}-{\mathbf{\lambda }}^{\top }(\mathbf{Ax}-\mathbf{b}).$$
关于$\mathbf{x}$，$\mathbf{\lambda }$的梯度为
$$\begin{array}{l}{\nabla }_{\mathbf{x}}L(\mathbf{x},\mathbf{\lambda })=\mathbf{Hx}+\mathbf{c}-{\mathbf{A}}^{\top }\mathbf{\lambda }\\ {\nabla }_{\mathbf{\lambda }}L(\mathbf{x},\mathbf{\lambda })=-\mathbf{Ax}+\mathbf{b}\end{array}$$
令${\nabla }_{\mathbf{x}}L(\mathbf{x},\mathbf{\lambda })=\mathbf{o}$且${\nabla }_{\mathbf{\lambda }}L(\mathbf{x},\mathbf{\lambda })=\mathbf{o}$，得线性方程组
$$\{\begin{array}{l}\mathbf{Hx}+\mathbf{c}-{\mathbf{A}}^{\top }\mathbf{\lambda }=\mathbf{o}\\ -\mathbf{Ax}+\mathbf{b}=\mathbf{o}\end{array},$$
等价地表示为
$$\left(\begin{array}{cc}\mathbf{H}& -{\mathbf{A}}^{\top }\\ -\mathbf{A}& \mathbf{O}\end{array}\right)\left(\begin{array}{c}\mathbf{x}\\ \mathbf{\lambda }\end{array}\right)=\left(\begin{array}{c}-\mathbf{c}\\ -\mathbf{b}\end{array}\right).$$
系数矩阵$\left(\begin{array}{cc}\mathbf{H}& -{\mathbf{A}}^{\top }\\ -\mathbf{A}& \mathbf{O}\end{array}\right)$称为拉格朗日矩阵。由$\mathbf{H}$对称正定且rank$\mathbf{A}=l$的假设，拉格朗日矩阵可逆，设
$${\left(\begin{array}{cc}\mathbf{H}& -{\mathbf{A}}^{\top }\\ -\mathbf{A}& \mathbf{O}\end{array}\right)}^{-1}=\left(\begin{array}{cc}\mathbf{Q}& -{\mathbf{R}}^{\top }\\ -\mathbf{R}& \mathbf{S}\end{array}\right),$$
根据
$$\left(\begin{array}{cc}\mathbf{H}& -{\mathbf{A}}^{\top }\\ -\mathbf{A}& \mathbf{O}\end{array}\right)\left(\begin{array}{cc}\mathbf{Q}& -{\mathbf{R}}^{\top }\\ -\mathbf{R}& \mathbf{S}\end{array}\right)=\mathbf{I}$$
算得
$$\{\begin{array}{l}\mathbf{HQ}+{\mathbf{A}}^{\top }\mathbf{R}=\mathbf{I}\\ -\mathbf{H}{\mathbf{R}}^{\top }-{\mathbf{A}}^{\top }\mathbf{S}=\mathbf{O}\\ -\mathbf{AQ}=\mathbf{O}\\ {\mathbf{AR}}^{\top }=\mathbf{I}\end{array}$$
由于$\mathbf{A}$行满秩，故${\mathbf{AH}}^{-1}{\mathbf{A}}^{\top }$可逆。$({\mathbf{AH}}^{-1}{\mathbf{A}}^{\top }{)}^{-1}{\mathbf{AH}}^{-1}$是${\mathbf{A}}^{\top }$的伪逆。解上列连立式得
$$\{\begin{array}{l}\mathbf{S}=-({\mathbf{AH}}^{-1}{\mathbf{A}}^{\top }{)}^{-1}\\ \mathbf{R}=-{\mathbf{SAH}}^{-1}\\ \mathbf{Q}={\mathbf{H}}^{-1}-{\mathbf{H}}^{-1}{\mathbf{A}}^{\top }\mathbf{R}\end{array}$$
于是，二次规划(1)的KKT点
$$\left(\begin{array}{c}{\mathbf{x}}_0\\ {\mathbf{\lambda }}_0\end{array}\right)=\left(\begin{array}{cc}\mathbf{Q}& -{\mathbf{R}}^{\top }\\ -\mathbf{R}& \mathbf{S}\end{array}\right)\left(\begin{array}{c}-\mathbf{c}\\ -\mathbf{b}\end{array}\right)=\left(\begin{array}{c}-\mathbf{Qc}+{\mathbf{R}}^{\top }\mathbf{b}\\ \mathbf{Rc}-\mathbf{Sb}\end{array}\right).$$
下列代码实现求解等式约束二次规划(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)中的正定矩阵$\mathbf{H}$，行满秩阵$\mathbf{A}$，向量$\mathbf{b}$和$\mathbf{c}$。
函数体内的第3行调用numpy.linalg的inv函数计算$\mathbf{H}$的逆阵${\mathbf{H}}^{-1}$，赋予H1。第4~6行分别计算
$$\begin{array}{l}\mathbf{S}=-({\mathbf{AH}}^{-1}{\mathbf{A}}^{\top }{)}^{-1}\\ \mathbf{R}=-{\mathbf{SAH}}^{-1}\\ \mathbf{Q}={\mathbf{H}}^{-1}-{\mathbf{H}}^{-1}{\mathbf{A}}^{\top }\mathbf{R}\end{array}$$
并赋予S，R和Q。第7、8行分别计算最优解和对应的拉格朗日乘子
$$\begin{array}{l}{\mathbf{x}}_0=-\mathbf{Qc}+{\mathbf{R}}^{\top }\mathbf{b}\\ {\mathbf{\lambda }}_0=\mathbf{Rc}-\mathbf{Sb}\end{array}$$
并赋予x0和lamd0。
例1用qlag函数求解下列二次规划
$$\{\begin{array}{l}\text{minimize}{x}_1^2+2x_2^2+x_3^2-2x_1{x}_2+x_3\\ \text{s.t. }{x}_1+x_2+x_3=4\\ 2x_1-x_2+x_3=2\end{array}.$$
解:本问题中，
$$\mathbf{H}=\left(\begin{array}{ccc}2& -2& 0\\ -2& 4& 0\\ 0& 0& 2\end{array}\right),\mathbf{A}=\left(\begin{array}{ccc}1& 1& 1\\ 2& -1& 1\end{array}\right),\mathbf{b}=\left(\begin{array}{c}4\\ 2\end{array}\right),\mathbf{c}=\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)$$
下列代码利用这些数据进行计算。

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]))

意味着最优解${\mathbf{x}}_0=\left(\begin{array}{c}\frac{21}{11}\\ \frac{43}{22}\\ \frac{3}{22}\end{array}\right)$，对应的拉格朗日乘子${\mathbf{\lambda }}_0=\left(\begin{array}{c}\frac{29}{11}\\ -\frac{15}{11}\end{array}\right)$。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！