準備
長方形の面積最大化
たとえば2変数の問題として、長方形の周囲長
を一定として、その面積が最大となる長方形の形状と面積はどのようになるかを考える。この場合、長方形の辺の長さを
とすると、問題は以下のように表せる。
(1) ![Rendered by QuickLaTeX.com \begin{equation*} \mathrm{maximize} \quad S(x, y) = xy \quad \mathrm{subject~to} \quad x+y = \frac{L}{2} \end{equation*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-a0d4babc57496c47bebcb716760623a4_l3.png)
これは以下のように代数的に簡単に解けて、答えは正方形とわかる。
(2) ![Rendered by QuickLaTeX.com \begin{gather*} S = x \left( \frac{L}{2} - x \right) = -x^2 + \frac{L}{2} x = - \left( x - \frac{L}{4} \right)^2 + \frac{L^2}{16} \\ \max S = \frac{L^2}{16} \quad \mathrm{for} \: x = y = \frac{L}{4} \end{gather*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-5533cd2518c18c7b45fcbdff0797df5a_l3.png)
ただし変数の数が増えたり、目的関数や制約条件が複雑になると、解析的に解くのが面倒になる。
Lgrangeの未定乗数法による解
解法から先に示す。Lagrangeの未定乗数法では、目的関数
に対して以下の問題となる。
(3) ![Rendered by QuickLaTeX.com \begin{gather*} \mathrm{maximize} \quad L(x, y, \lambda) = S(x, y) - \lambda g(x, y) = xy - \lambda \left( x + y - \frac{L}{2} \right) \\ \mathrm{subject~to} \quad S(x, y) = xy, \quad g(x, y) = x + y - \frac{L}{2} \end{gather*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-ea52b32ee50110726a5742b3b9c0d69c_l3.png)
を最大化するために、
で偏微分した以下の方程式を設定する。
(4) ![Rendered by QuickLaTeX.com \begin{gather*} \frac{\partial L}{\partial x} = 0, \quad \frac{\partial L}{\partial y} = 0, \quad \frac{\partial L}{\partial \lambda} = 0 \end{gather*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-26948297a2820748dc6690d1aa79ea0e_l3.png)
これを計算すると
(5) ![Rendered by QuickLaTeX.com \begin{gather*} y - \lambda = 0, \quad x - \lambda = 0, \quad x + y - \frac{L}{2}= 0 \\ \therefore \lambda = x = y = \frac{L}{4} \end{gather*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-a1a486b134a5d294378ea1f64b49ba70_l3.png)
Lagrangeの未定乗数法の一般形
一般には、変数
について、目的関数
を制約条件
の下で最大化/最小化する問題として与えられる。
(6) ![Rendered by QuickLaTeX.com \begin{align*} & \mathrm{maximize/minimize} \quad f(\boldsymbol{x}) \\ & \mathrm{subject~to} \quad g(\boldsymbol{x}) = 0 \end{align*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-73a28e767ebf7e66ec1d1adcc650d859_l3.png)
この等式制約条件付き最大化/最小化問題は、以下のように
を導入して、連立方程式として表現される。
(7) ![Rendered by QuickLaTeX.com \begin{align*} & L(\boldsymbol{x}, \lambda) = f(\boldsymbol{x}) - \lambda g(\boldsymbol{x}) \\ & \frac{\partial L(\boldsymbol{x}, \lambda)}{\partial x_1} = \cdots = \frac{\partial L(\boldsymbol{x}, \lambda)}{\partial x_n} = \frac{\partial L(\boldsymbol{x}, \lambda)}{\partial \lambda} = 0 \end{align*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-2fea609295ffdb7ecdb99f9e3cdffff9_l3.png)
例題
例題1:平面と円
平面
について、制約条件
の下での極値を求める。
(8) ![Rendered by QuickLaTeX.com \begin{align*} & \mathrm{minimize} \quad f(x, y) = x + y - 1 \\ & \mathrm{subject~to} \quad x^2 + y^2 = 1 \end{align*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-88e51473db10fa6e829758700d25015f_l3.png)
lagrangeの未定乗数を導入して問題を定式化すると以下のようになる。
(9) ![Rendered by QuickLaTeX.com \begin{equation*} L(x, y, \lambda) = x + y - 1 - \lambda (x^2 + y^2 - 1) \end{equation*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-45009b22775327753693c7332b8d1e8a_l3.png)
(10) ![Rendered by QuickLaTeX.com \begin{align*} &\dfrac{\partial L}{\partial x} = 1 - 2 \lambda x = 0 \\ &\dfrac{\partial L}{\partial y} = 1 - 2 \lambda y = 0 \\ &\dfrac{\partial L}{\partial \lambda} = - x^2 - y^2 + 1 = 0 \end{align*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-828ebce0cd639676d9f39797a9951668_l3.png)
この連立方程式を解くと以下のようになり、解として2つの極値を得るが、それらは最大値と最小値に相当する。
(11) ![Rendered by QuickLaTeX.com \begin{gather*} x = y = \frac{1}{2 \lambda} \quad \Rightarrow \quad \frac{1}{4 \lambda ^2} + \frac{1}{4 \lambda ^2} = 1 \quad \Rightarrow \quad \lambda = \pm \frac{1}{\sqrt{2}}\\ \therefore \; x = y = \pm \frac{1}{\sqrt{2}} \approx \pm 0.7071\\ \max f(x, y) = \sqrt{2} - 1 \approx 0.414\\ \min f(x, y) = -\sqrt{2} - 1 \approx -2.414 \end{gather*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-27359643ef4eb3d7b9033559a4aa52af_l3.png)
これを目的関数のコンターと制約条件の線で表すと以下の通り。
![](http://taustation.com/wp1/wp-content/uploads/2020/02/lagrange_multiplier2.png)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
|
import numpy as np import matplotlib.pyplot as plt import matplotlib.patches as patches x = np.linspace(-1, 1, 20) y = np.linspace(-1, 1, 20) x, y = np.meshgrid(x, y) z = x + y - 1 fig = plt.figure() ax = fig.add_subplot(111) cntr = ax.contour(x, y, z ,levels=[-2.5, -2.414, -2, -1.5, -1, -0.5, 0, 0.414, 0.5]) ax.clabel(cntr) crc = patches.Circle(xy=(0, 0), radius=1, fill=False) ax.add_patch(crc) s2 = 1 / np.sqrt(2) ax.plot((0, s2), (s2, s2), c='k', linestyle='dotted') ax.plot((s2, s2), (0, s2), c='k', linestyle='dotted') ax.set_xlim((-1, 1)) ax.set_ylim((-1, 1)) ax.set_xticks([-1, -0.5, 0.5, 1]) ax.set_yticks([-1, -0.5, 0.5, 1]) ax.spines['bottom'].set_position('zero') ax.spines['left'].set_position('zero') ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.text(0.6, -0.08, "0.707", c='r') ax.text(-0.23, 0.68, "0.707", c='r') ax.set_aspect('equal') plt.show() |
例題2:凸関数と直線
下に凸な関数
について、直線
の制約条件下での最小値を求める。
(12) ![Rendered by QuickLaTeX.com \begin{align*} & \mathrm{minimize} \quad f(x, y) = x^2 + y^2 \\ & \mathrm {subject~to} \quad x + y - 1 = 0 \end{align*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-5f1c2824bbfc3e67bd1a640e21414492_l3.png)
ここでlagrangeの未定乗数を導入して問題を定式化すると以下のようになる。
(13) ![Rendered by QuickLaTeX.com \begin{equation*} L(x, y, \lambda) = x^2 + y^2 - \lambda (x + y - 1) \end{equation*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-a3012ebdada269c931e5d1917ab5583c_l3.png)
(14) ![Rendered by QuickLaTeX.com \begin{align*} &\dfrac{\partial L}{\partial x} = 2x - \lambda = 0 \\ &\dfrac{\partial L}{\partial y} = 2y - \lambda = 0 \\ &\dfrac{\partial L}{\partial \lambda} = - x - y + 1 = 0 \end{align*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-2d7f1f22bcb1e52cb9c23d8693a16d2f_l3.png)
この連立方程式を解くと以下のようになり、解は最小値1つとなる。
(15) ![Rendered by QuickLaTeX.com \begin{gather*} x = y = \frac{\lambda}{2} \quad \Rightarrow \quad \frac{\lambda}{2} + \frac{\lambda}{2} = 1 \quad \Rightarrow \quad \lambda = 1\\ \therefore \; x = y = \frac{1}{2} \quad , \quad \min f(x, y) = \frac{1}{2} \end{gather*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-7bf840984fb44c99614ee795acb9ff49_l3.png)
これを目的関数のコンターと制約条件の線で表すと以下の通り。
![](http://taustation.com/wp1/wp-content/uploads/2020/02/lagrange_multiplier1.png)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
|
import numpy as np import matplotlib.pyplot as plt x = np.linspace(-1.5, 1.5, 40) y = np.linspace(-1.5, 1.5, 40) x, y = np.meshgrid(x, y) z = x * x + y * y fig = plt.figure() ax = fig.add_subplot(111) cntr = ax.contour(x, y, z, levels=[0.5, 1, 1.5, 2, 2.5]) ax.clabel(cntr) ax.plot([-0.5, 1.5], [1.5, -0.5]) ax.plot((0, .5), (.5, .5), c='k', linestyle='dotted') ax.plot((.5, .5), (0, .5), c='k', linestyle='dotted') ax.set_xticks([-1, -0.5, 0.5, 1]) ax.set_yticks([-1, -0.5, 0.5, 1]) ax.spines['bottom'].set_position('zero') ax.spines['left'].set_position('zero') ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.set_aspect('equal') plt.show() |
なお、これを3次元で表示すると以下のようになる。青い曲面が目的関数で、赤い直線が制約条件となる。最適化問題は、制約条件を満たす曲面上の点(図中、赤い放物線)の最小値を求めることになる。
![](http://taustation.com/wp1/wp-content/uploads/2020/02/lagrange_multiplier1_3d.png)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
|
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D xmin, xmax = -2, 2 ymin, ymax = -2, 2 x = np.linspace(xmin, xmax) y = np.linspace(ymin, ymax) x_paraboloid, y_paraboloid = np.meshgrid(x, y) z_paraboloid = x_paraboloid**2 + y_paraboloid**2 x_constraint = x.copy() y_constraint = 1 - x z_constraint = x_constraint**2 + y_constraint**2 for i, v in enumerate(y_constraint): if v <= ymin or v >= ymax: x_constraint[i] = np.nan y_constraint[i] = np.nan z_constraint[i] = np.nan fig = plt.figure() ax = fig.add_subplot(111, projection='3d') ax.plot_surface(x_paraboloid, y_paraboloid, z_paraboloid, alpha=0.5) ax.plot(x_constraint, y_constraint, z_constraint, color='r') ax.plot(x_constraint, y_constraint, 0, color='r') plt.show() |
幾何学的説明
式(7)は以下のように書ける。
(16) ![Rendered by QuickLaTeX.com \begin{gather*} \left[ \begin{array}{c} \dfrac{\partial f}{\partial x_1} \\ \vdots \\ \dfrac{\partial f}{\partial x_n} \end{array} \right] = \lambda \left[ \begin{array}{c} \dfrac{\partial g}{\partial x_1} \\ \vdots \\ \dfrac{\partial g}{\partial x_n} \end{array} \right] \\ g(x_1, \ldots, x_n) = 0 \end{gather*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-c297beeeabd64856118393aad789fe1a_l3.png)
さらにgradientで表すと
(17) ![Rendered by QuickLaTeX.com \begin{gather*} \nabla f = \lambda \nabla g \\ g(x_1, \ldots, x_n) = 0 \end{gather*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-c103c62ee3c25b77c372f3f484f709af_l3.png)
すなわちこの式の解
は、制約条件である
を満足し、その曲線上にある。さらに解の点において
の勾配ベクトルとゼロ平面上における
の勾配ベクトルが平行になる。これはゼロ平面上の解の点において制約条件の曲線と
のコンターのうち特定の曲線が接するのと同義であり、この点は停留点(stationary point)である。
![](http://taustation.com/wp1/wp-content/uploads/2020/05/lagrange_mlutiplier_image.png)
つまりこのような停留点を発見する手順は、制約条件を満たす(制約条件の線上にある)点のうち、その点において目的関数のgradientと制約条件の関数のgradientが平行となる点を求めるということになる。再度これを式で表すと、
(18) ![Rendered by QuickLaTeX.com \begin{gather*} \nabla f(\boldsymbol{x}) = \lambda \nabla g(\boldsymbol{x}) \\ g(\boldsymbol{x}) = 0 \end{gather*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-8629d2ac5ba1b0ff692b8327077f12f4_l3.png)
となるが、これをLangrange関数
と定義したうえで各変数で偏微分したものをゼロと置いた方程式を解くと表現している。未定乗数λは停留点における目的関数のgradientと制約条件の関数のgradientの比を表している。
λの符号
λの符号に意味があるかどうか。
たとえば、以下の制約条件付き最適化問題を考える。
(19) ![Rendered by QuickLaTeX.com \begin{align*} & \mathrm{minimize} \quad f(x, y) = x^2 + y^2 \\ & \mathrm {subject~to} \quad x + y - 2 = 0 \end{align*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-2f39e4d483f4c197a565d6c9f0a0b14d_l3.png)
(20) ![Rendered by QuickLaTeX.com \begin{gather*} L(x, y, \lambda) = x^2 + y^2 - \lambda(x + y - 2) \\ \frac{\partial L}{\partial x} = \frac{\partial L}{\partial y} = \frac{\partial L}{\partial \lambda} = 0 \\ 2x - \lambda = 2y - \lambda = -x - y + 2 = 0\\ \lambda = 2, \; x=y=1 \end{gather*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-38035eee938d28e51f89d8f29444931d_l3.png)
この問題の制約条件を変更して以下のようにした場合。
(21) ![Rendered by QuickLaTeX.com \begin{align*} & \mathrm{minimize} \quad f(x, y) = x^2 + y^2 \\ & \mathrm {subject~to} \quad - x - y + 2 = 0 \end{align*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-7f4656467fee9a1fb0aec16f34884081_l3.png)
(22) ![Rendered by QuickLaTeX.com \begin{gather*} L(x, y, \lambda) = x^2 + y^2 - \lambda(- x - y + 2) \\ \frac{\partial L}{\partial x} = \frac{\partial L}{\partial y} = \frac{\partial L}{\partial \lambda} = 0 \\ 2x + \lambda = 2y + \lambda = x + y - 2 = 0\\ \lambda = -2, \; x=y=1 \end{gather*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-e7e2a08ec03c34bf24238792db3db254_l3.png)
このように、制約条件の正負を反転するとλの符号が逆になるが解は変わらない。これを表したのが以下の図。
等式制約条件の場合、制約条件の線上でgradientが平行になりさえすればよいので、λ符号(制約条件式の正負)には拘らなくてよい。
![](http://taustation.com/wp1/wp-content/uploads/2020/02/lagrange_multiplier_lambda_sign.png)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
|
import numpy as np import matplotlib.pyplot as plt import matplotlib.patches as patch from mpl_toolkits.mplot3d import Axes3D xmin, xmax = -3, 3 ymin, ymax = -3, 3 xls, xle = -1, 3 yls, yle = 3, -1 xs = np.linspace(xmin, xmax) ys = np.linspace(ymin, ymax) xm, ym = np.meshgrid(xs, ys) zero_plain = np.full_like(xm, 0) z_prb = xm**2 + ym**2 z_pl1 = xm + ym - 2 z_pl2 = - xm - ym + 2 r = np.sqrt(2) t = np.linspace(-np.pi, np.pi, num=100) xc = r * np.cos(t) yc = r * np.sin(t) uqc = 2 * xc vqc = 2 * yc xql = np.linspace(xls, xle, num=80) yql = np.linspace(yls, yle, num=80) uql1 = np.full_like(xql, 1) vql1 = np.full_like(yql, 1) uql2 = np.full_like(xql, -1) vql2 = np.full_like(yql, -1) fig = plt.figure(figsize=(11, 4.8)) ax1 = fig.add_subplot(121, projection='3d') ax2 = fig.add_subplot(122) ax1.plot_surface(xm, ym, zero_plain, color='gray', alpha=0.5) ax1.plot_surface(xm, ym, z_prb, color='tab:blue', alpha=0.5) ax1.plot_surface(xm, ym, z_pl1, color='tab:orange', alpha=0.5) ax1.plot_surface(xm, ym, z_pl2, color='tab:green', alpha=0.5) ax1.plot([xls, yls], [xle, yle], [0, 0], c='k') ax1.set_xlabel("x") ax1.set_ylabel("y") ax2.plot(xc, yc, c='b') ax2.plot([xls, yls], [xle, yle], c='k') ax2.quiver(xc, yc, uqc, vqc, scale_units='xy', scale=4, color='tab:blue', alpha=0.5) ax2.quiver(xql, yql, uql1, vql1, scale_units='xy', scale=4, color='tab:orange', alpha=0.5) ax2.quiver(xql, yql, uql2, vql2, scale_units='xy', scale=4, color='tab:green', alpha=0.5) ax2.set_xlabel("x") ax2.set_ylabel("y") ax2.set_xlim(0, 2) ax2.set_ylim(0, 2) plt.show() |
停留点が極致とならない例
以下の最適化問題の解をみてみる。
(23) ![Rendered by QuickLaTeX.com \begin{align*} & \mathrm{maxmize} \quad f(x, y) = x^3 + y^3 \\ & \mathrm{subject~to} \quad g(x, y) = x - y \end{align*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-a2d0a50fb8c6317553498ad7a386e902_l3.png)
(24) ![Rendered by QuickLaTeX.com \begin{align*} &L(x, y, \lambda) = x^3 + y^3 - \lambda(x - y) \\ &\left\{ \begin{array}{c} \dfrac{\partial L}{\partial x} = 3x^2 - \lambda = 0\\ \\ \dfrac{\partial L}{\partial y} = 3y^2 + \lambda= 0\\ \\ \dfrac{\partial L}{\partial \lambda} = x - y = 0 \end{array} \right. \end{align*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-106d5aedd6907280da45eec5230f6152_l3.png)
第1式から第2式を引き、第3式を適用して、
(25) ![Rendered by QuickLaTeX.com \begin{gather*} \3x^2 - 3y^2 - 2\lambda = 0 \\ 3(x + y)(x - y) = 2\lambda = 0 \\ \therefore x = y = \lambda = 0 \end{gather*}](http://taustation.com/wp1/wp-content/ql-cache/quicklatex.com-b1069eac54c2329585e75de39b988e17_l3.png)
ここで
の曲面と制約条件
に対する曲面上の軌跡を描くと下図左のようになる。下図の右は
として曲線を表したもので、
で勾配は水平になっており、停留点ではあるが極大/極小となっていない。
![](http://taustation.com/wp1/wp-content/uploads/2020/02/lagrange-multiplier-stationary-point.png)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
|
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D x = np.linspace(-1, 1) y = np.linspace(-1, 1) xs, ys = np.meshgrid(x, y) zs = xs**3 + ys**3 t = np.linspace(-1, 1) xp = yp = t zp = xp**3 + yp**3 fig = plt.figure(figsize=(13, 4.8)) ax1 = fig.add_subplot(121, projection='3d') ax2 = fig.add_subplot(122) ax1.plot_surface(xs, ys, zs, alpha=0.5) ax1.plot(xp, yp, zp) ax1.set_xlabel("x") ax1.set_ylabel("y") ax1.set_xticks(np.arange(-1, 1.5, 0.5)) ax1.set_yticks(np.arange(-1, 1.5, 0.5)) ax2.plot(t, zp) ax2.plot((-2, 2), (0, 0), c='k') ax2.plot((0, 0), (-2, 2), c='k') ax2.set_xlim(-2, 2) ax2.set_ylim(-2, 2) ax2.spines['top'].set_visible(False) ax2.spines['right'].set_visible(False) ax2.set_xlabel("t") ax2.set_ylabel("z") ax2.set_aspect('equal') ax2.grid(True) plt.show() |
参考サイト
本記事をまとめるにあたって、下記サイトが大変参考になったことに感謝したい。