Rの関数:constrOptim {stats}

Rの関数から constrOptim {stats} を確認します。

関数 constrOptim とは

constrOptim 関数は、Adaptive Barrier Algorithm を利用して、線形不等式制約の条件下で目的関数を最小化（または最大化）するための最適化手法を提供する関数です。

指定された制約領域の境界に近づくにつれて、目的関数に対してペナルティとなる barrier term (以降、本投稿では障壁項)を加算することにより、探索プロセスが制約領域の外部へ逸脱しないように制御します。

最適化の内部処理においては、 optim {stats} を反復的に呼び出して計算を進めます。

制約を持たない単純な最適化を実施する optim と比較して、本関数はパラメータが特定の線形不等式を満たす必要がある場合に有用です。

本関数のソースコードを確認しますと、関数内部において R と dR という内部関数が定義されています。

当該内部関数は、本来の目的関数 f に対して、障壁項を加算する役割を担っています。

当該障壁項は、パラメータが制約の境界へ近づくにつれてペナルティを増大させ、探索が実行可能領域の外部へ逸脱する事態を未然に防ぎます。

外側の反復ループが進むたびに、前回探索されたパラメータ theta.old を基準として適応的にペナルティの形状が更新され、目的関数が着実に減少するよう保証する仕組みが組み込まれています。

探索が停滞するか反復上限に到達した際には、処理を終了し、最適解と各種収束情報を返す設計となっております。

関数 constrOptim の引数

args(constrOptim)

function (theta, f, grad, ui, ci, mu = 1e-04, control = list(), 
    method = if (is.null(grad)) "Nelder-Mead" else "BFGS", outer.iterations = 100, 
    outer.eps = 1e-05, ..., hessian = FALSE) 
NULL

theta
- 最適化を開始するための初期値となる数値ベクトルです。当該値は、指定されたすべての線形不等式制約を厳密に満たしている、すなわち制約領域の内部に完全に位置している必要があります。境界線上にある初期値を利用することはできません。
f
- 最適化の対象となる目的関数を指定します。当該関数は最初の引数としてパラメータベクトルを受け取り、スカラー値を返すように設計されている必要があります。
grad
- 目的関数の勾配（一階導関数）を計算する関数を指定します。こちらも最初の引数としてパラメータベクトルを受け取り、同パラメータと同等の長さのベクトルを返す必要があります。method 引数が "Nelder-Mead" に設定されている場合に限り、当該引数へ NULL を指定することが可能です。
ui
- 線形不等式制約の係数を示す行列を指定します。制約式の左辺を構成し、内部的には ui %*% theta - ci >= 0 という条件として評価されます。
ci
- 線形不等式制約の定数項を示すベクトルを指定します。前述した不等式の右辺に関連する値となります。
mu
- 障壁関数に対するペナルティの強さを調整するためのチューニングパラメータです。デフォルト値は 1e-04 となっております。当該値が大きくなるにつれて、制約の境界付近における目的関数の変化が急峻になります。また、control$fnscale が負の値（最大化問題）に設定されている場合は、内部の処理によって自動的に符号が反転する仕組みが備わっております。
control
- 内部で呼び出される optim 関数に対して渡される制御用パラメータのリストを指定します。最大化問題として解きたい場合は、当該リスト内に fnscale = -1 と設定します。
method
- 内部の optim 関数で使用される最適化アルゴリズムを指定する文字列です。grad 引数に NULL が指定された場合は自動的に "Nelder-Mead" が選択され、勾配関数が与えられた場合はデフォルトで "BFGS" が選ばれます。
outer.iterations
- 外側ループの最大反復回数を指定する整数値です。デフォルト値は 100 です。当該ループの中で適応的にペナルティが更新され、反復の都度 optim 関数が実行されます。
outer.eps
- 外側ループにおける相対収束許容誤差を示す非負の数値を指定します。デフォルト値は 1e-05 です。連続する反復間における目的関数値の変化量が当該許容誤差を下回った際、最適化が収束したと見なされて計算が終了します。
...
- 目的関数 f や勾配関数 grad、ならびに内部の optim に対して追加で渡すべき任意の引数を指定します。
hessian
- 最適解におけるヘッセ行列（二次導関数の行列）を計算し、出力結果に含めるかどうかを指定する論理値です。デフォルト値は FALSE となっております。

シミュレーションコード

関数 constrOptim の機能を理解するためのシミュレーションを実施します。

本シミュレーションでは、ある仮想的な2種類の物質（物質Aおよび物質B）の混合比率を決定するプロセスを想定します。

目的は、当該2つの物質を混合する際に発生する「エネルギーコスト関数」を最小化することです。

エネルギーコスト関数は二次関数として定義され、物質の混合には以下の3つの線形不等式制約が課されているものとします。

物質の合計量は4以下でなければならない（-x - y >= -4）
物質Aの量は物質Bの量以上でなければならない（x - y >= 0）
物質Bの量は0以上でなければならない（y >= 0）

