復習シリーズ2. 表題の内容について自分の理解でまとめます.

前回:

制約つき最適化問題
等式制約問題
不等式制約問題
実装

制約つき最適化問題

とは, 以下のように定式化されるものです.

$\mathcal{X} \subset \mathbb{R}^d$ 上の微分可能な目的関数 $f:\mathcal{X}\mapsto\mathbb{R}$ の最小値を求めよ. $$ \min_{\boldsymbol{x}\in\mathcal{X}} f(\boldsymbol{x}) \ \ \ \mathrm{subject \ to} \ \ \ \vec{g}(\boldsymbol{x})=\boldsymbol{0}, \ \vec{h}(\boldsymbol{x})\leq \boldsymbol{0} $$ $$ \vec{g}(\boldsymbol{x}) = ( g_1(\boldsymbol{x}),\ldots,g_m(\boldsymbol{x}))^{\mathrm{T}}$$ $$ \vec{h}(\boldsymbol{x}) = ( h_1(\boldsymbol{x}),\ldots,h_n(\boldsymbol{x}))^{\mathrm{T}}$$

ベクトルについては, 変数・定数は太字表記・関数は矢印表記ということにします. (あまり一般的ではないと思いますが, 自分が学習中に関数のベクトルに混乱したことがあったので…）

"subject to…" は「…という制約条件を持つ」という意味で, ベクトルの等式・不等式は要素ごとであるものとします.

この定式化では, $\boldsymbol{x}$ について, $ m $ 個の等式と $n$ 個の不等式による制約条件がかかっているという意味になります. （$m=0$ または $n=0$ でもかまいませんが, $m=n=0$ だと制約なしになります）

ある $\boldsymbol{x}_0$ が全ての制約条件 $\vec{g}(\boldsymbol{x}_0)=\boldsymbol{0}, \ \vec{h}(\boldsymbol{x}_0)\leq \boldsymbol{0} $ を満たすとき, $\boldsymbol{x}_0$ を実行可能解 feasible solution と呼びます. 実行可能解の中から目的関数を最小にするものを探そうというわけです.

結論から言うと, 制約つき最適化問題は, 等式・不等式問わず, 制約なし最適化問題に変換することで解くことができます.

そのため, 以下では等式・不等式制約問題の性質について触れつつ, そのいずれにも利用できる 二乗罰則法 を紹介します *1 . また, 不等式制約問題については, 対数バリア法 も紹介します.

その後, 二乗罰則法および対数バリア法を Julia で実装し, 簡単な不等式制約つき最適化問題を解く例を示します.

等式制約問題

まずは等式制約つき最適化問題から.

$$ \min_{\boldsymbol{x}\in\mathcal{X}} f(\boldsymbol{x}) \ \ \ \mathrm{subject \ to} \ \ \ \vec{g}(\boldsymbol{x})=\boldsymbol{0}$$ $$ \vec{g}(\boldsymbol{x}) = ( g_1(\boldsymbol{x}),\ldots,g_m(\boldsymbol{x}))^{\mathrm{T}}$$

前述の通り, 二乗罰則法という手法が知られています. そのアルゴリズムは以下のようなものです.

各 $k$ について, 次の $\boldsymbol{x}_{k+1} $を,

$\boldsymbol{x}_{k}$を初期値とした以下の制約なし最適化問題の解として選ぶ: $$ \boldsymbol{x}_{k+1} \leftarrow \underset{\boldsymbol{x} \in \mathcal{X}}{\mathrm{argmin}} \
\left\{ f(\boldsymbol{x}) + c_k ||\vec{g}(\boldsymbol{x})||^2 \right\} $$ $c_k$ は $k$ の増加に従って徐々に増やす: $$ 0 < c_1 < c_2 < \ldots $$

これは, $ \vec{g}(\boldsymbol{x})=\boldsymbol{0}$ という制約条件を, 「$ \vec{g}(\boldsymbol{x})=\boldsymbol{0}$ からできるだけ離れないようにする」と読み替えたものと言えるでしょう.

すなわち, 点 $\boldsymbol{x}$ と直線（あるいは超平面） $ \vec{g}(\boldsymbol{x})=\boldsymbol{0}$ からの距離に応じて増える罰則項 $ c_k ||\vec{g}(\boldsymbol{x})||^2$ を加味した $$ f(\boldsymbol{x}) + c_k ||\vec{g}(\boldsymbol{x})||^2 $$ を最小化するという制約なし最適化問題を解くということです.

ここで, 「$c_k$ は $k$ の増加に従って徐々に増やす」ことで, 始めはほぼ$f(\boldsymbol{x})$単体での最小化のような挙動をしながら, 徐々に制約条件$ \vec{g}(\boldsymbol{x})=\boldsymbol{0}$ を満たすように$\boldsymbol{x}_{k} $ が移動していくことを狙うアルゴリズムだといえます.

$c_k$をいかに決め, いかに増加させるかが難しい点です.

不等式制約問題

二乗罰則法は, 不等式制約問題:

$$ \min_{\boldsymbol{x}\in\mathcal{X}} f(\boldsymbol{x}) \ \ \ \mathrm{subject \ to} \ \ \ \vec{h}(\boldsymbol{x})\leq\boldsymbol{0}$$ $$ \vec{h}(\boldsymbol{x}) = ( h_1(\boldsymbol{x}),\ldots,h_n(\boldsymbol{x}))^{\mathrm{T}}$$

でも適用できます. 等式制約問題と同じ発想で,

各 $k$ について, 次の $\boldsymbol{x}_{k+1} $を,

$\boldsymbol{x}_{k}$を初期値とした以下の制約なし最適化問題の解として選ぶ: $$ \boldsymbol{x}_{k+1} \leftarrow \underset{\boldsymbol{x} \in \mathcal{X}}{\mathrm{argmin}} \
\left\{ f(\boldsymbol{x}) + c_k \max\left(0, \vec{h}(\boldsymbol{x})\right)^2 \right\} $$ $c_k$ は $k$ の増加に従って徐々に増やす: $$ 0 < c_1 < c_2 < \ldots $$ ただし, $\max\left(0, \vec{h}(\boldsymbol{x})\right)^2 = \displaystyle \displaystyle \sum_{i=1}^{n}\max\left(0, h_i (\boldsymbol{x})\right)^2 $ とする.

という形になります.

等式制約問題では, $ \vec{g}(\boldsymbol{x})<\boldsymbol{0}$ と $ \vec{g}(\boldsymbol{x})>\boldsymbol{0}$ の両方にペナルティを課す必要がありましたが, 今回は $ \vec{h}(\boldsymbol{x})\leq\boldsymbol{0}$ という条件であるため, $\vec{h}(\boldsymbol{x})<\boldsymbol{0}$ である分には何も問題がないので, $ \vec{g}(\boldsymbol{x})>\boldsymbol{0}$ 側にだけ罰則を課すために, 罰則項が $ \max \left( 0, \vec{h}(\boldsymbol{x}) \right)$ となるという違いが生じています.

このアルゴリズムでは, 途中の解 $\boldsymbol{x}_k$は一般に, 実行可能領域の外側にも来る（$ \vec{g}(\boldsymbol{x}_k)>\boldsymbol{0}$）ため, 外点法とも呼ばれます.

等式制約を厳密に満たすことはそもそも数値計算では困難 *2 なので, 等式制約問題では解が実行可能領域から少し外れていてもあまり問題にはならないと考えられますが, 不等式制約問題では, 「実行可能領域の内側に解があってくれないと困る」という要求はありえます.

そんなとき, 前述の外点法に対して 内点法 と呼ばれる, 「常に実行可能領域の内側に解が存在することを保証する」手法のひとつが 対数バリア法 で, 以下のようなアルゴリズムです.

各 $k$ について, 次の $\boldsymbol{x}_{k+1} $を,

$\boldsymbol{x}_{k}$を初期値とした以下の制約なし最適化問題の解として選ぶ: $$ \boldsymbol{x}_{k+1} \leftarrow \underset{\boldsymbol{x} \in \mathcal{X}}{\mathrm{argmin}} \
\left\{ f(\boldsymbol{x}) - c_k \log\right(-\vec{h}(\boldsymbol{x})\left) \right\} $$ $c_k$ は $k$ の増加に従って徐々に減らす: $$ c_1 > c_2 > \ldots > 0 $$ ただし, $\log\left(-\vec{h}(\boldsymbol{x})\right) =\displaystyle \sum_{i=1}^{n}\log\left(-h_i(\boldsymbol{x})\right) $ とする.

二乗罰則法での罰則項が $- c_k \log\left(-\vec{h}(\boldsymbol{x})\right) $ に変わりました. 障壁項と呼ぶことにします.

実数の範囲では $\log(x)$ は $x\leq 0$ では定義されず, $\underset{x\to +0}{\lim} \log(x) \to -\infty$ であることは（ここまで読んできた方なら）周知かと思います.

これを踏まえると, 上の障壁項は,

「実行可能領域の内側から境界に近づけば $$\underset{\vec{h}(\boldsymbol{x}) \to -0}{\lim} \left\{- c_k \log\left(-\vec{h}(\boldsymbol{x})\right) \right\} \to \infty$$ となり, 実行可能領域の境界上 $\boldsymbol{x}=\boldsymbol{0}$ および外側 $\boldsymbol{x}>\boldsymbol{0}$ では定義されない」

ようなものであることがわかります.

ここで, $\log(x)$ は $x=1$ で正負が入れ替わるので, $-\vec{h}(\boldsymbol{x})\geq \boldsymbol{1}$ のときはむしろ目的関数を減少させる項となってしまい, おかしくなるのでは？と思われるかもしれません.

ですが, $c_k$ は $k$ の増加に従って $ c_1 > c_2 > \ldots > 0 $ と減らしていくものであったことを思うと, $c_k$が $0$ に近づくとき, 障壁項は

${ \underset{c_k \to +0}{\lim}\left\{- c_k \log\left(-\vec{h}(\boldsymbol{x})\right) \right\} \to \left\{ \begin{array}{l} \infty \ \ (\vec{h}(\boldsymbol{x}) =\boldsymbol{0}) \\ 0 \ \ (\vec{h}(\boldsymbol{x}) \lt \boldsymbol{0}) \end{array} \right. }$

となるので, 問題なさそうです.

実装

それでは, Julia の時間です.

using LinearAlgebra

今回, 例として以下の最適化問題を数値的に解きます.

$$ \min_{x,y \in \mathbb{R}} 3 x^2 + 2 y^2 \ \mathrm{subject \ to} \ \ x+y \geq 1 $$

解析解は $(x,y) = (0.4,0.6)$ です.

目的関数と制約条件の定義と, 今回はステップ幅一定の単純な勾配法を用いることにします *3 .

function f(x)
    3x[1]^2+2x[2]^2
end

function h(x)
    1-(x[1]+x[2])
end

function ∇(f,x)
    h = 0.0001
    return [( f(x+[h,0]) - f(x-[h,0]) ) / 2h, 
            ( f(x+[0,h]) - f(x-[0,h]) ) / 2h]
end


function GradientDescent(f, init_x, ε=0.05)
    x = init_x
    m = 0.1
    while norm( ε.*∇(f,x) ) > m
            x = x - ε.*∇(f,x)
    end
    return x
end

まずは二乗罰則法です.

function SquarePenalty(f,h)
    x = (rand(2).+1)*20
    c = 2    #initialize
    ε = 2  # cの増加率
    m = 0.001  # 収束判定
    
    i = 1
    result = []
    append!(result, [x])
    print(i, ":", x, " c:", c, "\n")
    
    while i==1 || norm( result[i-1] - result[i] ) > m
        i += 1
        c = ε*c
        SquarePenaltyProblem(x) = f(x) + c * max( 0, h(x) )^2
        x_next = GradientDescent(SquarePenaltyProblem, x, 0.05)
        append!(result, [x_next])
        print(i, ":", x_next, " c:", c, "\n")
        x = x_next
    end 
    
    return result
end

各ループで, $\mathrm{argmin} \left\{ f(\boldsymbol{x}) + c_k \max\left(0, \vec{h}(\boldsymbol{x})\right)^2 \right\}$ にあたる SquarePenaltyProblem() を定義し直してから GradientDescent() に与えています.

次に, 対数バリア法を同様に実装します.

function LogBarrier(f,h)
    x = (rand(2).+1)*20
    c = 10    #initialize
    ε = 0.5  # cの減衰率
    m = 0.01  # 収束判定
    
    i = 1
    result = []
    append!(result, [x])
    print(i, ":", x, " c:", c, "\n")
    
    while i==1 || norm( result[i-1] - result[i] ) > m
        i += 1
        c = ε*c
        logBarrierProblem(x) = f(x) - c * log( -h(x) )
        x_next = GradientDescent(logBarrierProblem, x)
        append!(result, [x_next])
        print(i, ":", x_next, " c:", c, "\n")
        x = x_next
    end 
    
    return result
end

それぞれ実行します.

SquarePenalty(f,h)

1:[33.5958, 39.9659] c:2
2:[0.0634654, 0.728719] c:4
3:[0.276879, 0.65636] c:8
4:[0.276879, 0.65636] c:16

LogBarrier(f,h)

1:[27.8908, 33.7221] c:10
2:[0.583505, 1.54771] c:5.0
3:[0.507339, 1.22301] c:2.5
4:[0.440713, 1.06398] c:1.25
5:[0.370418, 0.913105] c:0.625
6:[0.314403, 0.785594] c:0.3125
7:[0.314403, 0.785594] c:0.15625

解析解の $(x,y) = (0.4,0.6)$ とはややズレていますが, 概ねその近くに収束した……と言って良いのかなと思います.

*1:ラグランジュ双対問題を利用する解き方など他にも色々あるのですが, 長くなるので今回は割愛します

*2:「二次元平面で厳密に直線上に点を乗せること」などをイメージすると……

*3:バックトラック直線探索を行ったところ, 対数バリア法でバリアを突き抜けた位置に点が移動するなどしたため…