Rでラグランジュの未定乗数法

ラグランジュの未定乗数法 とは

解きたい問題

等式制約のもとで関数の最大値・最小値を求める問題を考えます。

\[
\begin{aligned}
&\text{最大化(または最小化)} && f(\boldsymbol{x}) \\
&\text{制約条件} && g(\boldsymbol{x}) = c
\end{aligned}
\]

ここで、

• \(\boldsymbol{x} = (x_1, x_2, \dots, x_n)\)
  • 変数のベクトル
• \(f(\boldsymbol{x})\)
  • 最大化(または最小化)したい目的関数
• \(g(\boldsymbol{x})\)
  • 制約を表す関数
• \(c\)
  • 制約の値を表す定数

です。

制約が1本で変数が2つの場合、制約を \(y = h(x)\) の形に解いて \(f(x, h(x))\) を1変数関数として微分する方法(代入法)が使えます。

しかし、制約が複雑で解けない場合や、変数・制約が多い場合には代入法は現実的ではありません。

ラグランジュの未定乗数法は、制約を解かずに扱える一般的な方法になります。

幾何学的な考え方

\(f\) の等高線(値が一定の曲線)と、制約曲線 \(g(\boldsymbol{x}) = c\) を同じ平面に描いて考えます。

制約曲線の上を動きながら \(f\) を大きくしていくと、\(f\) の等高線が制約曲線に接する点で、それ以上 \(f\) を大きくできなくなります。

もし等高線が制約曲線を横切っていれば、制約曲線に沿って動くことで \(f\) をさらに増減できるためです。

等高線と制約曲線が接するとは、両者の法線ベクトル(曲線(曲面)の接線に直交する向きのベクトル)が平行になることです。

関数の勾配ベクトルは等高線に直交しますので、

\[
\nabla f(\boldsymbol{x}) = \lambda \, \nabla g(\boldsymbol{x})
\]

となる実数 \(\lambda\) が存在します。ここで、

• \(\nabla f\)
  • \(f\) の勾配ベクトル \(\left(\dfrac{\partial f}{\partial x_1}, \dots, \dfrac{\partial f}{\partial x_n}\right)\)
• \(\nabla g\)
  • \(g\) の勾配ベクトル \(\left(\dfrac{\partial g}{\partial x_1}, \dots, \dfrac{\partial g}{\partial x_n}\right)\)
• \(\lambda\)
  • ラグランジュ乗数(2本の勾配ベクトルの比例定数)

です。

ラグランジュ関数と停留条件

上の条件と制約条件をまとめて1つの式で表すために、ラグランジュ関数を次のように定義します。

\[
L(\boldsymbol{x}, \lambda) = f(\boldsymbol{x}) - \lambda \, \bigl( g(\boldsymbol{x}) - c \bigr)
\]

\(L\) を \(\boldsymbol{x}\) と \(\lambda\) のすべてで偏微分して 0 とおくと、

\[
\begin{aligned}
\frac{\partial L}{\partial x_i} &= \frac{\partial f}{\partial x_i} - \lambda \frac{\partial g}{\partial x_i} = 0 \quad (i = 1, \dots, n) \\
\frac{\partial L}{\partial \lambda} &= -\bigl( g(\boldsymbol{x}) - c \bigr) = 0
\end{aligned}
\]

となります。

上段の \(n\) 本の式が「勾配が平行」という条件 \(\nabla f = \lambda \nabla g\) であり、下段の式がもとの制約条件 \(g(\boldsymbol{x}) = c\) です。

つまり、制約付き問題を、\(n + 1\) 個の未知数 \((\boldsymbol{x}, \lambda)\) に関する \(n + 1\) 本の連立方程式に置き換えることができます。

その連立方程式を解いて得られる点を停留点(候補点)と呼びます。

制約が複数ある場合

制約が \(m\) 本ある場合は、制約ごとに乗数を用意します。

