ラグランジュの未定乗数法を理解できた気がするのでまとめる

2019-05-13

ラグランジュ(ラグランジェ)の未定乗数法は 一つ以上の束縛条件(制約関数)のもとで、目的関数の極値を求めるための最適化手法です。

制約条件が同値の場合にしか使えませんが、 KKT条件を用いることにより不等号を含む制約条件を扱えるようになります。

これについてはまた別の機会に記事にしたいと思います。

備考

同僚の Tsutomu Saito さんにレビューしていただきました。 斎藤さんだぞ!(これを言うのも最後かな..

ざっとみて違和感はないと言っていただきましたが、細かいところで変なこと書いてたら指摘ください。

目次

手順

目的関数 f と 制約関数 g がある前提で話を進めます。 なお、この記事では考えやすいように 2変数関数 の目的関数と制約関数を例に考えます。

制約関数 g に未定定数と呼ばれる変数(一般的には λ とする)をかけ、目的関数から引いたものを新たな関数 L とする

$$ L(x, y, λ) = f(x, y) - λg(x, y) $$

L をそれぞれの引数 で偏微分したものが 0 になるよう式を連立して解く

$$ \begin{eqnarray} \left\{ \begin{array}{l} \frac{∂L}{∂x} = \frac{∂f}{∂x} - λ\frac{∂g}{∂x} = 0 \\ \frac{∂L}{∂y} = \frac{∂f}{∂y} - λ\frac{∂g}{∂y} = 0 \\ \frac{∂L}{∂λ} = g(x, y) = 0 \\ \end{array} \right. \end{eqnarray} $$
In [ ]:
 

この連立方程式をとけば、極値の座標が求められます。

考え方

なぜこのような手順で極値が求められるのかを考えてみましょう。

  • 目的関数 f(x, y) は 単なる2変数関数なので です
  • 制約関数 g(x, y) = 0線(直線または曲線) です
    • (関数の結果を Z軸とすると)制約関数の Z=0のXY平面 と交わっている部分です
    • この線というのが、x, y の定義域ということになります。

基本例題

線を面に沿って描画してみたのが以下です。

In [1]:
%matplotlib notebook
In [2]:
import numpy as np
import sympy as sy

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
In [3]:
sy.init_printing()
In [4]:
x, y, l = sy.symbols('x y λ')
In [5]:
def purpose(x, y):
    return x ** 2 + y ** 2 + 3
In [6]:
def restrict(x, y):
    return x * 2 + 10 * y - 15
In [7]:
def gradient(x):
    return (-2 * x + 15) / 10
In [8]:
fx = purpose(x, y)
In [9]:
gx = restrict(x, y)
In [10]:
fx
Out[10]:
$\displaystyle x^{2} + y^{2} + 3$
In [11]:
gx
Out[11]:
$\displaystyle 2 x + 10 y - 15$
In [12]:
L = fx - l * gx
L
Out[12]:
$\displaystyle x^{2} + y^{2} - λ \left(2 x + 10 y - 15\right) + 3$
In [13]:
dx = L.diff(x)
dx
Out[13]:
$\displaystyle 2 x - 2 λ$
In [14]:
dy = L.diff(y)
dy
Out[14]:
$\displaystyle 2 y - 10 λ$
In [15]:
dl = L.diff(l)
dl
Out[15]:
$\displaystyle - 2 x - 10 y + 15$
In [16]:
p = sy.solve([dx, dy, dl])
p
Out[16]:
$\displaystyle \left\{ x : \frac{15}{52}, \ y : \frac{75}{52}, \ λ : \frac{15}{52}\right\}$
In [17]:
xs = np.arange(-3, 3, 0.1)
ys = np.arange(-3, 3, 0.1)
In [18]:
X, Y = np.meshgrid(xs, ys)
In [19]:
fig = plt.figure()
ax = Axes3D(fig)

ax.plot_surface(X, Y, purpose(X, Y), color='blue', alpha=0.5)
ax.plot(xs, gradient(xs), [purpose(x, gradient(x)) for x in xs], color="red")
ax.scatter([p[x]], [p[y]], [purpose(p[x], p[y])], color='orange')
plt.show()
In [20]:
fig = plt.figure()
ax = Axes3D(fig)

ax.contour(X, Y, purpose(X, Y), alpha=0.5)
ax.plot(xs, gradient(xs), [purpose(x, gradient(x)) for x in xs], color="red")
ax.scatter([p[x]], [p[y]], [purpose(p[x], p[y])], color='orange')
plt.show()

最後の図は 曲面を 等高線で表示してみました。次のセクションで使います。

式の意味

さて、ここからは式の意味を考えていきましょう。

一言で言うなら、求めるべき極値は 「 等高線接線 (接点)」です。 (もしくは共通の法線を持つ点)

等高線の 交点 は制約条件は満たしていますが 極値ではありません 。(上の図を見てもわかると思います)

備考

等高線とは同じ高さの座標(点)をつないだ線のことです。

実は等高線というものは無数に書くことができますが、 それでは非常に見づらいので普通は切りの良い単位で間隔を設定します.

今回は天下り的に接している位置に等高線がくる感じになりました。

つまり、接線の座標を求めれば良いのですね。では接線とはなんでしょうか?

接線 - Wikipedia より

曲線の接線(せっせん、英: tangent line, tangent)は、平面曲線に対しては、曲線上の一点が与えられたとき、その点において曲線に「ただ触れるだけ」の直線を意味する。ライプニッツは接線を、曲線上の無限に近い二点を通る直線として定義した[1]。 より具体的に解析幾何学において、与えられた直線が曲線 y = f(x) の x = c(あるいは曲線上の点 (c, f(c))における接線であるとは、その直線が曲線上の点 (c, f (c)) を通り、傾きが f の微分係数 f'(c) に等しいときに言う。 同様の定義は空間曲線やより高次のユークリッド空間内の曲線に対しても適用できる。

曲線と接線が相接する点は接点 (point of tangency) と言い、曲線との接点において接線は曲線と「同じ方向へ」進む。 その意味において接線は、接点における曲線の最適直線近似である。

備考

等高線なので当然ですが、XY平面における傾きを指しています。

簡単に言うと今回の3次元空間を上(Z方向)から覗き込んだ図です。

以下の 2条件を満たせば接線であるといえそうです。 それぞれの条件を表す3式も併せて記載します。これらを連立させて導いた座標が求めるべき座標となります。

等高線が平行

L を x で偏微分 した結果が 0 -- (1)

$$ \frac{∂L}{∂x} = 0 $$

L を展開すると

$$ \frac{∂f}{∂x} - λ\frac{∂g}{∂x} = 0\\ \frac{∂f}{∂x} = λ\frac{∂g}{∂x} $$

f の偏微分(x)は g の偏微分(x)の λ倍 であることがわかる

L を y で偏微分 した結果が 0 -- (2)

$$ \frac{∂L}{∂y} = 0 $$

L を展開すると

$$ \frac{∂L}{∂y} = \frac{∂f}{∂y} - λ\frac{∂g}{∂y} = 0 \\ \frac{∂f}{∂y} = λ\frac{∂g}{∂y} $$

f の偏微分(y)は g の偏微分(y)の λ倍 であることがわかる

ベクトルで表現すると $$ \left[ \begin{array}{rrr} \frac{∂f}{∂x} \\ \frac{∂f}{∂y} \end{array} \right] = \left[ \begin{array}{rrr} λ\frac{∂g}{∂x} \\ λ\frac{∂g}{∂y} \end{array} \right] $$

λ を括って以下のようになる

$$ \left[ \begin{array}{rrr} \frac{∂f}{∂x} \\ \frac{∂f}{∂y} \end{array} \right] = λ\left[ \begin{array}{rrr} \frac{∂g}{∂x} \\ \frac{∂g}{∂y} \end{array} \right] $$

もう一方のベクトルのスカラー(λ)倍で表現できているので、これらは平行だと言えるのですね。

未定乗数として定義した λ は 目的関数と制約関数の接線倍率だったのですね

制約関数上に共有点を持つ

座標が離れていては接しているとは言えないので、共有点をもつ必要があります。

L を λ で偏微分 した結果が 0 -- (3)

$$ \frac{∂L}{∂λ} = 0 $$

(偏微分により)目的関数が消え、制約関数のみが取り出される

$$ g(x, y) = 0 $$

応用例題

これらを踏まえ問題を解いてみましょう。

備考

以下の問題は 別のサイト様からお借りしました。私がごちゃごちゃ解説するよりもわかりやすいと思ったからです。

載せたらまずいのがあったらお手数ですが、 問い合わせTwitterのDM まで連絡ください。

詳しい解説はリンク先を参照ください。

長径2a,短径2bの楕円に内接する長方形の面積の最大値を求めよ.

引用元リンク:

ラグランジェの未定乗数法 [物理のかぎしっぽ]

In [1]:
%matplotlib inline
In [2]:
import math
import numpy as np
import sympy as sy

from matplotlib import pyplot as plt, patches as pat
In [3]:
sy.init_printing()
In [4]:
# イメージとしてはこんな感じ
fig = plt.figure()
ax = fig.add_subplot(111)

a = 3
b = 2
sq2 = math.sqrt(2)

el = pat.Ellipse(xy=(0, 0), width=a * 2, height=b * 2, alpha=0.5)
rec = pat.Rectangle(xy=(-sq2 * a / 2, -sq2 * b / 2), width=sq2 * a, height=sq2 * b, alpha=0.5, color='pink')

ax.add_patch(el)
ax.add_patch(rec)

plt.xlim(-a, a)
plt.ylim(-b, b)
plt.grid(True, color='black')
plt.tick_params(which='major', width=1)
plt.xticks([0])
plt.yticks([0])
plt.show()

長方形の面積は 4象限の 長方体の合計から目的関数 f は $$ f(x, y) = 4xy $$

楕円の方程式から 制約関数 g は $$ g(x, y) = \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 $$ $$ g(x, y) = \frac{x^2}{a^2} + \frac{y^2}{b^2} - 1 = 0 $$

In [5]:
a, b, x, y, l = sy.symbols('a b x y λ')
In [6]:
fx = 4 * x * y
fx
Out[6]:
$\displaystyle 4 x y$
In [7]:
gx = x **2 / a ** 2 + y ** 2 / b ** 2 - 1
gx
Out[7]:
$\displaystyle -1 + \frac{y^{2}}{b^{2}} + \frac{x^{2}}{a^{2}}$
In [8]:
L = fx - l * gx
L
Out[8]:
$\displaystyle 4 x y - λ \left(-1 + \frac{y^{2}}{b^{2}} + \frac{x^{2}}{a^{2}}\right)$
In [9]:
dx = L.diff(x)
dx
Out[9]:
$\displaystyle 4 y - \frac{2 x λ}{a^{2}}$
In [10]:
dy = L.diff(y)
dy
Out[10]:
$\displaystyle 4 x - \frac{2 y λ}{b^{2}}$
In [11]:
dl = L.diff(l)
dl
Out[11]:
$\displaystyle 1 - \frac{y^{2}}{b^{2}} - \frac{x^{2}}{a^{2}}$
In [12]:
points = sy.solve([dx, dy, dl], exclude=[a, b])
points
Out[12]:
$\displaystyle \left[ \left\{ x : - \frac{\sqrt{2} a}{2}, \ y : - \frac{\sqrt{2} b}{2}, \ λ : 2 a b\right\}, \ \left\{ x : - \frac{\sqrt{2} a}{2}, \ y : \frac{\sqrt{2} b}{2}, \ λ : - 2 a b\right\}, \ \left\{ x : \frac{\sqrt{2} a}{2}, \ y : - \frac{\sqrt{2} b}{2}, \ λ : - 2 a b\right\}, \ \left\{ x : \frac{\sqrt{2} a}{2}, \ y : \frac{\sqrt{2} b}{2}, \ λ : 2 a b\right\}\right]$
In [13]:
# 面積
sum(abs(p[x] * p[y]) for p in points)
Out[13]:
$\displaystyle 2 \left|{a b}\right|$

関数 f(x,y,z) = 8xyz の最大値を,条件 g(x,y,z) = x^2+y^2+z^2-1 = 0 のもとで求めよ。

引用元リンク:

ときわ台学/解析学/ラグランジュの未定係数法

この例題では 目的関数と制約関数を明示的に与えてくれていますが、 図形的な意味でいうと前の例題と似ていて球体の中に収まる直方体の最大値を求める問題いうことになります。

In [1]:
%matplotlib notebook
In [2]:
import math
import numpy as np
import sympy as sy

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
In [3]:
sy.init_printing()
In [4]:
# イメージとしてはこんな感じ
fig = plt.figure()
ax = Axes3D(fig)

r = 1
u, v = np.mgrid[0:2 * np.pi:20j, 0:2 * np.pi:60j]

x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_wireframe(x, y, z, alpha=0.1, color='black')

corner = math.sqrt(3) / 3

X, Y = np.meshgrid([-corner, corner], [-corner, corner])
Z = np.ones(X.shape) * corner
ax.plot_surface(X, Y, Z)
ax.plot_surface(X, Y, -Z)
ax.plot_surface(Z, X, Y)
ax.plot_surface(-Z, X, Y)
ax.plot_surface(X, Z, Y)
ax.plot_surface(X, -Z, Y)

plt.show()
In [5]:
x, y, z, l = sy.symbols('x y z λ')
In [6]:
fx = 8 * x * y * z
fx
Out[6]:
$\displaystyle 8 x y z$
In [7]:
gx = x ** 2 + y ** 2 + z ** 2 - 1
gx
Out[7]:
$\displaystyle x^{2} + y^{2} + z^{2} - 1$
In [8]:
L = fx - l * gx
L
Out[8]:
$\displaystyle 8 x y z - λ \left(x^{2} + y^{2} + z^{2} - 1\right)$
In [9]:
dx = L.diff(x)
dx
Out[9]:
$\displaystyle - 2 x λ + 8 y z$
In [10]:
dy = L.diff(y)
dy
Out[10]:
$\displaystyle 8 x z - 2 y λ$
In [11]:
dz = L.diff(z)
dz
Out[11]:
$\displaystyle 8 x y - 2 z λ$
In [12]:
dl = L.diff(l)
dl
Out[12]:
$\displaystyle - x^{2} - y^{2} - z^{2} + 1$
In [13]:
points = sy.solve([dx, dy, dz, dl])
# 極小値は無視する
[p for p in points if p[l]]
Out[13]:
$\displaystyle \left[ \left\{ x : - \frac{\sqrt{3}}{3}, \ y : - \frac{\sqrt{3}}{3}, \ z : - \frac{\sqrt{3}}{3}, \ λ : - \frac{4 \sqrt{3}}{3}\right\}, \ \left\{ x : - \frac{\sqrt{3}}{3}, \ y : - \frac{\sqrt{3}}{3}, \ z : \frac{\sqrt{3}}{3}, \ λ : \frac{4 \sqrt{3}}{3}\right\}, \ \left\{ x : - \frac{\sqrt{3}}{3}, \ y : \frac{\sqrt{3}}{3}, \ z : - \frac{\sqrt{3}}{3}, \ λ : \frac{4 \sqrt{3}}{3}\right\}, \ \left\{ x : - \frac{\sqrt{3}}{3}, \ y : \frac{\sqrt{3}}{3}, \ z : \frac{\sqrt{3}}{3}, \ λ : - \frac{4 \sqrt{3}}{3}\right\}, \ \left\{ x : \frac{\sqrt{3}}{3}, \ y : - \frac{\sqrt{3}}{3}, \ z : - \frac{\sqrt{3}}{3}, \ λ : \frac{4 \sqrt{3}}{3}\right\}, \ \left\{ x : \frac{\sqrt{3}}{3}, \ y : - \frac{\sqrt{3}}{3}, \ z : \frac{\sqrt{3}}{3}, \ λ : - \frac{4 \sqrt{3}}{3}\right\}, \ \left\{ x : \frac{\sqrt{3}}{3}, \ y : \frac{\sqrt{3}}{3}, \ z : - \frac{\sqrt{3}}{3}, \ λ : - \frac{4 \sqrt{3}}{3}\right\}, \ \left\{ x : \frac{\sqrt{3}}{3}, \ y : \frac{\sqrt{3}}{3}, \ z : \frac{\sqrt{3}}{3}, \ λ : \frac{4 \sqrt{3}}{3}\right\}\right]$
In [14]:
# 体積
(2 * sy.sqrt(3) / 3) ** 3
Out[14]:
$\displaystyle \frac{8 \sqrt{3}}{9}$

直径1の円に内接する三角形の中で, その周の長さが最大のものを求めよ.

引用元リンク:

極値, 最大最小問題, ラグランジュの未定乗数法 - 名古屋大学

In [1]:
%matplotlib inline
In [2]:
import math
import numpy as np
import sympy as sy

from matplotlib import pyplot as plt, patches as pat
In [3]:
sy.init_printing()
In [4]:
# イメージとしてはこんな感じ
fig = plt.figure()
ax = plt.axes()
ax.set_aspect('equal')

r = 0.5
sq2 = math.sqrt(2)

# 円の表示
circle = pat.Circle(xy=(0, 0), radius=r, alpha=0.5)
ax.add_patch(circle)

# 三角形の表示
p1 = sy.Point(r, 0)
p2 = sy.Point(np.cos(math.radians(120)) * r, np.sin(math.radians(120)) * r)
p3 = sy.Point(np.cos(math.radians(240)) * r, np.sin(math.radians(240)) * r)
t = sy.Triangle(p1, p2, p3)
ax.add_patch(plt.Polygon(t.vertices))

# 軸の表示
plt.xlim(-r, r)
plt.ylim(-r, r)
plt.grid(True, color='black')
plt.tick_params(which='major', width=1)
plt.xticks([0])
plt.yticks([0])
plt.show()

内接する三角形の角をそれぞれ x, y, z と置いて、それと向かい合う辺を X, Y, Z とすると求めるべき三辺の合計は $$ X + Y + Z $$ となる。

ここで、正弦定理より (直径R) $$ \frac{X}{\sin x} = R $$ $$ \frac{Y}{\sin y} = R $$ $$ \frac{Z}{\sin z} = R $$ であり、先程の式を sin の式に置き換えると $$ R(\sin x + \sin y + \sin z) $$ となり、今回の問では 直径は 1 なので、Rは消え、目的関数 f は 以下のようになる $$ f(x, y, z) = \sin x + \sin y + \sin z $$

三角形の内角の和(π)より次のような関係となる $$ x + y + z = π $$

z を 上記の式を使って表すと $$ f(x, y, z) = \sin x + \sin y + \sin (π - x - y) $$ $$ f(x, y, z) = \sin x + \sin y + \sin (x + y) $$

この場合制約関数はないので、ラグランジュ関数は次のようになり $$ L(x, y, z) = f(x, y, z) $$

停留点を求めるために偏微分する $$ \frac{∂L}{∂x} = \frac{∂f}{∂x} = 0 $$ $$ \frac{∂L}{∂y} = \frac{∂f}{∂y} = 0 $$

偏微分した結果で連立方程式を作る $$ \begin{eqnarray} \left\{ \begin{array}{l} \cos x + \cos (x + y) = 0 \\ \cos y + \cos (x + y) = 0 \\ \end{array} \right. \end{eqnarray} $$ これを解いて $$ \cos x = \cos y $$ となる。

これらが等しくなる x, y の適切な組み合わせは x = y

yx を代入して因数分解 $$ cosx + cos(2x) = 0 \\ cosx + 2cos^2x - 1 = 0 \\ (2cosx - 1)(cosx + 1) = 0 \\ $$

$$ \cos x = -1, 1/2 $$$$ x = π, π/3 $$$$ (x, y) = (π, π), (π/3, π/3)$$

3辺の和は $$ \sin (π/3) + \sin (π/3) + \sin (π/3) = \frac{3 \sqrt{3}}{2} $$

なお、「面積」の最大を求める場合は 上記の結果に加えて ヘロンの公式を用いて解けます。 (ラグランジュの未定乗数法は使ってない

https://ja.wikipedia.org/wiki/%E3%83%98%E3%83%AD%E3%83%B3%E3%81%AE%E5%85%AC%E5%BC%8F

./heron_formula.png

ヘロンの公式(ヘロンのこうしき、英:Heron's formula)は任意の三角形の3辺a, b, c の長さから面積 S を求める公式。

In [2]:
%matplotlib inline
In [3]:
import math
import numpy as np
import sympy as sy

from matplotlib import pyplot as plt, patches as pat
In [4]:
sy.init_printing()

内接する三角形の辺をそれぞれ a, b, c とおいて、三辺の合計とすると $$ a + b + c = 2s $$

この s を使って 面積 S を以下のように表せる $$ S = \sqrt{s(s-a)(s-b)(s-c)} $$

続いて s-a, s-b, s-c を相加相乗平均で表すと

$$ \frac{(s-a) + (s-b) + (s-c)}{3} \geqq \sqrt[3]{(s-a)(s-b)(s-c)} $$

相乗平均の最大条件 は $$ s - a = s - b = s - c $$ となり、式を変形して a = b = c .

以下の理由より 相乗平均の最大条件三角形の面積の最大条件 は等しい

  • 三角形の性質より、 s-a, s-b, s-c はいずれも正の自然数
  • s = a + b + c は、正三角形(a = b = c) の時に最大となる

つまり、面積の最大条件も同様に a = b = c のときと言える

In [4]:
# (ヘロンの公式で解いてもいいけどなんとなく)
# 円に内接する正三角形の面積は中心で分割すると鈍角が 120° の二等辺三角形 x 3 と考えられる
fig = plt.figure()
ax = plt.axes()
ax.set_aspect('equal')

r = 0.5
sq2 = math.sqrt(2)

# 円の表示
circle = pat.Circle(xy=(0, 0), radius=r, alpha=0.5)
ax.add_patch(circle)

# 三角形の表示
p1 = sy.Point(r, 0)
p2 = sy.Point(np.cos(math.radians(120)) * r, np.sin(math.radians(120)) * r)
p3 = sy.Point(np.cos(math.radians(240)) * r, np.sin(math.radians(240)) * r)
t = sy.Triangle(p1, p2, p3)
ax.add_patch(plt.Polygon(t.vertices))

plt.plot([0, r], [0, 0], color='black')
plt.plot([0, np.cos(math.radians(120)) * r], [0, np.sin(math.radians(120)) * r], color='black')
plt.plot([0, np.cos(math.radians(240)) * r], [0, np.sin(math.radians(240)) * r], color='black')
plt.show()
In [5]:
# 二等辺三角形 * 3
r = sy.Symbol('r')
r * sy.sqrt(3) * r * sy.Rational(1, 2) * sy.Rational(1, 2) * 3
Out[5]:
$\displaystyle \frac{3 \sqrt{3} r^{2}}{4}$
In [6]:
# ヘロンの公式を使うやり方 (何故かrの平方根が取れない)
a = b = c = sy.sqrt(3) * r
s = (a + b + c) / 2
sy.sqrt(s * (s-a) * (s-b) * (s-c))
Out[6]:
$\displaystyle \frac{3 \sqrt{3} \sqrt{r^{4}}}{4}$

備考

なお、円に内接はせず、三辺の和が一定になるようないわゆる等周問題の場合、 s が定数になるため、円に内接する場合と同様に a = b = c が最大となります。

これについては 面積最大が正三角形であることの証明 - YouTube がわかりやすいです。

x^2 +y^2 = 1 の条件のもとで,関数f(x,y)=x^3+y^3−3x−3yの最大値最小値を求めよ.

引用元リンク:

極値, 最大最小問題, ラグランジュの未定乗数法 - 名古屋大学

In [1]:
%matplotlib notebook
In [2]:
import math
import numpy as np
import sympy as sy

from matplotlib import pyplot as plt, patches as pat
from mpl_toolkits.mplot3d import Axes3D
In [3]:
sy.init_printing()
In [4]:
def f(x, y):
    return x ** 3 + y ** 3 - x * 3 - y * 3
In [5]:
fig = plt.figure()
ax = Axes3D(fig)
# 目的関数の表示(青)
xs = np.arange(-2, 2, 0.1)
ys = np.arange(-2, 2, 0.1)
X, Y = np.meshgrid(xs, ys)
ax.plot_wireframe(X, Y, f(X, Y))

# 制約関数の表示(赤)
r = 1
u = np.linspace(0, np.pi * 2, 100)
xs = np.sin(u) * r
ys = np.cos(u) * r
zs = f(xs, ys)
ax.plot(xs, ys, zs, color="red", linewidth=3)