はじめに
数理最適化をある程度学ぶと，パラメータ（外生変数）を含むような最適化問題を扱わなければならないときがある．そういうとき，定数を書き直してまた最適化計算に掛けるという作業はどうしても面倒になる．特に近似解法と呼ばれる類いのアルゴリズム（本稿でもペナルティ法を例に紹介する）を実装するとき，必ず近似問題の条件の厳しさを調節するパラメータが存在するだろう．それを何十回と反復的に解く事を考えれば，for文で回せないかなぁ，なんて思うのも当然である．
本稿では，下記のパラメータ $a\in\Re^t$ を含む数理最適化問題
$$
\begin{align}
\begin{array}{lcl}
{\rm P}(a): & \underset{x\in\Re^n}{\rm minimize} & f(x;a)\\
& {\rm subject\ to} & g_i(x;a)\leq 0\ (i=1,...,l)\\
& & h_j(x;a)=0\ (j=1,...,m)
\end{array}
\end{align}
$$
を MATLAB で解くためのプログラム作成方法を述べる．ここで，関数 $f(\cdot;a):\Re^n\to\Re,$ $g(\cdot;a):\Re^n\to\Re^l,$ $h(\cdot,a):\Re^n\to\Re^m$ は $a\in\Re^t$ をパラメータとする微分可能な関数である．ここでは，問題 ${\rm P}(a)$ は任意の $a$ に対して凸最適化問題となると仮定する．すなわち，局所的最適解が大域的最適解となるような問題を取り扱う．
少し余談であるが，この記事を書こうと思ったきっかけは，MATLAB を最適化計算で利用する初心者がパラメータを含む最適化問題を解こうとしたときにどう書けば良いかわからないとよく著者に聞くためである．
対象とする読者，難易度
本記事は数理最適化に関する知識を要する．なるべく，多くの人に理解していただけるよう，専門用語や数式の多用は避けたいところであるが，内容の都合上，若干ハードルが上がってしまうことをご了承願いたい．
- Matlab を使って最適化計算をしようと考えている人
- 数理最適化に関してある程度興味がある人
- 数理計画法をかじったことのある人
今後ニーズがあれば，数理最適化初学者向けの記事も書いてみようと思っている．
実装環境
- OS: Windows 10 Home
- Matlab version: 2016b
- 使用した Matlab ツール：
- fmincon（制約付き非線形計画問題ソルバ）
- fminunc（無制約非線形計画問題ソルバ）
パラメータを含む最適化問題の例
”はじめに”ではあまりにも最適化問題を抽象的に書いてしまったので，どういう問題を扱おうとしているのかわかりにくいところがある．例えば，下記の非線形計画問題
$$
\begin{align}
\begin{array}{cl}
\underset{x\in\Re^3}{\rm minimize} & f(x):=3x_1^2+2x_2^2+5x_3^2+x_1-x_2-x_3\\
{\rm subject\ to} & 2x_1-3x_2+x_3 \leq a_{t,1}\\
& -5x_1+3x_2-x_3 \leq a_{t,2}\\
\end{array}
\end{align}
$$
を考える．いま，$t=1,2,3$ が与えられてて，$a_{t,j}\ (j=1,2)$ は時刻 $t$ に関して変化する何らかの定数としよう．この非線形計画問題を解くなら，$a_{t,j}$ に関する3×2の行列Aを定義して
A = [-3 1; -1 4; -2 5];
f = @(x)(3*x(1)^2+2*x(2)^2+5*x(3)^2+x(1)-x(2)-x(3));
x0 = [0; 0; 0];
for t=1:3
fmincon(f,x0,[2 -3 1; -5 3 -1],[A(t,1); A(t,2)])
end
のように入力して解けば，それぞれの $t$ に関して解を出してくれるだろう．しかし，以下のように目的関数にパラメータを含む場合はどう扱えば良いだろうか．
$$
\begin{align}
\begin{array}{cl}
\underset{x\in\Re^3}{\rm minimize} & f(x):=c_{t,1}x_1^2+c_{t,2}x_2^2c_{t,3}x_3^2+x_1-x_2-x_3\\
{\rm subject\ to} & 2x_1-3x_2+x_3 \leq -3\\
& -5x_1+3x_2-x_3 \leq 1\\
\end{array}
\end{align}
$$
試しに
c = [3 2 5; 1 5 2; 3 1 2];
f = @(x,t)(c(t,1)*x(1)^2+c(t,2)*x(2)^2+c(t,3)*x(3)^2+x(1)-x(2)-x(3));
x0 = [0; 0; 0];
for t=1:3
fmincon(f,x0,[2 -3 1; -5 3 -1],[-3; 1])
end
と入力しても，
Not enough input arguments.
Error in @(x,t)(c(t,1)*x(1)^2+c(t,2)*x(2)^2+c(t,3)*x(3)^2+x(1)-x(2)-x(3))
Error in fmincon (line 536)
initVals.f = feval(funfcn{3},X,varargin{:});
Caused by:
Failure in initial objective function evaluation. FMINCON cannot continue.
と言われてしまう．fmincon の行を
fmincon(@(x)f,x0,[2 -3 1; -5 3 -1],[-3; 1])
としても，
Error using fmincon (line 684)
FMINCON requires all values returned by functions to be of data type double.
と弾かれてしまう．つまり，変数の文字を特定するだけでは x を指しているのか t を指しているのかわからず，関数の引数情報を入力してあげる必要がある．したがって，
fmincon(@(x)f(x,t),x0,[2 -3 1; -5 3 -1],[-3; 1])
とすれば，
Local minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.
<stopping criteria details>
ans =
0.6667
1.4255
-0.0567
Local minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.
<stopping criteria details>
ans =
0.6667
1.2174
-0.6812
Local minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.
<stopping criteria details>
ans =
0.6667
1.4737
0.0877
>>
と $t=1,2,3$ それぞれの場合の解を求めることができる．これを応用することで，次節のような解法も Matlab で実装することができる．
応用例：ペナルティ法を用いた最適化計算プログラムの作成
本節では，標題の本質をわかりやすくするため，${\rm P}(a)$ のパラメータ $a\in\Re^t$ を取り除いた問題
$$
\begin{align}
\begin{array}{lcl}
{\rm P}: & \underset{x\in\Re^n}{\rm minimize} & f(x)\\
& {\rm subject\ to} & g_i(x)\leq 0\ (i=1,...,l)\\
& & h_j(x)=0\ (j=1,...,m)
\end{array}
\end{align}
$$
を考えることにする．ペナルティ法は制約条件にペナルティパラメータと呼ばれる係数を掛け合わせた関数を目的関数に導入し，制約条件を含まない無制約非線形計画問題として取り扱い，元の制約条件を違反した分だけ罰則を与えるようにペナルティパラメータを大きくしながら逐次的に解く古くから知られた手法である．すなわち，${\rm P}$ に対するペナルティ（近似）問題はペナルティパラメータ $p>0$ を用いて以下のように書ける．
$$
\begin{array}{ccl}
\tilde{\rm P}(p): & \underset{x\in\Re^n}{\rm minimize} & \displaystyle f(x)+p\bigg(\sum_{i=1}^l \max\{g_i(x),0\}^q + \sum_{j=1}^m |h_j(x)|^q\bigg)
\end{array}
$$
ここで，$q>0$ は指数である（後述するが，一般には $q=2$ が用いられる）．実装方法としては，ある $q$ に固定してペナルティパラメータの初期値 $p_0$ を与える（はじめは $0\leq p\leq 1$ が望ましい）．$\tilde{\rm P}(p_0)$ を解き，制約条件を十分に満たしているなら終了し，満たしてなければ $p_0 < p$ なる $p$ を代入して，再び解く．このとき，前のステップで得られた解をアルゴリズムの初期点に与えておくと比較的良い収束性が得られる．これを繰り返すことで，$p$ が大きくなって行くにつれ，ペナルティ項（ $p$ が掛かっている項）の重みが増すため，$p$ の掛かっている項は小さくなるような解を見つける．すなわち，$\max\{g_i(x),0\}$ も $|h_j(x)|$ も0に収束することになり，元の制約条件を満たすような解が得られるという手法である．
さて，指数の話になるが，$q=1$ としたときは正確なペナルティ法（exact penalty method: EPM）と呼ばれ，有限回数で制約条件を「厳密に」満たす解が得られることが知られている．しかし，1-ノルム関数は連続的微分可能とはならないため，微分を使うアルゴリズムを直接使うことはできない．
一方，$q=2$ としたときは二乗ペナルティ法と呼ばれ，$g_i$ もしくは $h_j$ が微分可能なら，2-ノルムをとった関数も微分可能であるため，$q=1$ のときよりも比較的容易に解くことができる．しかし，正確なペナルティ法とは違い，解析的には $p\to\infty$ とすると元の問題を遵守する解へ収束する12という性質を持っている．したがって，実装するときは有限回数で終了し，近似解を得るに留まることになる．本稿では $q=2$ の場合を取り扱う．
テスト問題
$$
\begin{align}
\begin{array}{cl}
\underset{(x_1,x_2)\in\Re^2}{\rm minimize} & \displaystyle f(x):=\frac{1}{2}(2x_1^2+2x_1x_2+3x_2^2)-x_1-2x_2\\
{\rm subject\ to} & 3x_1-x_2+1\leq 0\\
& -2x_1+x_2+2\leq 0
\end{array}
\end{align}
$$
この問題は明らかに凸性を満たしている（目的関数は2次関数でヘッセ行列が正定値になることと，制約関数が一次関数であることから）．この問題のペナルティ問題を記述すると以下のようになる．
$$
\begin{array}{cl}
\underset{(x_1,x_2)\in\Re^2}{\rm minimize} & \displaystyle \frac{1}{2}(2x_1^2+2x_1x_2+3x_2^2)-x_1-2x_2+p(\max\{3x_1-x_2+1,0\}^2+\max\{-2x_1+x_2+2,0\}^2)
\end{array}
$$
この問題は $p>0$ をパラメータとし，$(x_1,x_2)$ を変数とする最適化問題である．また，この問題の最適解は凸性から最適性の必要条件（Karush-Kuhn-Tucker 条件）を用いて解析的に解くことができ，最適解は $(x_1,x_2)=(-3,-8)$ である．
プログラムの説明
便宜上，目的関数，制約関数はそれぞれ以下の式で簡略的に表している．
$$
\displaystyle f(x):=\frac{1}{2}x^\top Ax+b^\top x,\ g(x):=Dx+d\ (\leq 0)
$$
$$
A:=\left[\begin{array}{cc}
2 & 1\\
1 & 3
\end{array}\right],
b:=\left[\begin{array}{c}
-1\\
-2
\end{array}\right],
D:=\left[\begin{array}{cc}
3 & -1\\
-2 & 1
\end{array}\right],
d:=\left[\begin{array}{c}
1\\
2
\end{array}\right]
$$
アルゴリズムの終了条件は以下の式で与えている．
$$
\max\{3x_1-x_2+1,-2x_1+x_2+2\}<\varepsilon
$$
$\varepsilon>0$ は0に十分に近い数で，下のサンプルプログラムでは $\varepsilon:=1.0\times10^{-4}$ と与えている．
% 関数の係数を行列，ベクトルを用いて表現
A = [2 1; 1 3];
b = [-1; -2];
D = [3 -1; -2 1];
d = [1; 2];
% 目的関数
f = @(x) (1/2*x'*A*x+b'*x);
% 制約関数
g_1 = @(x) (D(1,:)*x+d(1));
g_2 = @(x) (D(2,:)*x+d(2));
% ペナルティ関数
P = @(x,p) (f(x)+p*(max(0,g_1(x))^2+max(0,g_2(x))^2));
% 逐次計算の最大反復数
Kmax = 20;
% ペナルティパラメータの初期値
p = 1;
% アルゴリズムの初期点を決める
x0 = [0;0];
xk = x0;
% 解の終了条件（許容誤差）を決める
eps = 1.0e-04;
% 計算
for k=1:Kmax
[xk,fvalk] = fminunc(@(x)P(x,p),xk);
% xk: 第kステップ目の解
% fvalk: 第kステップ目の目的関数値
% 終了条件（制約条件を十分に満たすかどうか）
if (max(g_1(xk),g_2(xk))<eps)
% 制約条件を十分に満たす解が見つかった場合の出力
fprintf('Solution found: [%.3f,%.3f]\n',xk(1),xk(2));
break
elseif (k<Kmax)
p = p * 10; % ペナルティパラメータの更新
elseif (k==Kmax)
% 指定の回数(Kmax回)以内で解けなかった場合の出力
fprintf('No solution found.\n');
end
end
>> example
Warning: Gradient must be provided for trust-region algorithm; using quasi-newton
algorithm instead.
> In fminunc (line 395)
In execution (line 31)
Local minimum found.
Optimization completed because the size of the gradient is less than
the default value of the optimality tolerance.
<stopping criteria details>
Warning: Gradient must be provided for trust-region algorithm; using quasi-newton
algorithm instead.
（中略...）
> In fminunc (line 395)
In execution (line 31)
Local minimum found.
Optimization completed because the size of the gradient is less than
the default value of the optimality tolerance.
<stopping criteria details>
Solution found: [-3.000,-8.000]
>>
上のように大域的最適解が求まった．ちなみに，ユーザーが目的関数の勾配を入力すれば Warning メッセージは消える上，計算速度も圧倒的に速くなる．
ペナルティ法を利用する上で注意したいのが，ペナルティパラメータを大きくしすぎると問題の性質が悪条件（ill-conditioned）になるため，解の精度はある程度の限界があるということである．上の例でも最終的に $p=1000000$ で終了したが，決して最適化問題として良い性質を持っているとは言い難い．それは元の問題とペナルティ問題でちょうど制約条件の境目となるところで，少しでも制約条件を違反すると飛躍的に数値が上がることから，制約条件の境目付近で折れ曲がったような概形をとるため，微分した関数を陽に与えない近似微分で解く場合，微分情報が著しく劣る．
実装のポイントとしては，関数ハンドルを如何にして使うかというところだと思う．
まとめ
実装ポイントをまとめると以下のようになる
> f = @(決定変数,外生変数) (目的関数)
で目的関数を記述し，ソルバで求解するときは目的関数のあとに引数を入力しておき@ハンドルを利用して決定変数を特定してあげるということである．
> fmincon(@(決定変数)f(決定変数,外生変数),...
と書くことで，決定変数をfmincon (fminunc等ツール)に認識させることができる．