\[
L(\boldsymbol{x}, \boldsymbol{\lambda}) = f(\boldsymbol{x}) - \sum_{j=1}^{m} \lambda_j \bigl( g_j(\boldsymbol{x}) - c_j \bigr)
\]

ここで、

• \(m\)
  • 制約の本数
• \(g_j(\boldsymbol{x}) = c_j\)
  • \(j\) 番目の制約条件
• \(\lambda_j\)
  • \(j\) 番目の制約に対応するラグランジュ乗数
• \(\boldsymbol{\lambda} = (\lambda_1, \dots, \lambda_m)\)
  • 乗数のベクトル

です。未知数は \(n + m\) 個、方程式も \(n + m\) 本になります。

ラグランジュ乗数の意味

最適値を制約の値 \(c\) の関数として \(f^{*}(c)\) と書くと、

\[
\frac{d f^{*}(c)}{d c} = \lambda
\]

が成り立ちます(包絡線定理)。ここで、

• \(f^{*}(c)\)
  • 制約の値が \(c\) のときの \(f\) の最適値
• \(\lambda\)
  • 最適解におけるラグランジュ乗数

です。

すなわち \(\lambda\) は、制約を少し緩めたとき(たとえば予算 \(c\) を1単位増やしたとき)に最適値がどれだけ変化するかを表す感度になります。

注意点

• 必要条件であって十分条件ではない
  • 連立方程式の解は極値の候補です。最大か最小か、あるいはどちらでもないかは、候補点での \(f\) の値の比較や、2階の条件(縁付きヘッセ行列)で判断します。
• 制約の勾配が 0 でないこと
  • 最適解で \(\nabla g(\boldsymbol{x}) \neq \boldsymbol{0}\) でなければ、上の方法で最適解を見つけられないことがあります(制約想定)。
• 等式制約が対象
  • 不等式制約 \(g(\boldsymbol{x}) \leq c\) には、ラグランジュの未定乗数法を拡張したKKT条件を使います。

Rコード

取り上げる例題

次の問題を扱います。

\[
\begin{aligned}
&\text{最大化・最小化} && f(x, y) = x y \\
&\text{制約条件} && g(x, y) = x^2 + y^2 = 1
\end{aligned}
\]

ここで、

• \(x, y\)
  • 変数
• \(f(x, y) = xy\)
  • 目的関数
• \(g(x, y) = x^2 + y^2\)
  • 制約を表す関数(単位円)
• 制約の値
  • \(c = 1\)

です。単位円の上で積 \(xy\) の最大値・最小値を求める問題です。

手計算による解

ラグランジュ関数は

\[
L(x, y, \lambda) = x y - \lambda \, (x^2 + y^2 - 1)
\]

です。ここで \(\lambda\) はラグランジュ乗数です。停留条件は

\[
\begin{aligned}
&\frac{\partial L}{\partial x} = y - 2 \lambda x = 0 \\
&\frac{\partial L}{\partial y} = x - 2 \lambda y = 0 \\
&x^2 + y^2 = 1
\end{aligned}
\]

となります。

第1式から \(y = 2\lambda x\) となりますので、第2式に代入すると \(x = 4 \lambda^2 x\) です。

\(x = 0\) の場合 \(y = 0\) となり制約を満たしませんので \(x \neq 0\) であり、\(\lambda^2 = 1/4\)、すなわち \(\lambda = \pm 1/2\) を得ます。

• \(\lambda = +1/2\) のとき \(y = x\) となり、\(x^2 + y^2 = 1\) から \((x, y) = \pm\left(\tfrac{1}{\sqrt{2}}, \tfrac{1}{\sqrt{2}}\right)\)、\(f = \tfrac{1}{2}\) (最大値)
• \(\lambda = -1/2\) のとき \(y = -x\) となり、\((x, y) = \pm\left(\tfrac{1}{\sqrt{2}}, -\tfrac{1}{\sqrt{2}}\right)\)、\(f = -\tfrac{1}{2}\) (最小値)

以降のRコードで、この結果を数値計算と図で確かめます。

問題の定義

# 目的関数
f <- function(x, y) x * y

