二次計画法とは
最大化もしくは最小化したい目的関数が二次関数であり、制約条件が一次関数である最適化問題のことを二次計画問題と呼び、その解法のことを二次計画法といいます。
具体的には
目的関数 $5 x_1 ^2+10x_2^2+50x_1+25x_2$ （最小化）
制約条件 $x_1≧1, x_2≧2$
のような問題です。
一般的には $x=(x1, x2, …)$ として、二次計画問題を
目的関数 $x ^ {T} Q x + c ^ {T} x $
制約条件 $Ax ≧ b$
のようにベクトルで表記します。
上の問題だと
目的関数 $ \begin {pmatrix} x_1 & x_2 \\ \end {pmatrix} \begin {bmatrix} 5 & 0 \\ 0 & 10 \\ \end {bmatrix} \begin {pmatrix} x_1 \\ x_2 \end {pmatrix} + \begin {pmatrix} 50 & 25 \end {pmatrix} \begin {pmatrix} x_1 \\ x_2 \end {pmatrix} $
制約条件 $\begin {bmatrix} 1 & 0 \\ 0 & 1 \\ \end {bmatrix} \begin {pmatrix} x_1 \\ x_2 \\ \end {pmatrix} ≧ \begin {pmatrix} 1 \\ 2 \\ \end {pmatrix} $
となります。
ライブラリの使い方
上記で説明した二次計画問題をRを使って解く方法をご紹介します。
今回ご紹介するのは、パッケージ：quadprogのsolve.QP()関数です。
quadprogのインストール
solve.QPを使用するためにパッケージquadprogをインストールします。
Rのコンソール画面で下記実行します。
install.packages("quadprog")Rstudioユーザーの場合は
Packes → Install からquadprogと入力し、インストールすると早いです。
solve.QP()の使い方
目的関数 $ \begin {pmatrix} x_1 & x_2 \\ \end {pmatrix} \begin {bmatrix} 5 & 0 \\ 0 & 10 \\ \end {bmatrix} \begin {pmatrix} x_1 \\ x_2 \end {pmatrix} + \begin {pmatrix} 50 & 25 \end {pmatrix} \begin {pmatrix} x_1 \\ x_2 \end {pmatrix} $
制約条件 $\begin {bmatrix} 1 & 0 \\ 0 & 1 \\ \end {bmatrix} \begin {pmatrix} x_1 \\ x_2 \\ \end {pmatrix} ≧ \begin {pmatrix} 1 \\ 2 \\ \end {pmatrix} $
を最小化する場合のソースコードを掲載します。
#ライブラリの読み込み
library(quadprog)
#行列の定義
Dmat <- matrix(c(5,0,0,10),2,2) #目的関数の行列Q
dvec <- c(50,25) #目的関数のベクトルc
Amat <- matrix(c(1,0,0,1),2,2) #制約条件行列
bvec <- c(1,2) #制約条件のベクトルb
#solve.QP()に対応するように式変形
Dmat <- 2*Dmat
dvec <- - dvec
Amat <- t(Amat)
#二次計画法をを適用
result <- solve.QP(Dmat,dvec,Amat,bvec=bvec)
#解の確認
resultsolve.QP()を使用する際の注意点
solve.QP()の使う際には、いくつかの注意点があります。
まず目的関数は、最小化を目的として
$ 1/2 x ^ {T} Q x – c ^ {T} x $
のように定義されています。
そのため、
目的関数 $ \begin {pmatrix} x_1 & x_2 \\ \end {pmatrix} \begin {bmatrix} 5 & 0 \\ 0 & 10 \\ \end {bmatrix} \begin {pmatrix} x_1 \\ x_2 \end {pmatrix} + \begin {pmatrix} 50 & 25 \end {pmatrix} \begin {pmatrix} x_1 \\ x_2 \end {pmatrix} $
の例では、目的関数を下のように整理する必要があります。
$1/2 \begin {pmatrix} x_1 & x_2 \\ \end {pmatrix} \begin {bmatrix} 10 & 0 \\ 0 & 20 \\ \end {bmatrix} \begin {pmatrix} x_1 \\ x_2 \end {pmatrix} – \begin {pmatrix} -50 & -25 \end {pmatrix} \begin {pmatrix} x_1 \\ x_2 \end {pmatrix} $
1/2をしている理由は、$x_1 x_2$をのような項を考える時に都合がよいからです。
また、目的関数の中に定数項がある場合は、定数項は解に影響を及ぼさないので省略して大丈夫です。
制約条件は
$ A^ {T} x ≧ b $
と定義されているため、必要な場合は転置を忘れないようにしてください。
上記の注意点は、ソースコードの式変形部分に盛り込まれています。
結果の見方
$x$の解の確認
> result$solution
[1] 1 2目的関数の最小値を確認
> result$value
[1] 145目的関数の中に定数項がある場合は、この値にプラスしてください。
最適化の中では、定数項は解に影響を及ぼさないので、省略されています。
Rのおすすめの本をご紹介
Rでのデータ分析を体系的に学べるおすすめの本です。
ちょっと古いですが、Rの有名な本です。
まとめ
いかがだったでしょうか。私は、学生時代に二次計画問題を解いていたので、Rでの解き方についてまとめてみました。
同じような方の役に立てれば幸いです。
また、ミスがあれば修正しますので、ツイッターなどで教えていただけるとありがたいです。
コメント