精华内容
下载资源
问答
  • 此提交是应请求编写的 - 作为处理线性最小二乘问题的工具,受可能存在秩不足的线性等式约束。 (当然,它也可以处理满秩约束的问题。)如果出现秩不足约束系统,它会测试约束的一致性。 我向 LSE 添加了一些其他...
  • 求解带线性等式约束的界约束优化问题的有效集算法,闫秀娟,王永丽,许多实际应用问题可以被表述为带界约束和单一线性等式约束的优化问题。本文介绍了一种新的求解此问题的方法。该方法利用零空间法
  • 等式约束下凸二次规划问题的改进拟Newton算法,王建芳,杨晓光,提出了拟Newton法求解凸二次规划问题的改进拟Newton法,对于等式约束下凸二次规划问题利用增广Lagrange函数将该约束问题转化为无约束问�
  • 线性等式约束优化的大规模信赖域方法 MATLAB代码 该项目包含文章《具有低维线性等式约束的大型拟牛顿信任区域方法》中算法的实现,JJ Brust,RF Marcia,CG Petra,2018年。 组成该项目的文件夹为: main:算法1和...
  • 根据广义乘子法的思想,将等式约束的凸二次规划转化为无约束问题,再利用正交校正共扼梯度法来求解,得到等式约束严格凸二次规划的新算法,不用求逆矩阵,这样可用来解大规模稀疏问题,数值结果表明:在微机486/33上就能解...
  • 等式和等式约束

    2021-04-18 16:19:11
    创建一个名为 x 的 4×6 优化...创建 x 的每行总和为 1 的等式。constrsum = sum(x,2) == 1constrsum =4x1 Linear OptimizationEquality array with properties:IndexNames: {{} {}}Variables: [1x1 struct] cont...

    创建一个名为 x 的 4×6 优化变量矩阵。

    x = optimvar('x',4,6);

    创建 x 的每行总和为 1 的等式。

    constrsum = sum(x,2) == 1

    constrsum =

    4x1 Linear OptimizationEquality array with properties:

    IndexNames: {{} {}}

    Variables: [1x1 struct] containing 1 OptimizationVariable

    See equality formulation with show.

    查看等式。

    show(constrsum)

    (1, 1)

    x(1, 1) + x(1, 2) + x(1, 3) + x(1, 4) + x(1, 5) + x(1, 6) == 1

    (2, 1)

    x(2, 1) + x(2, 2) + x(2, 3) + x(2, 4) + x(2, 5) + x(2, 6) == 1

    (3, 1)

    x(3, 1) + x(3, 2) + x(3, 3) + x(3, 4) + x(3, 5) + x(3, 6) == 1

    (4, 1)

    x(4, 1) + x(4, 2) + x(4, 3) + x(4, 4) + x(4, 5) + x(4, 6) == 1

    要在优化问题中包含等式,请使用圆点表示法将 Constraints 属性设置为 constrsum。

    prob = optimproblem;

    prob.Constraints.constrsum = constrsum

    prob =

    OptimizationProblem with properties:

    Description: ''

    ObjectiveSense: 'minimize'

    Variables: [1x1 struct] containing 1 OptimizationVariable

    Objective: [0x0 OptimizationExpression]

    Constraints: [1x1 struct] containing 1 OptimizationConstraint

    See problem formulation with show.

    同样,要在方程问题中包含等式,请使用圆点表示法将 Constraints 属性设置为 constrsum。

    eqnprob = eqnproblem;

    eqnprob.Equations.constrsum = constrsum

    eqnprob =

    EquationProblem with properties:

    Description: ''

    Variables: [1x1 struct] containing 1 OptimizationVariable

    Equations: [1x1 struct] containing 1 OptimizationEquality

    See problem formulation with show.

    展开全文
  • 等式约束的优化问题的求解方法 对数障碍 log barrier\text{log barrier}log barrier 首先, 介绍一下 log barrier\text{log barrier}log barrier 方法: min⁡f0(x)s.t.fi(x)≤0i=1⋯mAx=b⇔min...

    有等式约束的优化问题的求解方法

    对数障碍 log barrier \text{log barrier} log barrier

    首先, 介绍一下 log barrier \text{log barrier} log barrier 方法:

    min ⁡ f 0 ( x ) s.t. f i ( x ) ≤ 0 i = 1 ⋯ m A x = b ⇔ min ⁡ f 0 ( x ) + ∑ i = 1 m − ( 1 / t ) log ⁡ ( − f i ( x ) ) s.t. A x = b \begin{array}{lll} \min & f_{0}(x) &\\ \text{s.t.} & f_{i}(x) \leq 0 & i=1 \cdots m \\ & Ax = b & \\ \Leftrightarrow & & \\ \min & f_{0}(x)+\sum_{i=1}^{m}-(1 / t) \log \left(-f_{i}(x)\right) \\ \text{s.t.} & A x=b \end{array} mins.t.mins.t.f0(x)fi(x)0Ax=bf0(x)+i=1m(1/t)log(fi(x))Ax=bi=1m

    这是一个惩罚的思想. 具体来说就是 f i ( x ) → 0 f_{i}(x) \rightarrow 0 fi(x)0 时, p ∗ → + ∞ p^{*} \rightarrow+\infty p+, 会让目标函数根本无法得到最优解, 就好像有一个“棚栏”一样, 阻止 f i ( x ) f_{i}(x) fi(x) 趋近于 0, 让它只能小于0. 基于 log barrier \text{log barrier} log barrier 方法, 我们之后所有的约束都只讨论等式约束, 不等式约束用 log barrier \text{log barrier} log barrier 去掉. 即所有的凸问题都形如:
    min ⁡ f 0 ( x ) s.t.  A x = b \begin{array}{lll} \min & f_{0}(x) & \\ \text {s.t. } & Ax=b & \end{array} mins.t. f0(x)Ax=b

    由于 − ( 1 / t ) log ⁡ ( − u ) -(1 / t) \log (-u) (1/t)log(u) u u u 的单增凸函数, 上式中的目标函数是可微凸函数. 假定恰当的闭性条件成立,则可用 Newton 方法求解该问题.
    我们将函数
    ϕ ( x ) = − ∑ i = 1 m log ⁡ ( − f i ( x ) ) \phi(x)=-\sum_{i=1}^{m} \log \left(-f_{i}(x)\right) ϕ(x)=i=1mlog(fi(x))
    dom ⁡ ϕ = { x ∈ R n ∣ f i ( x ) < 0 , i = 1 , ⋯   , m } \operatorname{dom} \phi=\left\{x \in \mathbf{R}^{n} \mid f_{i}(x)<0, i=1, \cdots, m\right\} domϕ={xRnfi(x)<0,i=1,,m} 称为原问题的对数障碍函数或对数障碍. 其定义域是满足原问题的严格不等式约束的点集. 不管正参数 t t t 取什么值, 对于任意 i i i, 当 f i ( x ) → 0 f_{i}(x) \rightarrow 0 fi(x)0 时,对数障碍函数将趋于无穷大. 并且, 我们可以知道当 t t t 越来越大时, ∑ i = 1 m − ( 1 / t ) log ⁡ ( − f i ( x ) ) \sum_{i=1}^{m}-(1 / t) \log \left(-f_{i}(x)\right) i=1m(1/t)log(fi(x)) 会非常近似于指示函数( indicator function \text{indicator function} indicator function), 即当 f i ( x ) < 0 f_{i}(x)<0 fi(x)<0 时, ∑ i = 1 m − ( 1 / t ) log ⁡ ( − f i ( x ) ) = 0 \sum_{i=1}^{m}-(1 / t) \log \left(-f_{i}(x)\right) = 0 i=1m(1/t)log(fi(x))=0, 当 f i ( x ) = 0 f_{i}(x) = 0 fi(x)=0 时, ∑ i = 1 m − ( 1 / t ) log ⁡ ( − f i ( x ) ) = + ∞ \sum_{i=1}^{m}-(1 / t) \log \left(-f_{i}(x)\right) = +\infty i=1m(1/t)log(fi(x))=+

    有等式约束的优化问题定义

    min ⁡ f ( x ) s.t. A x = b (eq0) \begin{array}{ll} \min &f(x) \\ \text{s.t.} &A x = b \end{array} \tag{eq0} mins.t.f(x)Ax=b(eq0)

    计算它的 KKT \text{KKT} KKT 条件:
    x ∗ 是 最 优 解 , v ∗ 是 拉 格 朗 日 乘 子 { A x ∗ = b primal feasibility ∇ f ( x ∗ ) + A T v ∗ = 0 stationarity x^{*} 是最优解, v^{*} 是拉格朗日乘子 \\ \left\{ \begin{array}{l} A x^{*}=b \quad \text{primal feasibility} \\ \nabla f\left(x^{*}\right)+A^{T} v^{*}=0 \quad \text{stationarity} \end{array} \right. x,v{Ax=bprimal feasibilityf(x)+ATv=0stationarity

    解约束优化问题就不得不用到 KKT \text{KKT} KKT 条件了,解 KKT \text{KKT} KKT 条件(如果可解的话), 那么相当于一次性 把原问题和对偶问题的最优解都解出来了. 目前有一个从 KKT \text{KKT} KKT 条件出发得出的算法, 就叫做拉格朗日乘子法]

    当目标函数 f ( x ) f(x) f(x) 是二次函数时, KKT \text{KKT} KKT 条件是线性方程组

    当目标函数 f ( x ) f(x) f(x) 是二次函数时, KKT \text{KKT} KKT 条件是线性方程组.例如
    min 1 2 x T P x + q T x + r s.t. A x = b \begin{array}{ll} \text{min} &\frac{1}{2} x^{T} P x+q^{T} x+r \\ \text{s.t.} &A x = b \end{array} mins.t.21xTPx+qTx+rAx=b

    其中 P ∈ S + n , A P \in S_{+}^{n}, A PS+n,A 一般是不可逆的(这样 x x x 才会是有一个解集的), 否则可逆矩阵的方程组 有唯一解就谈不上要优化了.
    计算这个问题的 KKT \text{KKT} KKT 条件为:
    { A x ∗ = b P x ∗ + q + A T v ∗ = 0 ⇔ [ P A T A 0 ] ⋅ [ x ∗ v ∗ ] = [ − q b ] \left\{\begin{array}{ll}A x^{*} =b \\ P x^{*} +q+A^{T} v^{*}=0\end{array}\right. \Leftrightarrow \left[\begin{array}{cc}P & A^{T} \\ A & 0\end{array}\right] \cdot\left[\begin{array}{l}x^{*} \\ v^{*}\end{array}\right]=\left[\begin{array}{c}-q \\ b\end{array}\right] {Ax=bPx+q+ATv=0[PAAT0][xv]=[qb]
    我们会发现,这时候 KKT \text{KKT} KKT 条件也是一个线性的方程组, 主要原因就是 ∇ f ( x ∗ ) \nabla f\left(x^{*}\right) f(x) 是线性的, 这 个时候就可以用Gauss-Seidel等解线性方程组的方法求解最优值了. 这样的情况不需要什 么下降算法, 本质就是求解线性方程组. 当然, 你也可以转化为求无约束优化的对偶问题.

    当目标函数 f ( x ) f(x) f(x) 是高次函数时, KKT \text{KKT} KKT 条件不是线性的方程组

    ∇ f ( x ∗ ) \nabla f\left(x^{*}\right) f(x) 不是一个线性函数时,则有
    { A x ∗ = b ∇ f ( x ∗ ) + A T v ∗ = 0 \left\{\begin{array}{l}A x^{*}=b \\ \nabla f\left(x^{*}\right)+A^{T} v^{*}=0\end{array}\right. {Ax=bf(x)+ATv=0

    如果 ∇ f ( x ∗ ) \nabla f\left(x^{*}\right) f(x) 不是一个关于 x ∗ x^{*} x 的线性函数时, 自然会想到要找一个办法把它变成线性函数. 例如使用近似的思想,如果函数二次可微的话, 可以对函数进行二阶展开.
    假设现在有一个可行点 x k x^{k} xk, x k x^{k} xk 可以使上式 ( e q 0 ) (eq0) (eq0) 成立,注意, x k x^{k} xk 只是可行解,而不一定是最优解.
    现在我们构建一个问题 ( e q 1 ) (eq1) (eq1) ( e q 1 ) (eq1) (eq1): 在 x k x^{k} xk 的邻域附近寻找一个新的点 x k + d x^{k}+d xk+d,使得 f ( x k + d ) f(x^{k}+d) f(xk+d) 尽可能地小. 问题 ( e q 1 ) (eq1) (eq1) ( e q 1 ) (eq1) (eq1) 一定有解, 最坏情况下, d = 0 d=0 d=0, ( e q 1 ) (eq1) (eq1) 的解与 ( e q 0 ) (eq0) (eq0) 的解相同. 令 x k + 1 = x k + d x^{k+1} = x^{k} + d xk+1=xk+d, 然后我们再将 x k + 1 x^{k+1} xk+1 当做 x k x^{k} xk, 再在 x k + 1 x^{k+1} xk+1 邻域内寻找问题 ( e q 1 ) (eq1) (eq1) ( e q 1 ) (eq1) (eq1) 的解, 如此往复, 则我们就可以得到原问题 ( e q 0 ) (eq0) (eq0) 的解(或者近似解).

    min ⁡ d f ( x k + d ) s.t. A ( x k + d ) = b ⇔ min ⁡ d f ( x k + d ) s.t. A d = 0 因 为 A x k = b , A ( x k + d ) = b 所 以 A d = 0 (eq1) \begin{array}{ll} \underset{d}{\min} & f\left(x^{k}+d\right) \\ \text{s.t.} & A\left(x^{k}+d\right)=b \end{array} \Leftrightarrow \begin{array}{ll} \underset{d}{\min} & f\left(x^{k}+d\right) \\ \text {s.t.} &A d=0 \\ \end{array} \\ 因为A x^{k}=b, A\left(x^{k}+d\right)=b 所以A d=0 \\ \tag{eq1} dmins.t.f(xk+d)A(xk+d)=bdmins.t.f(xk+d)Ad=0Axk=b,A(xk+d)=bAd=0(eq1)
    则我们现在把原问题(求解全局最优解)等价为新的优化问题(在 x k x^{k} xk的邻域附近寻找一个新的点 x k + d x^{k}+d xk+d,使得 f ( x k + d ) f(x^{k}+d) f(xk+d) 尽可能地小).
    这是一个关于 d d d 的新的优化问题. 为了让原问题 KKT \text{KKT} KKT 线性, 先对问题 ( e q 1 ) (eq1) (eq1)目标函数做二阶泰勒展开
    f ( x k + d ) = f ( x k ) + ∇ f ( x k ) T d + 1 2 d T ( ∇ 2 f ( x k ) ) d f\left(x^{k}+d\right)=f\left(x^{k}\right)+\nabla f\left(x^{k}\right)^{T} d+\frac{1}{2} d^{T}\left(\nabla^{2} f\left(x^{k}\right)\right) d f(xk+d)=f(xk)+f(xk)Td+21dT(2f(xk))d
    则问题 ( e q 1 ) (eq1) (eq1)又被近似等价为问题 ( e q 2 ) (eq2) (eq2)
    min ⁡ d f ( x k ) + ∇ f ( x k ) T d + 1 2 d T ∇ 2 f ( x k ) T d s.t. A d = 0 (eq2) \begin{array}{ll} \underset{d}{\min} & f\left(x^{k}\right)+\nabla f\left(x^{k}\right)^{T} d+\frac{1}{2} d^{T} \nabla^{2} f\left(x^{k}\right)^{T} d \\ \text {s.t.} & A d=0 \end{array} \\ \tag{eq2} dmins.t.f(xk)+f(xk)Td+21dT2f(xk)TdAd=0(eq2)
    求解这个问题 ( e q 2 ) (eq2) (eq2) KKT \text{KKT} KKT 条件为:
    [ ∇ 2 f ( x k ) A T A 0 ] ⋅ [ d ∗ v ∗ ] = [ − ∇ f ( x k ) 0 ] \left[\begin{array}{cc}\nabla^{2} f\left(x^{k}\right) & A^{T} \\ A & 0\end{array}\right] \cdot\left[\begin{array}{l}d^{*} \\ v^{*}\end{array}\right]=\left[\begin{array}{c}-\nabla f\left(x^{k}\right) \\ 0\end{array}\right] [2f(xk)AAT0][dv]=[f(xk)0]
    便又回到解线性方程组了. 但是这里要注意的是, 我们这样近似解出来的是 d k d^{k} dk, 也就是用上一次的迭代点沿着 d d d 方向展开了得到的 d k d^{k} dk, 这一个过程我们进行了二阶展开这样的近似操作, 并不是求解出了原问题的最优解 但是我们可以发现, 解出来的 d k d^{k} dk 满足 A d k = 0 A d^{k}=0 Adk=0, 也就是说 x k + d k x^{k}+d^{k} xk+dk 起码是可行点,满足 A ( x k + d k ) = b A\left(x^{k}+d^{k}\right)=b A(xk+dk)=b .
    到这里我们就会思考,可不可以用这个方向去下降呢? 答案是最好先别. 因为 x k + d k x^{k}+d^{k} xk+dk 只 是可行点, 不一定相对于 x k x^{k} xk 是下降点. 所以我们可以给 d k d^{k} dk 乘上一个步长, 对步长做一次线搜索,再来当作下降方向来用. 由此我们得到带等式约束的牛顿法:
    { [ ∇ 2 f ( x k ) A T A 0 ] ⋅ [ d ∗ v ∗ ] = [ − ∇ f ( x k ) 0 ] d k = d ∗ α k = arg ⁡ min ⁡ α ≥ 0 f ( x k + α d k ) x k + 1 = x k + α k d k ,    x k + 1 始 终 满 足 A x k + 1 = b \left\{ \begin{array}{l} \left[\begin{array}{cc}\nabla^{2} f\left(x^{k}\right) & A^{T} \\ A & 0\end{array}\right] \cdot\left[\begin{array}{l}d^{*} \\ v^{*}\end{array}\right]=\left[\begin{array}{c}-\nabla f\left(x^{k}\right) \\ 0\end{array}\right] \\ d^{k}=d^{*} \\ \alpha^{k}=\arg \underset{\alpha \geq 0}{\min} f \left(x^{k}+\alpha d^{k}\right) \\ x^{k+1}=x^{k}+\alpha^{k} d^{k}, \; x^{k+1}始终满足A x^{k+1} = b \\ \end{array} \right. [2f(xk)AAT0][dv]=[f(xk)0]dk=dαk=argα0minf(xk+αdk)xk+1=xk+αkdk,xk+1Axk+1=b
    所以,对于带约束优化问题, KKT \text{KKT} KKT 条件还是蛮重要的

    拉格朗日乘子法Method of Lagrangian Multiplier, 拉格朗日法Lagrangian Method

    拉格朗日乘子法/拉格朗日法是通过某种方法将原变量和对偶变量看成一个大的变量集合体, 然后确定一个优化方向, 同时优化原变量和对偶变量

    { x k + 1 = x k − α k ( ∇ f ( x k ) + A T v k ) v k + 1 = v k + α k ( A x k − b ) ⇔ [ x k + 1 v k + 1 ] = [ x k v k ] + α k [ − ( ∇ f ( x k ) + A T v k ) A x k − b ] \left\{\begin{array}{l}x^{k + 1}=x^{k}-\alpha^{k}\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right) \\ v^{k+1}=v^{k}+\alpha^{k}\left(A x^{k}-b\right)\end{array}\right. \\ \Leftrightarrow \\ \left [\begin{array}{c} x^{k + 1} \\ v^{k+1} \end{array}\right ] = \left [\begin{array}{c} x^{k} \\ v^{k} \end{array}\right]+ \alpha^{k} \left [\begin{array}{c} -\left( \nabla f (x^{k})+A^{T} v^{k} \right) \\ A x^{k}-b \end{array}\right] {xk+1=xkαk(f(xk)+ATvk)vk+1=vk+αk(Axkb)[xk+1vk+1]=[xkvk]+αk[(f(xk)+ATvk)Axkb]

    第二个式子乘子的迭代部分, x k , v k x^{k},v^{k} xk,vk 肯定是可以换成上一步已经计算好的 x k + 1 , v k + 1 x^{k+1},v^{k+1} xk+1,vk+1. 很明显,拉格朗日寻找下一个更新点的方法如上式,找到一个方向 [ − ( ∇ f ( x k ) + A T v k ) A x k − b ] \left [\begin{array}{c} -\left( \nabla f (x^{k})+A^{T} v^{k} \right) \\ A x^{k}-b \end{array}\right] [(f(xk)+ATvk)Axkb], 再寻找合适的步长 α k \alpha^{k} αk, 然后进行更新.

    从鞍点角度理解拉格朗日法

    因为
    L ( x , v ) = f ( x ) + v T ( A x − b ) L(x, v)=f(x)+v^{T}(A x-b) L(x,v)=f(x)+vT(Axb)
    所以鞍点
    { ( x ∗ , v ∗ ) = arg ⁡ max ⁡ v    min ⁡ x    L ( x , v ) ( x ∗ , v ∗ ) = arg ⁡ min ⁡ x    max ⁡ v    L ( x , v ) \begin{aligned} \left \{ \begin{array}{ll} (x^{*}, v^{*}) &=\arg \underset{v}{\max} \; \underset{x}{\min} \; L(x, v) \\ (x^{*}, v^{*}) &=\arg \underset{x}{\min} \; \underset{v}{\max} \; L(x, v) \end{array} \right. \end{aligned} {(x,v)(x,v)=argvmaxxminL(x,v)=argxminvmaxL(x,v)
    我们知道,如果强对偶性成立,那么鞍点就是最优点. 此时
    x ∗ = arg ⁡ min ⁡ x L ( x , v ∗ ) ( 已 知 v ∗ 求 x ∗ ) v ∗ = arg ⁡ max ⁡ v L ( x ∗ , v ) ( 已 知 x ∗ 求 v ∗ ) 上 述 求 x ∗ , v ∗ 的 过 程 是 无 约 束 优 化 问 题 求 极 值 的 方 法 , 可 以 使 用 梯 度 下 降 或 者 梯 度 上 升 的 算 法 . x^{*}=\arg \underset{x}{\min} L\left(x, v^{*}\right)\left(\right. 已知 v^{*} 求 \left.x^{*}\right) \\ v^{*}=\arg \underset{v}{\max} L\left(x^{*}, v\right)\left(\right.已知 x^{*}求 \left.v^{*}\right) \\ 上述求 x^{*}, v^{*}的过程是无约束优化问题求极值的方法,可以使用梯度下降或者梯度上升的算法. x=argxminL(x,v)(vx)v=argvmaxL(x,v)(xv)x,v,使.

    如果已知某一个分量的最优点, 求另一个分量的最优点本质上就是求拉格朗日函数对于单个变量的最优解. 接下来分类讨论下:

    1. 当已知 v ∗ v^{*} v x ∗ x^{*} x 时:
      用梯度下降法迭代求解 L ( x , v ∗ ) L\left(x, v^{*}\right) L(x,v) . 因为 − ∇ L ( x , v ∗ ) = − ∇ f ( x ) − A T v ∗ -\nabla L\left(x, v^{*}\right)=-\nabla f(x)-A^{T} v^{*} L(x,v)=f(x)ATv, 所以我们的迭代算法为 x k + 1 = x k + α k ( − ∇ f ( x k ) − A T v ∗ ) x^{k+1}=x^{k}+\alpha^{k}\left(-\nabla f\left(x^{k}\right)-A^{T} v^{*}\right) xk+1=xk+αk(f(xk)ATv) 但是我们不可能真的已知 v ∗ v^{*} v 是什么, 所以我们用 v k v^{k} vk 近似. 这便是拉格朗日法的关于 x x x 的 迭代部分: x k + 1 = x k − α k ( ∇ f ( x k ) + A T v k ) x^{k+1}=x^{k}-\alpha^{k}\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right) xk+1=xkαk(f(xk)+ATvk)
    2. 当已知 x ∗ x^{*} x v ∗ v^{*} v 时:
      因为我们求极大, 所以上面的迭代方向就是梯度方向,同样地 x ∗ x^{*} x x k x^{k} xk 替代. v k + 1 = v k + α k ( A x k − b ) v^{k+1}=v^{k}+\alpha^{k}\left(A x^{k}-b\right) vk+1=vk+αk(Axkb)
      拉格朗日法实际上是在求 L ( x , v ) L(x, v) L(x,v) 的鞍点, 来得出原问题的最优值.

    从罚函数的角度理解拉格朗日法(最好理解的方式)

    我们把带等式约束的原问题的 KKT \text{KKT} KKT 条件:
    { A x ∗ = b ∇ f ( x ∗ ) + A T v ∗ = 0 \left\{\begin{array}{l}A x^{*}=b \\ \nabla f\left(x^{*}\right)+A^{T} v^{*}=0\end{array}\right. \\ {Ax=bf(x)+ATv=0
    改写成一个罚函数的样子, 去当作一个新的优化问题来解:
    min ⁡ P ( x , v ) = 1 2 ∥ A x − b ∥ 2 2 + 1 2 ∥ ∇ f ( x ) + A T v ∥ 2 2 ( 一 般 为 非 凸 问 题 ) (eq3) \min P(x, v)=\frac{1}{2}\|A x-b\|_{2}^{2}+\frac{1}{2}\left\|\nabla f(x)+A^{T} v\right\|_{2}^{2} \quad (一般为非凸问题) \\ \tag{eq3} minP(x,v)=21Axb22+21f(x)+ATv22()(eq3)
    特点: 这个新的无约束优化问题的解就是满足 KKT \text{KKT} KKT 条件的点. 某种意义上说,相当于把原问题转化成了无约束优化问题, 我们可以去求这个新问题的目标函数的负梯度方向,然后梯度下降求最优解:

    − ∇ P ( x k , v k ) = [ − ∇ x P ( x k , v k ) , − ∇ v P ( x k , v k ) ] = − [ A T ( A x k − b ) + ∇ 2 f ( x k ) ( ∇ f ( x k ) + A T v k ) A ( ∇ f ( x k ) + A T v k ) ] \begin{array}{ll} &-\nabla P\left(x^{k}, v^{k}\right) \\ &=\left[\begin{array}{l}-\nabla_{x} P\left(x^{k}, v^{k}\right), -\nabla_{v} P\left(x^{k}, v^{k}\right)\end{array}\right] \\ &=-\left[ \begin{array}{c} A^{T}\left(A x^{k}-b\right)+\nabla^{2} f\left(x^{k}\right)\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right) \\ A\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right) \end{array} \right] \end{array}\\ P(xk,vk)=[xP(xk,vk),vP(xk,vk)]=[AT(Axkb)+2f(xk)(f(xk)+ATvk)A(f(xk)+ATvk)]

    拉格朗日方法中迭代方向的正确性证明

    我们选择 − ∇ P ( x k , v k ) -\nabla P\left(x^{k}, v^{k}\right) P(xk,vk) 为梯度方向,如果拉格朗日法所确定的迭代方向 d k d^{k} dk
    − ∇ P ( x k , v k ) -\nabla P\left(x^{k}, v^{k}\right) P(xk,vk) 夹角小于 π 2 \frac{\pi}{2} 2π,那么我们认为迭代方向 d k d^{k} dk 能够使函数值下降.
    d k = [ − ( ∇ f ( x k ) + A T v k ) A x k − b ] d^{k}=\left [ \begin{array}{c} -\left( \nabla f (x^{k})+A^{T} v^{k} \right) \\ A x^{k}-b \end{array} \right] dk=[(f(xk)+ATvk)Axkb]

    ( d k ) T ( − ∇ P ( x k , v k ) ) = ( ∇ f ( x k ) + A T v k ) T A T ( A x k − b ) + ( ∇ f ( x k ) + A T v k ) T ∇ 2 f ( x k ) ( ∇ f ( x k ) + A T v k ) − ( A T x k − b ) T A ( ∇ f ( x k ) + A T v k ) = ( ∇ f ( x k ) + A T v k ) T ∇ 2 f ( x k ) ( ∇ f ( x k ) + A T v k ) \begin{aligned} &(d^{k})^{T} \left(-\nabla P(x^{k}, v^{k})\right)\\ &=\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right)^{T} A^{T}\left(A x^{k}-b\right) \\ &\qquad + \left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right)^{T} \nabla^{2} f\left(x^{k}\right) \left(\nabla f\left(x^{k}\right)+A^{T} v^{k} \right) \\ &\qquad - (A^{T}x^{k}-b)^{T}A(\nabla f\left(x^{k}\right)+A^{T} v^{k})\\ &=\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right)^{T} \nabla^{2} f\left(x^{k}\right)\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right) \end{aligned} (dk)T(P(xk,vk))=(f(xk)+ATvk)TAT(Axkb)+(f(xk)+ATvk)T2f(xk)(f(xk)+ATvk)(ATxkb)TA(f(xk)+ATvk)=(f(xk)+ATvk)T2f(xk)(f(xk)+ATvk)

    上边这个式子(二次型), 当 ∇ 2 f ( x k ) ≻ 0 \nabla^{2} f\left(x^{k}\right) \succ 0 2f(xk)0 ( ∇ f ( x k ) + A T v k ) \left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right) (f(xk)+ATvk) 不全为 0 时, 会大于等于0 . 即: { ∇ f ( x k ) + A T v k ≠ 0 ∇ 2 f ( x k ) ≻ 0 \left\{\begin{array}{l}\nabla f\left(x^{k}\right)+A^{T} v^{k} \neq 0 \\ \nabla^{2} f\left(x^{k}\right) \succ 0\end{array}\right. {f(xk)+ATvk=02f(xk)0 时, d k d^{k} dk P ( u k , v k ) P\left(u^{k}, v^{k}\right) P(uk,vk) 的一个下降方向(和负梯度方向同向). 而上面这个条件是成立的. 因为原问题凸, 所以Hessian矩阵半正定. 而当随着迭代 x k → x ∗ x^{k} \rightarrow x^{*} xkx 的时候, ∇ 2 f ( x ∗ ) = 0 \nabla^{2} f\left(x^{*}\right)=0 2f(x)=0, 那么大部分情况下 ∇ 2 f ( x k ) \nabla^{2} f\left(x^{k}\right) 2f(xk) 都是正定的. 另外, 只要 v k , x k v^{k}, x^{k} vk,xk 还没有到达最优解, ∇ f ( x k ) + A T v k ≠ 0 \nabla f\left(x^{k}\right)+A^{T} v^{k} \neq 0 f(xk)+ATvk=0 就一直满足, 当 v ∗ → v ∗ , x k → x ∗ v^{*} \rightarrow v^{*}, x^{k} \rightarrow x^{*} vv,xkx 时, ∇ f ( x k ) + A T v k = 0 \nabla f\left(x^{k}\right)+A^{T} v^{k}=0 f(xk)+ATvk=0, 此时 x k + 1 = x k , v k + 1 = v k x^{k+1}=x^{k}, v^{k+1}=v^{k} xk+1=xk,vk+1=vk,算法也就停止更新了. 对拉格朗日法进行理论分析很有意义, 但是这个算法实际运行起来速度很慢, 所以为了提高一下效率, 诞生了增广拉格朗日法.

    增广拉格朗日法

    增广拉格朗日法是交替优化原变量和对偶变量,并以对偶变量优化为主
    增广拉格朗日法相比与拉格朗日法, 会构建一个新的函数, 称为增广拉格朗日函数. 增广拉格朗日函数是由拉格朗日函数加上一个罚项得到的, 记为:
    L c ( x , v ) = f ( x ) + v T ( A x − b ) + c 2 ∥ A x − b ∥ 2 2 L_{c}(x, v)=f(x)+v^{T}(A x-b)+\frac{c}{2}\|A x-b\|_{2}^{2} Lc(x,v)=f(x)+vT(Axb)+2cAxb22
    它实际上是下面这个优化问题 ( e q 4 ) (eq4) (eq4) 的拉格朗日函数:
    min ⁡ f ( x ) + c 2 ∥ A x − b ∥ 2 2 s.t. A x = b (eq4) \begin{array}{ll} \min & f(x)+\frac{c}{2}\|A x-b\|_{2}^{2} \\ \text{s.t.} & A x=b \end{array} \tag{eq4} mins.t.f(x)+2cAxb22Ax=b(eq4)
    也就是把原问题的等式约束变为一个惩罚项再放到目标函数里面去, 得到的优化问题的拉格朗日函数, 被称为原问题拉格朗日函数的增广拉格朗日函数.
    接下来很容易可以验证, 上面 ( e q 4 ) (eq4) (eq4) 这个问题和开头的标准问题 ( e q 0 ) (eq0) (eq0) 的最优解相同, 两个问题是等价的.
    对于问题 ( e q 0 ) (eq0) (eq0), 其拉格朗日函数为: L ( x , v ) = f ( x ) + v T ( A x − b ) L(x, v)=f(x)+v^{T}(A x-b) L(x,v)=f(x)+vT(Axb), 假设拉格朗日函数的最优解为 ( x ∗ , v ∗ ) \left(x^{*}, v^{*}\right) (x,v), 则 L ( x ∗ , v ∗ ) = 0 L\left(x^{*}, v^{*}\right)=0 L(x,v)=0, 即: ∇ x { f ( x ∗ ) + v ∗ T ( A x ∗ − b ) } = 0 \nabla_{x}\left\{f\left(x^{*}\right)+v^{* T}\left(A x^{*}-b\right)\right\}=0 x{f(x)+vT(Axb)}=0
    对于问题 ( e q 4 ) (eq4) (eq4), 其拉格朗日函数为: L c ( x , v ) = f ( x ) + v T ( A x − b ) + c 2 ∥ A x − b ∥ 2 2 L_{c}(x, v)=f(x)+v^{T}(A x-b)+\frac{c}{2}\|A x-b\|_{2}^{2} Lc(x,v)=f(x)+vT(Axb)+2cAxb22
    假设拉格朗日函数的最优解为 ( x 1 ∗ , v 1 ∗ ) \left(x_{1}^{*}, v_{1}^{*}\right) (x1,v1), 则 L c ( x 1 ∗ , v 1 ∗ ) = 0 L_{c}\left(x_{1}^{*}, v_{1}^{*}\right)=0 Lc(x1,v1)=0, 即:
    ∇ x { f ( x 1 ∗ ) + v 1 ∗ T ( A x 1 ∗ − b ) } + c A T ( A x 1 ∗ − b ) = 0 \nabla_{x}\left\{f\left(x_{1}^{*}\right)+v_{1}^{* T}\left(A x_{1}^{*}-b\right)\right\}+c A^{T}\left(A x_{1}^{*}-b\right)=0 x{f(x1)+v1T(Ax1b)}+cAT(Ax1b)=0
    又因为 A x 1 ∗ − b = 0 A x_{1}^{*}-b=0 Ax1b=0, 所以实际上上式为 ∇ x { f ( x 1 ∗ ) + v 1 ∗ T ( A x 1 ∗ − b ) } = 0 \nabla_{x}\left\{f\left(x_{1}^{*}\right)+v_{1}^{* T}\left(A x_{1}^{*}-b\right)\right\}=0 x{f(x1)+v1T(Ax1b)}=0
    这和 ( e q 0 ) (eq0) (eq0) 的结论是一样的. 所以说 ( x ∗ , v ∗ ) = ( x 1 ∗ , v 1 ∗ ) \left(x^{*}, v^{*}\right)=\left(x_{1}^{*}, v_{1}^{*}\right) (x,v)=(x1,v1), 即 ( e q 0 ) ⇔ ( e q 4 ) (eq0) \Leftrightarrow(eq4) (eq0)(eq4) .
    既然两个问题是等价的,也就是它们各自的拉格朗日函数的最优值是等价的. 那么, 自然可 以想到用 ∇ L c ( x , v ) \nabla L_{c}(x, v) Lc(x,v) 去替换掉 ∇ L ( x , v ) \nabla L(x, v) L(x,v) .
    使用拉格朗日法来解决问题 ( e q 4 ) (eq4) (eq4)(问题 ( e q 0 ) (eq0) (eq0)的增广拉格朗日函数就是问题 ( e q 4 ) (eq4) (eq4)的拉格朗日函数)
    { x k + 1 = x k − α k ∇ x L c ( x k , v k ) = x k − α k ( ∇ f ( x k ) + A T v k + c A T ( A x k − b ) ) v k + 1 = v k + α k ∇ v L c ( x k , v k ) = v k + α k ( A x k − b ) \left\{\begin{array}{l}x^{k+1}=x^{k}-\alpha^{k} \nabla_{x} L_{c}\left(x^{k}, v^{k}\right)=x^{k}-\alpha^{k}\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}+c A^{T}\left(A x^{k}-b\right)\right) \\ v^{k+1}=v^{k}+\alpha^{k} \nabla_{v} L_{c}\left(x^{k}, v^{k}\right)=v^{k}+\alpha^{k}\left(A x^{k}-b\right)\end{array}\right. {xk+1=xkαkxLc(xk,vk)=xkαk(f(xk)+ATvk+cAT(Axkb))vk+1=vk+αkvLc(xk,vk)=vk+αk(Axkb)
    但是上面这组迭代式子, 有个东西没有利用起来, 就是第一个式子会更新好的 x k + 1 x^{k+1} xk+1 . 在第二式子计算 v k + 1 v^{k+1} vk+1 时用的还是 x k x^{k} xk 的信息, 这是把它替换成新的 x k + 1 x^{k+1} xk+1, 那么整个算法变为:
    { x k + 1 = x k − α k ( ∇ f ( x k ) + A T v k + c A T ( A x k − b ) ) v k + 1 = v k + α k ( A x k + 1 − b ) \left\{\begin{array}{l}x^{k+1}=x^{k}-\alpha^{k}\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}+c A^{T}\left(A x^{k}-b\right)\right) \\ v^{k+1}=v^{k}+\alpha^{k}\left(A x^{k+1}-b\right)\end{array}\right. {xk+1=xkαk(f(xk)+ATvk+cAT(Axkb))vk+1=vk+αk(Axk+1b)
    如果增广拉格朗日函数 L c ( x , v ) L_{c}(x, v) Lc(x,v) 性质比较理想或者比较容易被优化的话, 不一定非要用梯度下降法, 比如二阶可微也可以用牛顿法, 等等. 所以统一改成用 arg ⁡ min ⁡ \arg \min argmin 的形式来表达上面这 个算法为:
    { x k + 1 = arg ⁡ min ⁡ x L c ( x , v k ) v k + 1 = v k + α k ( A x k + 1 − b ) \left\{\begin{array}{l}x^{k+1}=\arg \underset{x}{\min} L_{c}\left(x, v^{k}\right) \\ v^{k+1}=v^{k}+\alpha^{k}\left(A x^{k+1}-b\right)\end{array}\right. {xk+1=argxminLc(x,vk)vk+1=vk+αk(Axk+1b)
    我们还可以把第二个式子的步长替换为罚参数, 得到最终的增广拉格朗日法:
    { x k + 1 = arg ⁡ min ⁡ x L c ( x , v k ) v k + 1 = v k + c ( A x k + 1 − b ) \left\{\begin{array}{l}x^{k+1}=\arg \underset{x}{\min} L_{c}\left(x, v^{k}\right) \\ v^{k+1}=v^{k}+c\left(A x^{k+1}-b\right)\end{array}\right. {xk+1=argxminLc(x,vk)vk+1=vk+c(Axk+1b)

    因此最终的增广拉格朗日法 Augmented Lagrangian Method \text{Augmented Lagrangian Method} Augmented Lagrangian Method
    x k + 1 = arg ⁡ min ⁡ x L c ( x , v k ) v k + 1 = v k + c ( A x k + 1 − b ) \begin{array}{l} x^{k+1}=\arg \underset{x}{\min} L_{c}\left(x, v^{k}\right) \\ v^{k+1}=v^{k}+c\left(A x^{k+1}-b\right) \end{array} xk+1=argxminLc(x,vk)vk+1=vk+c(Axk+1b)
    步长 α k \alpha^{k} αk 都是要自己找的, 比如精确搜索或非精确搜索.

    注:这里 c c c 要自己选定, 但是选定 c c c 肯定比去线搜索 α \alpha α 简单多了. 特别的,如果选择 c = 1 c=1 c=1, 基本上可以保证增广拉格朗日算法收敘. 而且, 从收敘性分析可以知道, c c c 也可以 是递增序列,不用担心不收敘.

    算法性质:

    (1). 若 v = v ∗ v=v^{*} v=v,则 ∀ c > 0 , x ∗ = arg ⁡ min ⁡ x L c ( x , v ∗ ) \forall c>0, x^{*}=\arg \underset{x}{\min} L_{c}\left(x, v^{*}\right) c>0,x=argxminLc(x,v), 只要 v → v ∗ , c v \rightarrow v^{*}, c vv,c 的值只要大于0就可以得到最优解 x ∗ x^{*} x
    (2). 若 c → + ∞ c \rightarrow+\infty c+,则 ∀ v , x ∗ = arg ⁡ min ⁡ x L c ( x , v ) \forall v, \quad x^{*}=\arg \underset{x}{\min} L_{c}(x, v) v,x=argxminLc(x,v), 步长 c c c 可以选择递增序列 { c k } \left\{c^{k}\right\} {ck}, 此时一定收敛.假设 c = + ∞ c=+\infty c=+

    参考资料

    展开全文
  • 等式约束优化问题等式约束的基本定义消除等式约束基于feasible初值的牛顿法基于infeasible初值的牛顿法通过解原对偶问题解释原理求解KKT系统实例总结 等式约束的基本定义 等式约束问题 min⁡f(x),s.t.Ax=b,A∈Rp...

    等式约束的基本定义

    • 等式约束问题
      min ⁡ f ( x ) , s . t . A x = b , A ∈ R p × n , r a n k ( A ) = p \min{f(x)},s.t. Ax=b,A\isin R^{p\times n},rank(A)=p minf(x),s.t.Ax=bARp×n,rank(A)=p

      f是凸函数,二次可微,假设 p ∗ p^* p是存在且确定的, p < n p<n p<n,即约束条件小于变量维数,可以存在多解,否则可能会无解

    • 等式约束的对偶问题
      拉格朗日函数: L ( x , ν ) = f ( x ) + ν T ( A x − b ) L(x,\nu)=f(x)+\nu^T(Ax-b) L(x,ν)=f(x)+νT(Axb),设最优点 x ∗ x^* x

      由此得到KKT条件方程组 { A x ∗ = b ∇ f ( x ∗ ) + A T ν ∗ = 0 \begin{cases}Ax^*=b \\ \nabla f(x^*)+A^T\nu^*=0 \end{cases} {Ax=bf(x)+ATν=0

      x ∗ 是 最 优 点 , ν ∗ 是 对 偶 最 优 点 , p ∗ = inf ⁡ x { f ( x ) ∣ A x = b } = f ( x ∗ ) x^*是最优点,\nu^*是对偶最优点,p^*=\inf_x\{f(x) | Ax=b\}=f(x^*) xνp=infx{f(x)Ax=b}=f(x)

    • 例子:二次规划问题

      { min ⁡ 1 2 x T P x + q T x + r A x = b    ⟹    { A x = b P x + q + A T ν ∗ = 0 \begin{cases}\min{\frac12x^TPx+q^Tx+r} \\ Ax=b \end{cases}\implies\begin{cases}Ax=b \\ Px+q+A^T\nu^*=0 \end{cases} {min21xTPx+qTx+rAx=b{Ax=bPx+q+ATν=0
      在这里插入图片描述
      前一个左边的矩阵 ( n + p ) × ( n + p ) (n+p)\times(n+p) (n+p)×(n+p),复杂度是 O ( ( n + p ) 3 ) O((n+p)^3) O((n+p)3)
      二次规划问题仅求导一次就可以实现KKT条件,所以

      • 牛顿法(迭代步需要求二阶导):利用二次函数(存在海森矩阵)去逼近普通函数(下文的目标)
      • 梯度下降法:利用一次函数(存在导数)去逼近普通函数

    消除等式约束

    • 通过变量替换实现等式消除

      • 变量替换 x ∗ = F z ∗ + x ^ x^*=Fz^*+\hat{x} x=Fz+x^

        x ^ \hat{x} x^是任意的一个特解,那么可以实现变换 { x ∣ A x = b } = { F z + x ^ ∣ z ∈ R ( n − p ) } \{x|Ax=b\}=\{Fz+\hat{x}|z\isin R^{(n-p)}\} {xAx=b}={Fz+x^zR(np)}

        于 是 得 到 A ( F z + x ^ ) = b    ⟹    A F z = 0    ⟹    A p × n F n × ( n − p ) = 0 于是得到A(Fz+\hat{x})=b\implies AFz=0\implies A_{p\times n}F_{n\times(n-p)}=0 A(Fz+x^)=bAFz=0Ap×nFn×(np)=0

        于是A的行组成空间与F的列组成空间正交,即A的零空间就是F的值空间。

        min ⁡ f ( x ) , s . t . A x = b    ⟺    min ⁡ z F ^ ( z ) = f ( F z + x ^ ) \min{f(x)},s.t. Ax=b \iff \min_z \hat{F}(z)=f(Fz+\hat{x}) minf(x),s.t.Ax=bminzF^(z)=f(Fz+x^)【把等式约束加到目标函数里,变成无约束问题】

      • 求解对偶最优解 ν ∗ \nu^* ν

        原问题的KKT系统: ∇ f ( x ∗ ) + A T ν ∗ = 0 \nabla f(x^*)+A^T\nu^*=0 f(x)+ATν=0,同时左乘矩阵A(实现满秩运算)

        得到 ν ∗ = − ( A A T ) − 1 A ∇ f ( x ∗ ) \nu^*=-(AA^T)^{-1}A\nabla f(x^*) ν=(AAT)1Af(x)

        验证目标: ∇ f ( x ∗ ) − A T ( A A T ) − 1 A ∇ f ( x ∗ ) = 0 \nabla f(x^*)-A^T(AA^T)^{-1}A\nabla f(x^*)=0 f(x)AT(AAT)1Af(x)=0

        左乘一个满秩矩阵
        在这里插入图片描述
        基于KKT系统,存在u使得 ∇ f ( x ∗ ) + A T u = 0 \nabla f(x^*)+A^T u=0 f(x)+ATu=0,同时左乘 F T F^T FT,再因为 A F = 0 AF=0 AF=0
        于是推导得 F T ∇ f ( x ∗ ) = 0 F^T\nabla f(x^*)=0 FTf(x)=0;左乘矩阵A可以发现等式恒成立,综上,得证。

      • 原问题转变成
        在这里插入图片描述

    • 通过对偶问题实现等式消除

      L ( x , ν ) = f ( x ) + ν T ( A x − b ) L(x,\nu)=f(x)+\nu^T(Ax-b) L(x,ν)=f(x)+νT(Axb)

      g ( ν ) = inf ⁡ x L ( x , ν ) = − b T ν + inf ⁡ x ( f ( x ) + ν T A x ) g(\nu)=\inf_xL(x,\nu)=-b^T\nu+\inf_x(f(x)+\nu^TAx) g(ν)=infxL(x,ν)=bTν+infx(f(x)+νTAx)

      = − b T ν − sup ⁡ x ( − ( A T ν ) T x − f ( x ) ) = − b T ν − f ∗ ( − A T ν ) =-b^T\nu-\sup_x(-(A^T\nu)^Tx-f(x))=-b^T\nu-f^*(-A^T\nu) =bTνsupx((ATν)Txf(x))=bTνf(ATν)【共轭函数-第三章凸函数】

      对偶问题: max ⁡ ( − b T ν − f ∗ ( − A T ν ) ) , g ( ν ∗ ) = p ∗ \max{(-b^T\nu-f^*(-A^T\nu))},g(\nu^*)=p^* max(bTνf(ATν)),g(ν)=p

    • 例子

      • 变量替换
        在这里插入图片描述
        重点是找到特解 x ^ \hat x x^
        在这里插入图片描述
      • 对偶问题
        在这里插入图片描述
        在这里插入图片描述

    基于feasible初值的牛顿法

    • 牛顿法:利用二次导数逼近【设 △ x n t = v \triangle x_{nt}=v xnt=v

      问题变成: min ⁡ v f ^ ( x + v ) = f ( x ) + ∇ f ( x ) v + 1 2 v T ∇ 2 f ( x ) v , s . t . A ( x + v ) = b \min_v\hat f(x+v)=f(x)+\nabla f(x)v+\frac12v^T\nabla^2f(x)v, \quad s.t. A(x+v)=b minvf^(x+v)=f(x)+f(x)v+21vT2f(x)v,s.t.A(x+v)=b

      设x是一个可行解( A x = b , A v = 0 Ax=b,Av=0 Ax=b,Av=0),基于KKT系统,存在w使得 ∇ v f ^ ( x + v ) + A T w = 0 \nabla_v \hat f(x+v)+A^Tw=0 vf^(x+v)+ATw=0

      ∇ f ^ ( x + v ) = 0 \nabla \hat f(x+v)=0 f^(x+v)=0带入KKT条件得: { ∇ f ( x ) + ∇ 2 f ( x ) v + A T w = 0 A v = 0 \begin{cases}\nabla f(x)+\nabla^2f(x)v+A^Tw=0 \\ Av=0 \end{cases} {f(x)+2f(x)v+ATw=0Av=0
      在这里插入图片描述
      (演变成求解该矩阵方程,求 △ x n t = v , 对 偶 w \triangle x_{nt}=v,对偶w xnt=vw
      当A=0时,该等式问题就退化到无约束问题

    • Newton decrement

      • 性质
        在这里插入图片描述

      • 部分证明

        λ ( x ) = ( △ x n t T ∇ 2 f ( x ) △ x n t ) 1 2 \lambda(x)=(\triangle x_{nt}^T\nabla^2f(x)\triangle x_{nt})^{\frac12} λ(x)=(xntT2f(x)xnt)21,衡量 x 到 x ∗ x到x^* xx的距离,可以作为二阶逼近下的停止条件。

        f ^ ( x + △ x ) = f ( x ) + ∇ f ( x ) T △ x + 1 2 △ x T ∇ 2 f ( x ) △ x \hat f(x+\triangle x)=f(x)+\nabla f(x)^T\triangle x+\frac12\triangle x^T\nabla^2f(x)\triangle x f^(x+x)=f(x)+f(x)Tx+21xT2f(x)x

        H = ∇ 2 f ( x ) , g = ∇ f ( x ) , △ x = △ x n t = v H=\nabla^2f(x),g=\nabla f(x),\triangle x=\triangle x_{nt}=v H=2f(x)g=f(x)x=xnt=v, KKT系统: { H △ x + A T w = − g A △ x = 0 \begin{cases}H\triangle x+A^Tw=-g \\ A\triangle x=0 \end{cases} {Hx+ATw=gAx=0

        第一个式子左右两边左乘 △ x T \triangle x^T xT,得到 △ x T H △ x = − △ x T g \triangle x^TH\triangle x=-\triangle x^Tg xTHx=xTg

        带回得 f ^ ( x + △ x ) = f ( x ) − 1 2 △ x T ∇ 2 f ( x ) △ x \hat f(x+\triangle x)=f(x)-\frac12\triangle x^T\nabla^2f(x)\triangle x f^(x+x)=f(x)21xT2f(x)x

        f ( x ) − inf ⁡ v { f ^ ( x + v ) ∣ A ( x + v ) = b } = f ( x ) − f ^ ( x + △ x ) = 1 2 λ 2 ( x ) f(x)-\inf_v\{\hat f(x+v)|A(x+v)=b\}=f(x)-\hat f(x+\triangle x)=\frac12\lambda^2(x) f(x)infv{f^(x+v)A(x+v)=b}=f(x)f^(x+x)=21λ2(x)

      • 结论

        • 误差 f ( x ) − p ∗ f(x)-p^* f(x)p不断变小
        • 方向导数恒小于0
    • 算法步骤
      在这里插入图片描述

    • 牛顿法=消除等式+牛顿法
      【把等式约束加入(x可行域内)一起做牛顿下降法】,等价于【把变量消除掉,再去做牛顿法】

      结论: min ⁡ x f ( x ) , s . t . A x = b    ⟺    min ⁡ z f ~ ( z ) = f ( F z + x ^ ) , s . t . A x ^ = b \min_x{f(x)},s.t. Ax=b\iff \min_z\tilde f(z)=f(Fz+\hat x),s.t. A\hat x=b minxf(x),s.t.Ax=bminzf~(z)=f(Fz+x^),s.t.Ax^=b

      • 两种问题对比
        在这里插入图片描述

      • 证明

        ∇ f ~ ( z ) = F T ∇ f ( F z + x ^ ) , ∇ 2 f ~ ( z ) = F T ∇ 2 f ( F z + x ^ ) F \nabla \tilde f(z)=F^T\nabla f(Fz+\hat x),\nabla^2\tilde f(z)=F^T\nabla^2 f(Fz+\hat x)F f~(z)=FTf(Fz+x^)2f~(z)=FT2f(Fz+x^)F【加入等式约束】

        牛顿迭代步: △ z n t = − ∇ 2 f ~ ( z ) − 1 ∇ f ~ ( z ) = − ( F T ∇ 2 f ( F z + x ^ ) F ) − 1 F T ∇ f ( F z + x ^ ) \triangle z_{nt}=-\nabla^2 \tilde f(z)^{-1}\nabla \tilde f(z)=-(F^T\nabla^2 f(Fz+\hat x)F)^{-1}F^T\nabla f(Fz+\hat x) znt=2f~(z)1f~(z)=(FT2f(Fz+x^)F)1FTf(Fz+x^)

        X = F z + x ^    ⟹    △ x n t = F △ z n t X=Fz+\hat x\implies\triangle x_{nt}=F\triangle z_{nt} X=Fz+x^xnt=Fznt

        一阶KKT条件: ∇ f ( x + △ x n t ) + A T w = 0 \nabla f(x+\triangle x_{nt})+A^Tw=0 f(x+xnt)+ATw=0 ,得到 w = − ( A A T ) − 1 A ( ∇ f ( x ) + ∇ 2 f ( x ) △ x n t ) w=-(AA^T)^{-1}A(\nabla f(x)+\nabla^2f(x)\triangle x_{nt}) w=(AAT)1A(f(x)+2f(x)xnt)

        二阶KKT条件: ∇ 2 f ( x ) △ x n t + A T w + ∇ f ( x ) = 0 \nabla^2 f(x)\triangle x_{nt}+A^Tw+\nabla f(x)=0 2f(x)xnt+ATw+f(x)=0 ,得到 [ F T A ] [ ∇ 2 f ( x ) △ x n t + A T w + ∇ f ( x ) ] = 0 \begin{bmatrix}F^T\\A\end{bmatrix}\begin{bmatrix} \nabla^2 f(x)\triangle x_{nt}+A^Tw+\nabla f(x) \end{bmatrix}=0 [FTA][2f(x)xnt+ATw+f(x)]=0

        根据 F T A T w = 0 F^TA^Tw=0 FTATw=0,发现矩阵方程成立,于是得到 λ ~ 2 ( z ) = λ 2 ( x ) \tilde \lambda^2(z)=\lambda^2(x) λ~2(z)=λ2(x),即迭代逼近的效果相同

    基于infeasible初值的牛顿法

    初始点infeasible,但是希望 x + △ x n t x+\triangle x_{nt} x+xnt满足二阶KKT条件 { ∇ 2 f ( x ) △ x n t + A T w + ∇ f ( x ) = 0 A ( x + △ x n t ) = b \begin{cases}\nabla^2 f(x)\triangle x_{nt}+A^Tw+\nabla f(x)=0\\A(x+\triangle x_{nt})=b\end{cases} {2f(x)xnt+ATw+f(x)=0A(x+xnt)=b

    通过解原对偶问题解释原理

    • 设置迭代步 r ( x , v ) r(x,v) r(x,v)满足KKT条件,希望迭代步 r ( x , v ) r(x,v) r(x,v)不断逼近0
      在这里插入图片描述
      注:雅可比矩阵
      在这里插入图片描述
    • 算法步骤
      在这里插入图片描述
    • 下降性质
      infeasible初值的不一定一直保持下降,虽然函数不一定下降,但是残差一定下降
      在这里插入图片描述在这里插入图片描述

    求解KKT系统

    • KKT系统
      在这里插入图片描述
    • 求解方法
      • L D L T LDL^T LDLT分解
        B = P L D L T P T B=PLDL^TP^T B=PLDLTPT,P是置换矩阵,L是下三角矩阵,D是块矩阵
        在这里插入图片描述
        复杂度: O ( 1 3 ( n + p ) 3 ) O(\frac13(n+p)^3) O(31(n+p)3)

      • 消元法
        在这里插入图片描述
        在这里插入图片描述

        结果: v = − H − 1 ( g + A T w ) , w = ( A H − 1 A T ) − 1 ( h − A H − 1 g ) v=-H^{-1}(g+A^Tw),w=(AH^{-1}A^T)^{-1}(h-AH^{-1}g) v=H1(g+ATw),w=(AH1AT)1(hAH1g)

        时间复杂度 O ( n 3 + p 3 ) O(n^3+p^3) O(n3+p3)

    实例

    • 例子&如何降低计算
      在这里插入图片描述

      • 等式约束牛顿法
        在这里插入图片描述在这里插入图片描述

      • 对偶问题+无约束牛顿法
        在这里插入图片描述在这里插入图片描述

      • infeasible初值的牛顿法
        在这里插入图片描述在这里插入图片描述

      结论:无论哪种方法,(设D是正对角阵)总要计算 A D A T w = h ADA^Tw=h ADATw=h

    • Network flow optimization
      在这里插入图片描述在这里插入图片描述

    • Analytic center of linear matrix inequality
      在这里插入图片描述在这里插入图片描述

    总结

    • 基本定义: min ⁡ f ( x ) , s . t . A x = b \min{f(x)},s.t. Ax=b minf(x),s.t.Ax=b

    • 直接解KKT条件(例子-二次规划函数)
      { A x ∗ = b ∇ f ( x ∗ ) + A T ν ∗ = 0 \begin{cases}Ax^*=b \\ \nabla f(x^*)+A^T\nu^*=0 \end{cases} {Ax=bf(x)+ATν=0

    • 消除等式约束

      • 通过变量替换
        • x = F z + x ^ x=Fz+\hat x x=Fz+x^【通解+特解】
        • 求解对偶变量 ν ∗ = − ( A A T ) − 1 A ∇ f ( x ∗ ) \nu^*=-(AA^T)^{-1}A\nabla f(x^*) ν=(AAT)1Af(x)【LSE逼近】
      • 通过对偶问题
        max ⁡ ( − b T ν − f ∗ ( − A T ν ) ) , g ( ν ∗ ) = p ∗ \max{(-b^T\nu-f^*(-A^T\nu))},g(\nu^*)=p^* max(bTνf(ATν)),g(ν)=p
    • 基于feasible初值的牛顿法

      • 二阶展开,先做二阶展开,带入KKT求导得到极值,再套上迭代量含义。(一阶展开,将迭代量带入KKT条件,再做一阶展开)
      • Newton decrement: λ ( x ) = ( △ x n t T ∇ 2 f ( x ) △ x n t ) 1 2 \lambda(x)=(\triangle x_{nt}^T\nabla^2f(x)\triangle x_{nt})^{\frac12} λ(x)=(xntT2f(x)xnt)21
        可以衡量 x 到 x ∗ x到x^* xx的距离,可以作为二阶逼近下的停止条件
        • 误差 f ( x ) − p ∗ f(x)-p^* f(x)p不断变小
        • 方向导数恒小于0
      • 牛顿法=消除等式+牛顿法
        min ⁡ x f ( x ) , s . t . A x = b    ⟺    min ⁡ z f ~ ( z ) = f ( F z + x ^ ) , s . t . A x ^ = b \min_x{f(x)},s.t. Ax=b\iff \min_z\tilde f(z)=f(Fz+\hat x),s.t. A\hat x=b minxf(x),s.t.Ax=bminzf~(z)=f(Fz+x^),s.t.Ax^=b
    • 基于infeasible初值的牛顿法

      • 基本定义:

        初始点infeasible,但是希望 x + △ x n t x+\triangle x_{nt} x+xnt满足二阶KKT条件 { ∇ 2 f ( x ) △ x n t + A T w + ∇ f ( x ) = 0 A ( x + △ x n t ) = b \begin{cases}\nabla^2 f(x)\triangle x_{nt}+A^Tw+\nabla f(x)=0\\A(x+\triangle x_{nt})=b\end{cases} {2f(x)xnt+ATw+f(x)=0A(x+xnt)=b

      • 通过解原对偶问题解释原理

        • 设置迭代步 r ( x , v ) r(x,v) r(x,v)满足KKT条件
        • 希望迭代步 r ( x , v ) r(x,v) r(x,v)不断逼近0
      • 下降性质
        infeasible初值的不一定一直保持下降,虽然函数不一定下降,但是残差一定下降

      • 具体的算法流程
        在这里插入图片描述

    • KKT系统方程的求解

      • L D L T LDL^T LDLT分解

        普通矩阵 B = P L D L T P T B=PLDL^TP^T B=PLDLTPT,P是置换矩阵,L是下三角矩阵,D是块矩阵

        复杂度: O ( 1 3 ( n + p ) 3 ) O(\frac13(n+p)^3) O(31(n+p)3)

      • 消元法

        v = − H − 1 ( g + A T w ) , w = ( A H − 1 A T ) − 1 ( h − A H − 1 g ) v=-H^{-1}(g+A^Tw),w=(AH^{-1}A^T)^{-1}(h-AH^{-1}g) v=H1(g+ATw),w=(AH1AT)1(hAH1g)

        时间复杂度 O ( n 3 + p 3 ) ∼ O ( n 3 ) O(n^3+p^3)\sim O(n^3) O(n3+p3)O(n3)

    • 三种求解牛顿法的实例

      • 等式约束牛顿法
      • 对偶问题+无约束牛顿法
      • infeasible初值的牛顿法
    展开全文
  • 等式约束优化

    千次阅读 2018-07-04 16:10:21
    介绍等式约束优化的求解。 等式约束优化问题 minf(x)minf(x)\min f(x) s.t.Ax=bs.t.Ax=bs.t.\quad Ax=b 其中fff为二次可微凸函数,假设等式约束少于变量数,并且等式约束互相独立。假定存在一个最优解x⋆x...

    介绍等式约束优化的求解。

    等式约束优化问题

    minf(x) min f ( x )

    s.t.Ax=b s . t . A x = b

    其中 f f 为二次可微凸函数,假设等式约束少于变量数,并且等式约束互相独立。假定存在一个最优解x,并用 p p ⋆ 表示其最优值,即:

    p=inf{f(x)|Ax=b}=f(x) p ⋆ = inf { f ( x ) | A x = b } = f ( x ⋆ )

    由KKT条件,其最优解的重要条件是满足:

    Ax=bf(x)+ATv=0 A x ⋆ = b ▽ f ( x ⋆ ) + A T v ⋆ = 0

    对于求解等式约束问题有两种方法:

    1. 任何等式约束优化问题都可以通过消除等式约束转化为等价的无约束问题。
    2. 使用对偶方法解决。

    很多时候,直接处理等式约束比转化为无约束问题要好,这是因为转化之后可能会破坏问题的结构。

    等式约束凸二次规划

    minf(x)=(1/2)xTPx+qTx+r min f ( x ) = ( 1 / 2 ) x T P x + q T x + r

    s.t.Ax=b s . t . A x = b

    此问题的最优性条件为:

    Ax=bPx+q+ATv=0 A x ⋆ = b P x ⋆ + q + A T v ⋆ = 0

    可以将其写成矩阵形式:

    [PAAT0][xv]=[qb] [ P A T A 0 ] [ x ⋆ v ⋆ ] = [ − q b ]

    这个矩阵称为KKT矩阵,接下来会经常用到。

    消除等式约束

    我们以参数化可行集的形式表示等式约束:

    {x|Ax=b}={Fz+x^} { x | A x = b } = { F z + x ^ }

    其中 x^ x ^ 为任意特解, F F A的零空间的任意矩阵,可以消除等式约束为:

    minf^(z)=f(Fz+x^) min f ^ ( z ) = f ( F z + x ^ )

    这里的变量 z z 没有约束,利用它的解z可以确定等式约束问题的解 x=Fz+x^ x ⋆ = F z ⋆ + x ^

    对偶方法求解等式约束

    可得约束问题的对偶函数为:

    g(v)=bTv+infx(f(x)+vTAx)=bTxsupx((ATv)Txf(x))=bTvf(ATv)(21)(22)(23) (21) g ( v ) = − b T v + inf x ( f ( x ) + v T A x ) (22) = − b T x − sup x ( ( − A T v ) T x − f ( x ) ) (23) = − b T v − f ⋆ ( − A T v )

    因此,对偶问题为:

    maxbTvf(ATv) max − b T v − f ⋆ ( − A T v )

    Slater条件成立,则强对偶性成立,即 g(v)=p g ( v ⋆ ) = p ⋆

    等式约束的Newton方法

    讨论扩展的Newton方法,与之前无约束类似,但初始点必须可行(即满足 Ax=b A x = b ),并且需要保证Newton方向是可行的方向,即 AΔx