# 制約を表す関数
g <- function(x, y) x^2 + y^2

# 制約の値
c0 <- 1

# 目的関数の勾配ベクトル
grad_f <- function(x, y) c(y, x)

# 制約を表す関数の勾配ベクトル
grad_g <- function(x, y) c(2 * x, 2 * y)

等高線と制約曲線を描く

まず、\(f\) の等高線と制約曲線(単位円)を重ねて描きます。

円の上を動いたときに \(f\) が最大・最小になる点は、等高線が円に接する点になります。

library(ggplot2)

# 等高線を描くための格子データ
grid_df <- expand.grid(
  x = seq(-1.5, 1.5, length.out = 301),
  y = seq(-1.5, 1.5, length.out = 301)
) |>
  transform(z = f(x, y))

# 制約曲線(単位円)のデータ
theta_seq <- seq(0, 2 * pi, length.out = 401)
circle_df <- data.frame(
  x = sqrt(c0) * cos(theta_seq),
  y = sqrt(c0) * sin(theta_seq)
)

p_base <- ggplot() +
  geom_contour(
    data = grid_df,
    aes(x = x, y = y, z = z, colour = after_stat(level)),
    breaks = seq(-1.2, 1.2, by = 0.2)
  ) +
  geom_path(data = circle_df, aes(x = x, y = y), linewidth = 1.2) +
  scale_colour_gradient2(
    low = "blue", mid = "grey70", high = "red", midpoint = 0
  ) +
  coord_equal() +
  labs(
    title = "f(x, y) = xy の等高線と制約曲線 x² + y² = 1",
    x = "x",
    y = "y",
    colour = "f の値"
  ) +
  theme_bw()

p_base

黒い太線が制約曲線、色付きの細線が \(f\) の等高線です。

等高線と制約曲線が接すること無く横切っている場合、制約曲線に沿って動くことで \(f\) をさらに増減できてしまうことが確認できます。

停留条件の連立方程式をニュートン法で解く

停留条件は、未知数 \((x, y, \lambda)\) についての3本の連立方程式

\[
\boldsymbol{F}(x, y, \lambda) =
\begin{pmatrix}
y - 2 \lambda x \\
x - 2 \lambda y \\
x^2 + y^2 - 1
\end{pmatrix}
= \boldsymbol{0}
\]

です。

ここで \(\boldsymbol{F}\) は3つの停留条件を縦に並べたベクトル値関数です(第3成分は \(\partial L / \partial \lambda\) の符号を反転したもので、解は変わりません)。

ニュートン法では、現在の値 \(\boldsymbol{v}_k = (x_k, y_k, \lambda_k)\) から次の値を

\[
\boldsymbol{v}_{k+1} = \boldsymbol{v}_k - J(\boldsymbol{v}_k)^{-1} \, \boldsymbol{F}(\boldsymbol{v}_k)
\]

で更新します。ここで \(J\) は \(\boldsymbol{F}\) のヤコビ行列

\[
J(x, y, \lambda) =
\begin{pmatrix}
-2\lambda & 1 & -2x \\
1 & -2\lambda & -2y \\
2x & 2y & 0
\end{pmatrix}
\]

です。

# 停留条件を並べたベクトル値関数
stationary_fn <- function(v) {
  x <- v[1]
  y <- v[2]
  lambda <- v[3]
  c(
    y - 2 * lambda * x,
    x - 2 * lambda * y,
    g(x, y) - c0
  )
}

# ヤコビ行列
stationary_jac <- function(v) {
  x <- v[1]
  y <- v[2]
  lambda <- v[3]
  matrix(
    c(
      -2 * lambda, 1, -2 * x,
      1, -2 * lambda, -2 * y,
      2 * x, 2 * y, 0
    ),
    nrow = 3,
    byrow = TRUE
  )
}