エネルギーコスト関数の理論上の最小値は物質Aが3、物質Bが2の場合ですが、当該組み合わせは制約条件を満たしません。それゆえ、関数 constrOptim を用いて、制約領域内における最適な混合比率を探索します。

以下にシミュレーション用のRコードを記します。

# 必要なパッケージの読み込み
library(ggplot2)

# 目的関数の定義（エネルギーコスト関数）
# 物質Aの量をx、物質Bの量をyとします
cost_function <- function(theta) {
  x <- theta[1]
  y <- theta[2]
  (x - 3)^2 + (y - 2)^2
}

# 目的関数の勾配関数の定義
grad_cost_function <- function(theta) {
  x <- theta[1]
  y <- theta[2]
  c(2 * (x - 3), 2 * (y - 2))
}

# 制約条件の設定
# 条件1: -x - y >= -4 (合計量は4以下)
# 条件2: x - y >= 0  (物質Aの量は物質Bの量以上)
# 条件3: y >= 0      (物質Bの量は0以上)
ui_matrix <- matrix(c(-1, -1,
                       1, -1,
                       0,  1), ncol = 2, byrow = TRUE)
ci_vector <- c(-4, 0, 0)

# 初期値の設定
# 制約領域の内部にある値を指定する必要があります。
# x = 1, y = 0.5 を初期値とします。
initial_values <- c(1, 0.5)

# 関数 constrOptim の実行
optim_result <- constrOptim(
  theta = initial_values,
  f = cost_function,
  grad = grad_cost_function,
  ui = ui_matrix,
  ci = ci_vector
)

# 最適解の取得
optimal_values <- optim_result$par

# 結果の出力
cat("最適解における物質Aの量: ", optimal_values[1], "\n")
cat("最適解における物質Bの量: ", optimal_values[2], "\n")
cat("最小化されたエネルギーコスト: ", optim_result$value, "\n")

最適解における物質Aの量:  2.499992 
最適解における物質Bの量:  1.500008 
最小化されたエネルギーコスト:  0.5

エネルギーコストを最小化する最適な混合比率は、物質Aが約2.5、物質Bが約1.5であり、到達した最小エネルギーコストは0.5であるとの結果となりました。

目的関数であるエネルギーコスト関数は、物質Aの量が3、物質Bの量が2の際に理論上の最小値0をとるように定義されています。

しかしながら、「物質の合計量は4以下でなければならない」という第一の制約条件が存在するため、理論上最適な組み合わせである合計量5（= 3＋2）は、許容される領域の完全に外部へ位置しています。

それゆえ、最適化アルゴリズムは制約領域の内部ならびに境界において、理論上の最適座標（3, 2）へ最も近い地点を探索する手順へ移行します。

幾何学的には、制約の境界線である「物質A ＋物質B ＝ 4」の直線上で、座標（3, 2）からの距離が最短となる直交点は（2.5, 1.5）となります。

当該地点におけるエネルギーコストを実際に計算いたしますと、(2.5 - 3)^2 + (1.5 - 2)^2 = 0.25 + 0.25 = 0.5 となり、シミュレーションが提示した最小値と理論値が一致していることを確認できます。

# シミュレーション結果の可視化のためのデータ作成
# 等高線を描画するためのグリッドデータ
grid_data <- expand.grid(
  x = seq(0, 5, length.out = 100),
  y = seq(0, 5, length.out = 100)
)
grid_data$cost <- with(grid_data, (x - 3)^2 + (y - 2)^2)

# 制約領域を示す多角形(三角形)の頂点データ
# 各方程式の交点から頂点を算出
polygon_data <- data.frame(
  x = c(0, 4, 2),
  y = c(0, 0, 2)
)

# ggplot2による可視化
ggplot() +
  geom_contour_filled(data = grid_data, aes(x = x, y = y, z = cost), alpha = 0.8) +
  geom_polygon(data = polygon_data, aes(x = x, y = y), 
               fill = "white", alpha = 0.5, color = "black", linewidth = 1) +
  geom_point(aes(x = initial_values[1], y = initial_values[2]), 
             color = "blue", size = 4) +
  geom_point(aes(x = optimal_values[1], y = optimal_values[2]), 
             color = "red", size = 4) +
  geom_segment(aes(x = initial_values[1], y = initial_values[2], 
                   xend = optimal_values[1], yend = optimal_values[2]), 
               arrow = arrow(length = unit(0.3, "cm")), color = "black", linewidth = 1) +
  annotate("text", x = initial_values[1], y = initial_values[2] - 0.2, 
           label = "初期値", color = "blue", size = 5) +
  annotate("text", x = optimal_values[1], y = optimal_values[2] + 0.2, 
           label = "最適解", color = "red", size = 5) +
  labs(title = "エネルギーコスト関数の等高線と制約領域における最適化",
       x = "物質Aの量", 
       y = "物質Bの量", 
       fill = "エネルギーコスト") +
  theme_minimal()