# ニュートン法(収束しなければ NA を返す)
newton_solve <- function(v, tol = 1e-12, max_iter = 100) {
  for (k in seq_len(max_iter)) {
    fv <- stationary_fn(v)
    if (max(abs(fv)) < tol) {
      return(v)
    }
    v <- tryCatch(
      v - solve(stationary_jac(v), fv),
      error = function(e) rep(NA_real_, 3)
    )
    if (anyNA(v)) {
      break
    }
  }
  rep(NA_real_, 3)
}

# 乱数で初期値を50組つくり、それぞれからニュートン法を実行する
set.seed(20261021)
start_list <- lapply(seq_len(50), function(i) runif(3, min = -1.5, max = 1.5))
sol_mat <- do.call(rbind, lapply(start_list, newton_solve))
colnames(sol_mat) <- c("x", "y", "lambda")

# 収束した解を丸めて重複を除く
stat_df <- sol_mat |>
  as.data.frame() |>
  na.omit() |>
  round(6) |>
  unique()

# 目的関数の値と、最大・最小の区別を加える
stat_df$f_value <- f(stat_df$x, stat_df$y)
stat_df$type <- ifelse(
  stat_df$f_value == max(stat_df$f_value), "最大",
  ifelse(stat_df$f_value == min(stat_df$f_value), "最小", "")
)

# 表示用に並べ替える
stat_df <- stat_df[order(-stat_df$f_value, stat_df$x), ]
rownames(stat_df) <- NULL

knitr::kable(
  stat_df,
  col.names = c("x", "y", "λ", "f の値", "判定"),
  digits = 4
)

# 1/√2の確認
1 / sqrt(2)

[1] 0.7071068

乱数で選んだ初期値から出発しても、停留点は4個に集まります。

得られた \(x, y, \lambda, f\) の値は、手計算の結果(\(\lambda = \pm 1/2\)、\(f = \pm 1/2\))と一致しています。

勾配ベクトルが平行になっていることを確かめる

ラグランジュの未定乗数法の要点は、最適解で \(\nabla f = \lambda \nabla g\) が成り立つことでした。

各停留点で \(\nabla f\) と \(\nabla g\) を計算し、確かめます。

# 各停留点で勾配ベクトルを計算する
grad_df <- do.call(rbind, lapply(seq_len(nrow(stat_df)), function(i) {
  x <- stat_df$x[i]
  y <- stat_df$y[i]
  lambda <- stat_df$lambda[i]
  gf <- grad_f(x, y)
  gg <- grad_g(x, y)
  data.frame(
    x = x,
    y = y,
    fx = gf[1],
    fy = gf[2],
    gx = gg[1],
    gy = gg[2],
    # lambda 倍した制約の勾配との差(0 になれば平行)
    diff_norm = sqrt(sum((gf - lambda * gg)^2))
  )
}))

knitr::kable(
  grad_df,
  col.names = c("x", "y", "∂f/∂x", "∂f/∂y", "∂g/∂x", "∂g/∂y", "|∇f - λ∇g|"),
  digits = 6
)

最後の列 \(|\nabla f - \lambda \nabla g|\) は、4つの停留点のすべてで 0 になります。

等高線の図に、各停留点での勾配ベクトルを重ねます。

赤の矢印が \(\nabla f\)、青の矢印が \(\nabla g\) です(見やすいように長さを 0.3 倍にしています)。

# 矢印の長さを調整する倍率
arrow_scale <- 0.3

p_base +
  geom_segment(
    data = grad_df,
    aes(x = x, y = y, xend = x + arrow_scale * fx, yend = y + arrow_scale * fy),
    colour = "red",
    linewidth = 0.9,
    arrow = arrow(length = unit(0.15, "cm"))
  ) +
  geom_segment(
    data = grad_df,
    aes(x = x, y = y, xend = x + arrow_scale * gx, yend = y + arrow_scale * gy),
    colour = "blue",
    linewidth = 0.9,
    linetype = "dashed",
    arrow = arrow(length = unit(0.15, "cm"))
  ) +
  geom_point(data = grad_df, aes(x = x, y = y), size = 2.5) +
  labs(
    title = "停留点での勾配ベクトル(赤: ∇f、青の破線: ∇g)",
    colour = "f の値"
  )

4つの停留点のどこでも、赤と青の矢印が同じ直線上に並んでいます。

\(\lambda > 0\) の点(最大点)では同じ向き、\(\lambda < 0\) の点(最小点)では逆向きです。

こちらが「\(\nabla f = \lambda \nabla g\)」の幾何学的な意味になります。

制約曲線上で f を1変数関数として見る

制約曲線は \(x = \cos\theta,\ y = \sin\theta\) (\(\theta\) は角度)とパラメータ表示できます。

このとき \(f = \cos\theta \sin\theta = \tfrac{1}{2}\sin 2\theta\) となり、制約付きの問題が \(\theta\) の1変数関数の極値問題になります。

# 制約曲線上の f を角度 theta の関数として計算する
theta_df <- data.frame(theta = theta_seq) |>
  transform(f_value = f(cos(theta), sin(theta)))

# 停留点の角度(0 から 2π の範囲にそろえる)
stat_df$theta <- atan2(stat_df$y, stat_df$x) %% (2 * pi)

ggplot(theta_df, aes(x = theta, y = f_value)) +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_line(linewidth = 1) +
  geom_point(
    data = stat_df,
    aes(x = theta, y = f_value, colour = type),
    size = 3
  ) +
  scale_x_continuous(
    breaks = seq(0, 2 * pi, by = pi / 2),
    labels = c("0", "π/2", "π", "3π/2", "2π")
  ) +
  labs(
    title = "制約曲線上での f の値(θ は角度)",
    x = "θ",
    y = "f の値",
    colour = "判定"
  ) +
  theme_bw()

ニュートン法で求めた4つの停留点が、この曲線の山の頂点・谷の底とちょうど重なっています。

ラグランジュの未定乗数法は、制約を解かなくても、この1変数化と同じ候補点を探索できます。

ラグランジュ乗数が表す感度を確かめる

最後に、\(\lambda\) が「制約の値 \(c\) を変えたときの最適値の変化率」であることを数値で確かめます。

制約を \(x^2 + y^2 = c\) に変えたときの最大値を \(f^{*}(c)\) とします。

円の半径は \(\sqrt{c}\) ですので、\(x = \sqrt{c}\cos\theta,\ y = \sqrt{c}\sin\theta\) として \(\theta\) について optimize() で最大化します。

# 制約の値が cc のときの最大値を数値的に求める
f_star <- function(cc) {
  radius <- sqrt(cc)
  optimize(
    function(theta) f(radius * cos(theta), radius * sin(theta)),
    interval = c(0, pi / 2),
    maximum = TRUE,
    tol = 1e-10
  )$objective
}

# 差分近似で f*(c) の c に関する微分を求める
step_h <- 1e-4
df_dc <- (f_star(c0 + step_h) - f_star(c0 - step_h)) / (2 * step_h)

# 最大点に対応するラグランジュ乗数
lambda_max <- stat_df$lambda[stat_df$type == "最大"][1]

data.frame(
  項目 = c("df*/dc の数値微分", "最大点のラグランジュ乗数 λ"),
  値 = c(df_dc, lambda_max)
) |>
  knitr::kable(digits = 6)

数値微分で求めた \(d f^{*}/dc\) は、最大点の \(\lambda = 1/2\) と一致します。

手計算でも、\(x^2 + y^2 = c\) のときの最大値は \(f^{*}(c) = c/2\) ですので、\(d f^{*}/dc = 1/2 = \lambda\) です。

まとめ

• 制約付きの最適化は、ラグランジュ関数 \(L = f - \lambda (g - c)\) の停留条件(連立方程式)に置き換えられる。
• 最適解では、等高線と制約曲線が接し、\(\nabla f\) と \(\nabla g\) が平行になる。
• 得られた停留点は候補であり、最大・最小の判定は \(f\) の値の比較などで別に行う。
• \(\lambda\) は、制約を緩めたときの最適値の変化率を表す。

以上です。