点と直線の距離～媒介変数による愚直な解法

2020-06-21 / tau / コメントする

点Pと直線lの距離を、媒介変数を通して愚直に求める。具体的には点Pから直線lへの垂線の足Hと点Pの距離となる。

直点をパラメーター表示し、点H(x_H, y_H)に対応するパラメーターをt_Hとする。

(1) $\begin{equation*} \left\{ \begin{align} x_H &= u t_H + x_0\\ y_H &= v t_H + y_0 \end{align} \right. \end{equation*}$

このとき、直線に沿う方向のベクトルと直線に直角なベクトルを与えられた変数で表示すると以下の通り。

(2) $\begin{equation*} \boldsymbol{d} = \left[\begin{array}{c} u \\ v \end{array}\right] \quad , \quad \boldsymbol{h} = \left[\begin{array}{c} x_H - x_P \\ y_H - y_P \end{array}\right] \end{equation*}$

これらのベクトルが直交する条件として内積をゼロとする。

(3) $\begin{gather*} \boldsymbol{d} \cdot \boldsymbol{h} = 0 \\ u(x_H - x_P) + v(y_H - y_P) = 0 \\ u(u t_H + x_0 - x_P) + v(v t_H + y_0 - y_P) = 0 \\ (u^2 + v^2)t_H + ux_0 + vy_0 - ux_P - vy_P = 0 \end{gather*}$

これから点Hに対する直線のパラメーターが得られる。

(4) $\begin{equation*} t_H = \frac{u(x_P - x_0) + v(y_P - y_0)}{u^2 + v^2} \end{equation*}$

このパラメーターを使ってHの座標を表示すると以下の通り。

(5) $\begin{equation*} \left\{ \begin{align} x_H &= u t_H + x_0 = \frac{u(u x_P + v y_P) + v(v x_0 - u y_0)}{u^2 + v^2} \\ y_H &= v t_H + y_0 = \frac{v(u x_P + v y_P) - u(v x_0 - u y_0)}{u^2 + v^2} \\ \end{align} \right. \end{equation*}$

点Pと直線lの距離はベクトルhの大きさとなる。

(6) $\begin{equation*} |\boldsymbol{h}|^2 = (x_H - x_P)^2 + (y_H- y_P)^2 \end{equation*}$

以下、それぞれの項を計算していく。

(7) $\begin{equation*} \left\{ \begin{align} x_H - x_P &= \frac{v^2(x_0 - x_P) + uv(y_P-y_0)}{u^2 + v^2}\\ y_H - y_P &= \frac{u^2(y_0 - y_P) + uv(x_P-x_0)}{u^2 + v^2}\\ \end{align} \reigh. \end{equation*}$

(8) $\begin{equation*} \left\{ \begin{align} (x_H - x_P)^2 &= \frac{v^4(x_0 - x_P)^2 + u^2v^2(y_P-y_0)^2 + 2uv^3(x_0 - x_P)(y_P - y_0)}{(u^2 + v^2)^2}\\ (y_H - y_P)^2 &= \frac{u^4(y_0 - y_P)^2 + u^2v^2(x_P-x_0)^2 + 2u^3v(y_0 - y_P)(x_P - x_0)}{(u^2 + v^2)^2}\\ \end{align} \reigh. \end{equation*}$

これらの結果から以下を導く。

(9) $\begin{align*} |\boldsymbol{h}|^2 &= (x_H - x_P)^2 + (y_H - y_P)^2 \\ &= \frac{ \begin{array}{l} \phantom{+}v^2(u^2+v^2)(x_0 - x_P)^2 \\ + u^2(u^2+v^2)(y_0 - y_P)^2 \\ - 2uv(u^2+v^2)(x_0 - x_P)(y_0 - y_P) \end{array} }{(u^2 + v^2)^2} \\ &= \frac{v^2(x_0 - x_P)^2 + u^2(y_0 - y_P)^2 - 2uv(x_0 - x_P)(y_0 - y_P)} {u^2 + v^2} \\ &= \frac{[v(x_0 - x_P) - u(y_0 - y_P)]^2}{u^2 + v^2} \\ &= \frac{(v x_0 - u y_0 - v x_P + v y_P)^2}{u^2 + v^2} \end{align*}$

ここで直線の一般式とパラメーター表示を比較すると

(10) $\begin{gather*} \frac{x - x_0}{u} = \frac{y - y_0}{v} \\ vx - uy - v x_0 + u y_0 = 0 \\ ax + by + c = 0 \\ \therefore v = a , -u = b , -v x_0 + u y_0 = c \end{gather*}$

この結果を先の式に代入して以下を得る。

(11) $\begin{gather*} |\boldsymbol{h}|^2 = \frac{(a x_P + b x_P + c)^2}{a^2 + b^2} \\ |\boldsymbol{h}| = \frac{|a x_P + b x_P + c|}{\sqrt{a^2 + b^2}} \end{gather*}$

Lasso回帰の理解

2020-05-06 / tau / コメントする

定義

Ridge回帰は単純な多重回帰の損失関数に対してL2正則化項を加え、多重共線性に対する正則化を図った。Lasso解析はこれに対してL1正則化項を加えて最小化する(正則化の意味についてはこちら)。

(1) $\begin{align*} L &= \frac{1}{2} \sum_{i=1}^n ( y_i - \hat{y}_i )^2 + \alpha (|w_1| + \cdots + |w_m|) \\ &= \frac{1}{2} \sum_i ( y_i - w_0 - w_1 x_{1i} - \cdots - w_m x_{mi} )^2 + \alpha (|w_1| + \cdots + |w_m|) \end{align*}$

L1正則化の意味

準備

L2正則化は各重み係数が全体として小さくなるように制約がかかったが、L1正則化では値がゼロとなる重み係数が発生する。このことを確認する。

係数wを求めるためには損失関数Lを最小化すればよいが、Ridge回帰とは異なりL1正則化項は通常の解析的な微分はできない。

(2) $\begin{align*} \frac{\partial L}{\partial w_k} &= - \sum_i x_{ki} ( y_i - w_0 - w_1 x_{1i} - \cdots - w_m x_{mi} ) + \alpha \frac{\partial |w_k|}{\partial w_k} \\ &= - \sum_i x_{ki}y_i + w_0 \sum_i x_{ki} + \sum_{j \ne k} w_j \sum_i x_{ji} x_{ki} + w_k \sum_i {x_{ki}}^2 + \alpha \frac{\partial |w_k|}{\partial w_k} \\ &= 0 \end{align*}$

ここで $\frac{\partial |w_k|}{\partial w_k}=|w_k|'$ と表し、左辺のw_k以外に関わる項をM_k、w_kの係数となっている2乗和をS_kkと表す。

(3) $\begin{equation*} M_k + w_k S_{kk} + \alpha |w_k|' = 0 \end{equation*}$

場合分け

ここで|w_k|’についてはw_kの符号によって以下の値をとる。

(4) $\begin{equation*} |w_k|' = \left\{ \begin{array}{rl} -1 & (w_k < 0) \\ 1 & (w_k > 0) \end{array} \end{equation*}$

これらを式(3)に適用する。まずw_k < 0に対しては

(5) $\begin{gather*} w_k < 0 \quad \rightarrow \quad M_k + w_k S_{kk} - \alpha = 0 \\ -M_k + \alpha < 0 \quad \rightarrow \quad w_k = \frac{-M_k + \alpha}{S_{kk}} \end{gather*}$

またw_k > 0に対しては、

(6) $\begin{gather*} w_k > 0 \quad \rightarrow \quad M_k + w_k S_{kk} + \alpha = 0 \\ -M_k - \alpha > 0 \quad \rightarrow \quad w_k = \frac{-M_k - \alpha}{S_{kk}} \end{gather*}$

以上をまとめると、

(7) $\begin{equation*} w_k = \left\{ \begin{array}{ll} \dfrac{-M_k - \alpha}{S_{kk}} & (M_k < -\alpha) \\ \\ \dfrac{-M_k + \alpha}{S_{kk}} & (M_k > \alpha) \\ \end{array} \right. \end{equation*}$

劣微分の導入

式(7)で−α ≤ M_k ≤ αについては得られていない。M_k → ±αについてそれぞれの側から極限を計算すると0となるのでその間も0でよさそうだが、その保証はない。

ここでこちらのサイトのおかげで”劣微分(subdifferential)”という考え方を知ることができた。|w_k|’についてw_k = 0では解析的に微分不可能だが、その両側から極限をとった微分係数の範囲の集合を微分係数とするという考え方のようだ。

(8) $\begin{equation*} \frac{d |x|}{dx} = \left\{ \begin{array}{cl} -1 & (x < 0) \\ \left[ -1, 1 \right] & (x = 0) \\ 1 & (x > 0) \end{array} \right. \end{equation*}$

そこで、w_k = 0に対してこの劣微分を適用してみる。

(9) $\begin{gather*} w_k = 0 \quad \rightarrow \quad M_k +w_k S_{kk} + \alpha \left[ -1, 1 \right] = \left[ M_k - \alpha , M_k + \alpha \right] = 0\\ M_k - \alpha \le 0 \le M_k + \alpha \quad \rightarrow \quad -\alpha \le M_k \le \alpha \quad \rightarrow \quad w_k = 0 \end{gather*}$

以上のことから、重みw_kについて以下のようになり、−α≤M_k≤αの範囲ではw_k = 0となることがわかる。

(10) $\begin{equation*} w_k = \left\{ \begin{array}{cl} \dfrac{-M_k - \alpha}{S_{kk}} & (M_k < -\alpha \quad)\\ \\ 0 & (-\alpha \le M_k \le \alpha) \\ \\ \dfrac{-M_k + \alpha}{S_{kk}} & (M_k > \alpha) \end{array} \right. \end{equation*}$

すなわちL1正則化の場合、ハイパーパラメータαは重み係数の大きさを制限すると同時に重み係数がゼロとなるような効果も持ち、αが大きいほど多くの重み係数がゼロとなりやすい。

参考サイト

本記事をまとめるにあたって、下記サイトが大変参考になったことに感謝したい。

Lassoを数式から実装まで(理論編)～Miidas Research

ラグランジュの未定乗数法～不等式制約条件

2020-05-05 / tau / コメントする

概要

不等式制約条件を持つLagrangeの未定乗数法は以下のように表示される(等式条件の場合はこちらを参照)。

(1) $\begin{align*} & \mathrm{maximize/minimize} \quad f(\boldsymbol{x}) \\ & \mathrm{subject~to} \quad g(\boldsymbol{x}) \le 0 \end{align*}$

この場合、停留点が制約条件の範囲外にあれば等式条件と同じ問題となり、範囲内にあれば制約条件なしの通常の極値問題となる。

例題

停留点が制約条件の範囲内の場合

停留点が制約条件の境界内にある次のような問題を考える。

(2) $\begin{align*} & \mathrm{minimize} \quad f(x, y) = \left(x - \frac{1}{2} \right)^2 + \left(y - \frac{1}{2} \right)^2 +10 \\ & \mathrm{subject~to} \quad g(x, y) =x^2 + y^2 - 4\le 0 \\ &L(x, y, \lambda) = (x - 1)^2 + (y - 1)^2 + 10 - \lambda (x^2 + y^2 - 4) \end{align*}$

この問題を解くのに、まず停留点が制約条件の範囲内にあるかをチェックする。上の式で $\lambda=0$ と置いて最適化問題を解く。

(3) $\begin{gather*} \frac{\partial f}{\partial x} = 2\left(x - \frac{1}{2} \right) = 0 \\ \frac{\partial f}{\partial y} = 2\left(y - \frac{1}{2} \right) = 0 \\ \Downarrow \\ x = y = \frac{1}{2} \end{gather*}$

上記の解は制約条件 $g(x, y) \le 0$ を満足する。この停留点が極値だということが分かっていれば、これが問題の解ということになる。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def f(x, y, a, b):
    return (x - a)**2 + (y - b)**2 + 10

a, b = 1/2, 1/2
x = np.linspace(-3, 6)
y = np.linspace(-3, 6)
x, y = np.meshgrid(x, y)

t = np.linspace(-np.pi, np.pi)
xg = 2 * np.cos(t)
yg = 2 * np.sin(t)

x4min = 6 / np.sqrt(13)
y4min = 4 / np.sqrt(13)
fmin = f(x4min, y4min, a, b)
x4max = - 6 / np.sqrt(13)
y4max = - 4 / np.sqrt(13)
fmax = f(x4max, y4max, a, b)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.plot_surface(x, y, f(x, y, a, b), alpha=0.5)
ax.plot(xg, yg, 0, c='r')
ax.plot(xg, yg, f(xg, yg, a, b), c='r', linestyle='dashed')

ax.scatter(a, b, 10, c='b')
ax.scatter(a, b, 0, c='b')
ax.plot([a, a], [b, b], [0, 10], c='b', linestyle='dotted')

ax.plot([x4min, x4min], [y4min, y4min], [0, fmin], c='r', linestyle='dotted')
ax.plot([x4max, x4max], [y4max, y4max], [0, fmax], c='r', linestyle='dotted')

ax.set_xlabel("x")
ax.set_ylabel("y")

print("f_max={:8.4f}, for x={:7.4f}, y={:7.4f}".format(fmax, x4max, y4max))
print("f_min={:8.4f}, for x={:7.4f}, y={:7.4f}".format(fmin, x4min, y4min))

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

def f(x, y, a, b):

return (x - a)**2 + (y - b)**2 + 10

a, b = 1/2, 1/2

x = np.linspace(-3, 6)

y = np.linspace(-3, 6)

x, y = np.meshgrid(x, y)

t = np.linspace(-np.pi, np.pi)

xg = 2 * np.cos(t)

yg = 2 * np.sin(t)

x4min = 6 / np.sqrt(13)

y4min = 4 / np.sqrt(13)

fmin = f(x4min, y4min, a, b)

x4max = - 6 / np.sqrt(13)

y4max = - 4 / np.sqrt(13)

fmax = f(x4max, y4max, a, b)

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

ax.plot_surface(x, y, f(x, y, a, b), alpha=0.5)

ax.plot(xg, yg, 0, c='r')

ax.plot(xg, yg, f(xg, yg, a, b), c='r', linestyle='dashed')

ax.scatter(a, b, 10, c='b')

ax.scatter(a, b, 0, c='b')

ax.plot([a, a], [b, b], [0, 10], c='b', linestyle='dotted')

ax.plot([x4min, x4min], [y4min, y4min], [0, fmin], c='r', linestyle='dotted')

ax.plot([x4max, x4max], [y4max, y4max], [0, fmax], c='r', linestyle='dotted')

ax.set_xlabel("x")

ax.set_ylabel("y")

print("f_max={:8.4f}, for x={:7.4f}, y={:7.4f}".format(fmax, x4max, y4max))

print("f_min={:8.4f}, for x={:7.4f}, y={:7.4f}".format(fmin, x4min, y4min))

plt.show()

停留点が制約条件の範囲外の場合

今度は目的関数の停留点が制約条件の境界の範囲外にある問題を考えてみる。

(4) $\begin{align*} & \mathrm{minimize} \quad f(x, y) = (x - 3)^2 + (y-2)^2 +10 \\ & \mathrm{subject~to} \quad g(x, y) =x^2 + y^2 - 4\le 0 \\ &L(x, y, \lambda) = (x - 3)^2 + (y - 2)^2 + 10 - \lambda (x^2 + y^2 - 4) \end{align*}$

まず停留点が境界条件の範囲内にあるとして、 $\lambda = 0$ として問題を解いてみる。

(5) $\begin{gather*} \frac{\partial f}{\partial x} = 2(x - 3) = 0 \\ \frac{\partial f}{\partial y} = 2(y - 2) = 0 \\ \Downarrow \\ x = 3 , \; y = 2 \end{gather*}$

ところがこの解は制約条件 $g(x, y) \le 0$ を満足しない。解は制約条件の境界上にあることになるので、 $\lambda \ne 0$ の条件でLagrangeの未定乗数法によって解く。

(6) $\begin{gather*} \left\{ \begin{array}{c} \dfrac{\partial L}{\partial x} = 2(x - 3) - 2\lambda x = 0\\ \\ \dfrac{\partial L}{\partial y} = 2(y - 2) - 2\lambda y = 0\\ \\ \dfrac{\partial L}{\partial \lambda} = -(x^2 + y^2 - 4) = 0 \end{array} \right. \quad \Rightarrow \quad \left\{ \begin{array}{c} (1 - \lambda) x = 3\\ (1 - \lambda) y = 2\\ x^2 + y^2 = 4 \end{array} \right. \end{gather*}$

$\lambda \ne 1$ として、

(7) $\begin{gather*} x = \frac{3}{1 - \lambda} , \; y = \frac{2}{1 - \lambda} , \; 1 - \lambda = \pm \frac{\sqrt{13}}{2} \\ \left\{ \begin{array}{rll} \lambda &= 1 \mp \dfrac{\sqrt{13}}{2} & \fallingdotseq -0.8028 , \; 2.8028 \\ x &= \pm \dfrac{6}{\sqrt{13}} & \fallingdotseq 1.6641 , \; -1.6641 \\ y &= \pm \dfrac{4}{\sqrt{13}} & \fallingdotseq 1.1094 , \; -1.1094 \\ f(x, y) &= 17 \pm \dfrac{52}{\sqrt{13}} & \fallingdotseq 12.5778 , \; 41.4222 \end{array} \right. \end{gather*}$

いずれの解も制約条件 $g(x, y)=0$ を満たしており、このうち最小値となる解と目的関数の値は以下の通り。

(8) $\begin{gather*} (x, y, f(x, y)) = \left( \frac{6}{\sqrt{13}} , \frac{4}{\sqrt{13}} , 17 - \dfrac{52}{\sqrt{13}} \right) \fallingdotseq (1,6641, 1.1094, 12.5778) \end{gather*}$

λの符号について

有界でない例

以下、制約条件が $g(x, y) \le 0$ で有界でない場合を考える。

最小化問題でgradientが逆向き(極値が制約条件外)

まず、目的関数を最小化する場合を考える。以下の問題では目的関数fの極値が制約条件の範囲外にあり、fとgのgradientがgの境界上で逆向きになる。

(9) $\begin{align*} &\mathrm{minimize} \quad f(x, y) = (x - 1)^2 + (y - 1)^2 \\ &\mathrm{subject~to} \quad g(x, y) = x + y \le 0 \\ &L(x, y, z) = (x - 1)^2 + (y - 1)^2 - \lambda(x + y) \end{align*}$

停留点は制約条件の境界上となることから、λ≠0として未定乗数法で解いて以下を得る。λ<0となり、gradientが逆向きであることが現れている(∇f = λ∇g)。

(10) $\begin{equation*} x = 0, \; y = 0, \; \lambda = -2 \end{equation*}$

最小化問題でgradientが同じ向き(極値が制約条件内)

次の問題では極値が制約条件の範囲内にあり、fとgのgradientがgの境界上で同じ向きとなる。

(11) $\begin{align*} &\mathrm{minimize} \quad f(x, y) = (x + 1)^2 + (y + 1)^2 \\ &\mathrm{subject~to} \quad g(x, y) = x + y \le 0 \\ &L(x, y, z) = (x + 1)^2 + (y + 1)^2 - \lambda(x + y) \end{align*}$

この場合、fの極値は制約条件(緑色)の範囲内にあり、fとgのgradientがgの境界上で同じ向きになる。停留点は制約条件内にあることから、λ=0として制約条件の効果をなくし、単純にfの極値問題を解いて以下を得る(∇f =0)。

(12) $\begin{equation*} x = -1, \; y = -1, \; \lambda = 0 \end{equation*}$

最大化問題でgradientが同じ向き(極値が制約条件範囲外)

次に目的関数を最大化する場合を考える。以下の問題では目的関数fの極値が制約条件の範囲外にあり、fとgのgradientがgの境界上で同じ向きになる。

(13) $\begin{align*} &\mathrm{minimize} \quad f(x, y) = - (x - 1)^2 - (y - 1)^2 \\ &\mathrm{subject~to} \quad g(x, y) = x + y \le 0 \\ &L(x, y, z) = - (x - 1)^2 - (y - 1)^2 - \lambda(x + y) \end{align*}$

停留点は制約条件の境界上となることから、λ≠0として未定乗数法で解いて以下を得る。λ>0となり、gradientが同じ向きであることが現れている(∇f = λ∇g)。

(14) $\begin{equation*} x = 0, \; y = 0, \; \lambda = 2 \end{equation*}$

最大化問題でgradientが逆向き(極値が制約条件範囲内)

次の問題では極値が制約条件の範囲内にあり、fとgのgradientがgの境界上で逆向きとなる。

(15) $\begin{align*} &\mathrm{minimize} \quad f(x, y) = - (x + 1)^2 - (y + 1)^2 \\ &\mathrm{subject~to} \quad g(x, y) = x + y \le 0 \\ &L(x, y, z) = - (x + 1)^2 - (y + 1)^2 - \lambda(x + y) \end{align*}$

この場合、fの極値は制約条件(緑色)の範囲内にあり、fとgのgradientがgの境界上で逆向きになる。停留点は制約条件内にあることから、λ=0として制約条件の効果をなくし、単純にfの極値問題を解いて以下を得る(∇f =0)。

(16) $\begin{equation*} x = -1, \; y = -1, \; \lambda = 0 \end{equation*}$

λの符号について

以上の4つのパターンをまとめると、以下のようになる。

	∇fと∇gが同方向	∇fと∇gが逆方向
fを最小化	制約条件内 λ = 0, g(x, y) < 0	制約条件外(境界上) λ < 0, g(x, y) = 0
fを最大化	制約条件外(境界上) λ > 0, g(x, y) = 0	制約条件内 λ = 0, g(x, y) < 0

たとえばfを最大化する場合の条件をまとめて書くと、

(17) $\begin{equation*} g(x, y) \le 0 , \quad \lambda \ge 0 , \quad \lambda g(x, y) = 0 \end{equation*}$

これがKKT条件と呼ばれるもの。

ただし上記の場合分けは、最小化問題か最大化問題かによってλの符号が異なっている。さらに未定乗数を求める関数をここではL = f − λgとしたが、L = f + λgとしている場合もあり、このときはまたλの不等号が逆になってややこしい。

KKT条件では、一般にλ ≥ 0とされるが、元の問題が最小化か最大化か、不等式制約条件を正とするか負とするか、未定乗数法の関数形がf − λgかf + λgか、設定条件によってこの符号が反転する。逆に言えば、もともとのKKT条件がλ ≥ 0と示されているので、最小化(あるいは最大化)に対して−∇gとしたり、L = f ± gを使い分けたり、g ⋚ 0の設定を選んだりしている節がある。

KKT条件

不等式制約条件付きの最適化問題において、Lagrangeの未定乗数法を適用する場合の条件は、一般にKKT条件(Karush-Kuhn-Tucker condition)として示されるが、問題設定の形式を明示する必要がある(あるいは最小化問題と最大化問題でλの不等号を反転させた形で示す)。

最大化問題を基本とし、不等式制約条件は0又は負とする。ラグランジュ関数の制約条件項の符号はマイナスとする。

(18) $\begin{align*} &\mathrm{maximize} \quad f(\boldsymbol{x}) \\ &\mathrm{subject~to} \quad g(\boldsymbol{x}) \le 0 \\ &\rightarrow \; \mathrm{maximize} \quad L(\boldsymbol{x}, \lambda) = f(\boldsymbol{x}) - \lambda g(\boldsymbols{x}) \end{align*}$

この時、ラグランジュ関数を最大化するKKT条件は

(19) $\begin{gather*} g(x, y) \le 0 \\ \lambda \ge 0\\ \lambda g(x, y) = 0 \end{gather*}$

最小化問題の場合には、最大化問題となるように目的関数を反転する。

(20) $\begin{align*} &\mathrm{minimize} \quad f(\boldsymbol{x}) \\ &\mathrm{subject~to} \quad g(\boldsymbol{x}) \le 0 \\ &\rightarrow \; \mathrm{maximize} \quad L(\boldsymbol{x}, \lambda) = - f(\boldsymbol{x}) - \lambda g(\boldsymbols{x}) \end{align*}$

あるいは元のままの関数形として、最小化問題とKKT条件を以下のように設定してもよい。

(21) $\begin{gather*} \mathrm{minimize} \quad L(\boldsymbol{x}, \lambda) = f(\boldsymbol{x}) - \lambda g(\boldsymbols{x}) \\ g(x, y) \le 0 \\ \lambda \le 0\\ \lambda g(x, y) = 0 \end{gather*}$

参考サイト

本記事をまとめるにあたって、下記サイトが大変参考になったことに感謝したい。

arg max

2020-04-29 / tau / コメントする

maxが関数の最大値を意味するのに対して、arg maxは関数が最大値をとる場合の定義域を意味する。

(1) $\begin{align*} &\max x(4-x) = 4\\ &\arg \max x(4-x) = 2 \end{align*}$

定義域を指定する場合。

(2) $\begin{align*} &{\arg \max}_{-1 \le x \le 0} \, x(x + 1)(x - 1) = -\frac{1}{\sqrt{3}}\\ &\max_{-1 \le x \le 0} x(x + 1)(x - 1) = \frac{2\sqrt{3}}{9} \end{align*}$

本来arg maxは関数が最大値をとる定義域の集合を表す。

(3) $\begin{equation*} {\arg \max}_{0 \le x \le 4\pi} \, \cos x = \{ 0, 2\pi, 4\pi\} \end{equation*}$

行列式の定義

2020-04-23 / tau / コメントする

定義式

(1) $\begin{align*} |\boldsymbol{A}| &= \sum_{\sigma \in S_n} {\rm sgn}(\sigma) \prod_{i=1}^n a_{i\sigma(i)} \\ &= \sum_{\sigma \in S_n} {\rm sgn}(\sigma) a_{1\sigma(1)} a_{2\sigma(2)} \cdots a_{n\sigma(n)} \end{align*}$

計算例

次数2の場合

(2) $\begin{align*} |\boldsymbol{A}| &= \left| \begin{array}{cc} a_{11} & a_{12} \\ a_{21} & a_{22} \end{array} \right| \\ &= {\rm sgn}(\sigma_1)a_{1\sigma_1(1)}a_{2\sigma_1(2)} + {\rm sgn}(\sigma_2)a_{1\sigma_2(1)}a_{2\sigma_2(2)} \end{align*}$

ここで、

(3) $\begin{equation*} \sigma_1 = \left( \begin{array}{cc} 1 & 2 \\ 1 & 2 \end{array} \right) ,\quad \sigma_2 = \left( \begin{array}{cc} 1 & 2 \\ 2 & 1 \end{array} \right) \end{equation*}$

(4) $\begin{equation*} {\rm sgn}(\sigma_1) = 1 ,\quad {\rm sgn}(\sigma_2) = -1 \end{equation*}$

(5) $\begin{equation*} \sigma_1(1) = 1 ,\; \sigma_1(2) = 2 ,\; \sigma_2(1) = 2 ,\; \sigma_2(2) = 1 \end{equation*}$

したがって行列式の値は、

(6) $\begin{equation*} |\boldsymbol{A}| = a_{11}a_{22} - a_{12}a_{21} \end{equation*}$

次数3の場合

(7) $\begin{align*} |\boldsymbol{A}| =& \left| \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array} \right| \\ =& {\rm sgn}(\sigma_1)a_{1\sigma_1(1)}a_{2\sigma_1(2)}a_{3\sigma_1(3)} + {\rm sgn}(\sigma_2)a_{1\sigma_2(1)}a_{2\sigma_2(2)}a_{3\sigma_2(3)}\\ &{\rm sgn}(\sigma_3)a_{1\sigma_3(1)}a_{2\sigma_3(2)}a_{3\sigma_3(3)} + {\rm sgn}(\sigma_4)a_{1\sigma_4(1)}a_{2\sigma_4(2)}a_{4\sigma_4(3)}\\ &{\rm sgn}(\sigma_5)a_{1\sigma_5(1)}a_{2\sigma_5(2)}a_{3\sigma_5(3)} + {\rm sgn}(\sigma_6)a_{1\sigma_6(1)}a_{2\sigma_6(2)}a_{4\sigma_6(3)} \end{align*}$

ここで、

(8) $\begin{gather*} \sigma_1 = \left( \begin{array}{ccc} 1 & 2 & 3 \\ 1 & 2 & 3 \end{array} \right) ,\quad \sigma_2 = \left( \begin{array}{ccc} 1 & 2 & 3 \\ 1 & 3 & 2 \end{array} \right) \\ \sigma_3 = \left( \begin{array}{ccc} 1 & 2 & 3 \\ 2 & 1 & 3 \end{array} \right) ,\quad \sigma_4 = \left( \begin{array}{ccc} 1 & 2 & 3 \\ 2 & 3 & 1 \end{array} \right) \\ \sigma_5 = \left( \begin{array}{ccc} 1 & 2 & 3 \\ 3 & 1 & 2 \end{array} \right) ,\quad \sigma_6 = \left( \begin{array}{ccc} 1 & 2 & 3 \\ 3 & 2 & 1 \end{array} \right) \end{gather*}$

(9) $\begin{equation*} \begin{array}{ll} {\rm sgn}(\sigma_1) = \phantom{-}1 ,& {\rm sgn}(\sigma_2) = -1 \\ {\rm sgn}(\sigma_3) = -1 ,& {\rm sgn}(\sigma_4) = \phantom{-}1 \\ {\rm sgn}(\sigma_5) = \phantom{-}1 ,& {\rm sgn}(\sigma_6) = -1 \end{array} \end{equation*}$

(10) $\begin{gather*} \sigma_1(1) = 1 ,\; \sigma_1(2) = 2 ,\; \sigma_1(3) = 3 \\ \sigma_2(1) = 1 ,\; \sigma_2(2) = 3 ,\; \sigma_2(3) = 2 \\ \sigma_3(1) = 2 ,\; \sigma_3(2) = 1 ,\; \sigma_3(3) = 3 \\ \sigma_4(1) = 2 ,\; \sigma_2(2) = 3 ,\; \sigma_4(3) = 1 \\ \sigma_5(1) = 3 ,\; \sigma_5(2) = 1 ,\; \sigma_5(3) = 2 \\ \sigma_6(1) = 3 ,\; \sigma_6(2) = 2 ,\; \sigma_6(3) = 1 \end{gather*}$

したがって行列式の値は、

(11) $\begin{align*} |\boldsymbol{A}| &= a_{11}a_{22}a_{33} - a_{11}a_{23}a_{32} \\ &- a_{12}a_{21}a_{33} + a_{12}a_{23}a_{31} \\ &+ a_{13}a_{21}a_{32} - a_{13}a_{22}a_{31} \end{align*}$

ベクトル・行列を含む微分

2020-04-19 / tau / 20件のコメント

記号の定義

以下の記号で統一的に定義しておく。ベクトルは原則として列ベクトル表示を標準とする。

(1) $\begin{equation*} \boldsymbol{x} = \left[ \begin{array}{c} x_1 \\ \vdots \\ x_n \\ \end{array} \right] \end{equation*}$

(2) $\begin{equation*} \boldsymbol{X} = \left[ \begin{array}{ccc} x_{11} & \ldots & x_{1n} \\ \vdots & x_{ij} & \vdots \\ x_{m1} & \ldots & x_{mn} \end{array} \right] \end{equation*}$

(3) $\begin{equation*} f(\boldsymbol{x}) = f(x_1, \ldots, x_m) \end{equation*}$

(4) $\begin{equation*} \boldsymbol{f}(x) = \left[ \begin{array}{c} f_1(x) \\ \vdots \\ f_n(x) \end{array} \right] \end{equation*}$

(5) $\begin{equation*} \boldsymbol{f}(\boldsymbol{x}) =\ \left[ \begin{array}{c} f_1(x_1, \ldots, x_n) \\ \vdots \\ f_m(x_1, \ldots, x_n) \end{array} \right] \end{equation*}$

(6) $\begin{equation*} \boldsymbol{F}(x) = \left[ \begin{array}{ccc} F_{11}(x) & \ldots & F_{1n}(x) \\ \vdots & F_{ij} & \vdots \\ F_{m1}(x) & \ldots & F_{mn}(x) \\ \end{array} \right] \end{equation*}$

ベクトル・行列をスカラーで微分

これらは素直にベクトル・行列の要素を微分すればよい。

(7) $\begin{equation*} \frac{d \boldsymbol{f}(x)}{dx} = \left[ \begin{array}{c} \dfrac{d f_1(x)}{dx} \\ \vdots \\ \dfrac{d f_n(x)}{dx} \end{array} \right] \end{equation*}$

(8) $\begin{equation*} \frac{d \boldsymbol{F}(x)}{dx} = \left[ \begin{array}{ccc} \dfrac{dF_{11}(x)}{dx} & \ldots & \dfrac{dF_{1n}(x)}{dx} \\ \vdots & \dfrac{dF_{ij}(x)}{dx} & \vdots \\ \dfrac{dF_{m1}(x)}{dx} & \ldots & \dfrac{dF_{mn}(x)}{dx} \end{array} \right] \end{equation*}$

スカラーをベクトルで微分

スカラーを $\mathbb{R}^n$ のベクトルで微分すると、同じ次数のベクトルになる。

(9) $\begin{equation*} \frac{df(\boldsymbol{x})}{d\boldsymbol{x}} = \left[ \begin{array}{c} \dfrac{\partial f}{\partial x_1} \\ \vdots \\ \dfrac{\partial f}{\partial x_m} \end{array} \right] \end{equation*}$

これは便宜的に偏微分係数を要素とするベクトルを導入して以下のように考えるとよい。

(10) $\begin{equation*} \frac{d}{d\boldsymbol{x}}f(\boldsymbol{x}) = \left[ \begin{array}{c} \dfrac{\partial}{\partial x_1} \\ \vdots \\ \dfrac{\partial}{\partial x_m} \end{array}\right] f(\boldsymbol{x}) \end{equation*}$

スカラーを行列で微分

スカラーを $\mathbb{R}^m\times\mathbb{R}^n$ の行列で微分すると、同じ次元・次数の行列になる。

(11) $\begin{equation*} \frac{df(\boldsymbol{X})}{d\boldsymbol{X}} = \left[ \begin{array}{ccc} \dfrac{\partial f}{\partial x_{11}} & \ldots & \dfrac{\partial f}{\partial x_{1n}} \\ \vdots & \dfrac{\partial f}{\partial x_{ij}} & \vdots \\ \dfrac{\partial f}{\partial x_{m1}} & \ldots & \dfrac{\partial f}{\partial x_{mn}} \\ \end{array} \right] \end{equation*}$

これは便宜的に以下のように考えるとよい。

(12) $\begin{equation*} \frac{d}{d\boldsymbol{X}} f(\boldsymbol{X}) = \\ \left[ \begin{array}{ccc} \dfrac{\partial}{\partial x_{11}} & \ldots & \dfrac{\partial}{\partial x_{1n}} \\ \vdots & \dfrac{\partial}{\partial x_{ij}} & \vdots \\ \dfrac{\partial}{\partial x_{m1}} & \ldots & \dfrac{\partial}{\partial x_{mn}} \\ \end{array} \right] f(\boldsymbol{X}) \end{equation*}$

ベクトルをベクトルで微分

この場合、微分する変数側を行ベクトルとするか、微分される関数側を行ベクトルとするか2通りの表現があるが、ここでは関数側を行ベクトルとする。

(13) $\begin{equation*} \frac{d\boldsymbol{f}(\boldsymbol{x})^T}{d\boldsymbol{x}} = \left[ \begin{array}{ccc} \dfrac{\partial f_1}{\partial x_1} & \ldots & \dfrac{\partial f_n}{\partial x_1} \\ \vdots & \dfrac{\partial f_j}{\partial x_i} & \vdots \\ \dfrac{\partial f_1}{\partial x_m} & \ldots & \dfrac{\partial f_n}{\partial x_m} \\ \end{array} \right] \end{equation*}$

これは(10)で導入した偏微分係数ベクトルを導入して、便宜的に以下のように考えるとよい。

(14) $\begin{equation*} \frac{d}{d\boldsymbol{x}} \boldsymbol{f}(\boldsymbol{x})^T = \left[ \begin{array}{c} \dfrac{\partial}{\partial x_1} \\ \vdots \\ \dfrac{\parial}{\partial x_m} \end{array} \right] [ f_1(\boldsymbol{x}) \; \ldots \; f_n(\boldsymbol{x}) ] = \left[ \begin{array}{ccc} \dfrac{\partial f_1}{\partial x_1} & \cdots & \dfrac{\partial f_n}{\partial x_1} \\ & \dfrac{\partial f_j}{\partial x_i} &\\ \dfrac{\partial f_1}{\partial x_m} & \cdots & \dfrac{\partial f_n}{\partial x_m} \end{array} \right] \end{equation*}$

公式

一般形

単位行列

ベクトルを同じベクトルで微分すると、単位ベクトルではなく単位行列になる。

(15) $\begin{equation*} \frac{d\boldsymbol{x}}{d\boldsymbol{x}} = \boldsymbol{I} \end{equation*}$

合成関数

スカラーの合成関数と似ているが、イメージと積の順番が逆で、この順番は変えられない。

(16) $\begin{align*} \frac{df(\boldsymbol{u}(\boldsymbol{x}))}{d\boldsymbol{x}} = \frac{d\boldsymbol{u}(\boldsymbol{x})^T}{d\boldsymbol{x}} \frac{df(\boldsymbol{u})}{d\boldsymbol{u}} \\ \end{align*}$

これは以下のように確認できる。

(17) $\begin{align*} \frac{df}{dx_i} &= \frac{\partial f}{\partial u_1}\frac{\partial u_1}{\partial x_i} + \cdots + \frac{\partial f}{\partial u_j}\frac{\partial u_j}{\partial x_i} + \cdots + \frac{\partial f}{\partial u_n}\frac{\partial u_n}{\partial x_i} \\ \rightarrow \frac{df(\boldsymbol{u}(\boldsymbol{x}))}{d\boldsymbol{x}} &= \left[ \begin{array}{ccc} \dfrac{\partial u_1}{\partial x_1} & \cdots & \dfrac{\partial u_n}{\partial x_1} \\ \vdots && \vdots \\ \dfrac{\partial u_1}{\partial x_m} & \cdots & \dfrac{\partial u_n}{\partial x_m} \\ \end{array} \right] \left[ \begin{array}{c} \dfrac{\partial f}{\partial u_1} \\ \vdots \\ \dfrac{\partial f}{\partial u_n} \end{array} \right] \\ &= \left[ \begin{array}{c} \dfrac{\partial}{\partial x_1} \\ \vdots \\ \dfrac{\partial}{\partial x_m} \end{array} \right] [u_1 \; \cdots \; u_n] \left[ \begin{array}{c} \dfrac{\partial}{\partial u_1} \\ \vdots \\ \dfrac{\partial}{\partial u_n} \end{array} \right] f(\boldsymbol{u}) \end{align*}$

積の微分

行列の積のスカラーによる微分

(18) $\begin{equation*} \frac{d\boldsymbol{(FG)}}{dx} = \frac{d\boldsymbol{F}}{dx}\boldsymbol{G} + \boldsymbol{F} \frac{d\boldsymbol{G}}{dx} \end{equation*}$

これは素直に次のように確認できる。

(19) $\begin{align*} \frac{d(\boldsymbol{FG})}{dx} &= \frac{d}{dx}\left[\sum_{j=1}^m f_{ij} g_{jk}\right] = \left[\sum_{j=1}^m \left(\frac{df_{ij}}{dx} g_{jk} + f_{ij} \frac{dg_{jk}}{dx} \right)\right] \\ &= \left[ \sum_{j=1}^m \frac{df_{ij}}{dx} g_{jk} \right] + \left[ \sum_{j=1}^m f_{ij} \frac{dg_{jk}}{dx} \right] \end{align*}$

一次～二次形式の微分

Axの形式

ベクトルAxをベクトルxで微分するのに(13)の考え方で計算する。

ベクトルxをn次、行列Aをm×nとすると、ベクトルAx次数はm次となる。微分する際に微分係数ベクトルはn×1の列ベクトル、Axは転置して1×mの行ベクトルとなり、結果はn×mの行列になる。その結果は元のAの転置行列となる。

(20) $\begin{equation*} \frac{d\left( (\boldsymbol{Ax})^T \right)}{d\boldsymbol{x}} = \boldsymbol{A}^T \end{equation*}$

(21) $\begin{equation*} \begin{align} & \frac{d}{d\boldsymbol{x}} \left[ \begin{array}{c} a_{11}x_1 + \cdots + a_{1j}x_j + \cdots + a_{1n}x_n \\ \vdots \\ a_{i1}x_1 + \cdots + a_{ij}x_j + \cdots + a_{in}x_n \\ \vdots \\ a_{m1}x_1 + \cdots + a_{mj}x_j + \cdots + a_{mn}x_n \\ \end{array} \right]^T \\ &= \left[ \begin{array}{c} \dfrac{\partial}{\partial x_1} \\ \vdots \\ \dfrac{\partial}{\partial x_j} \\ \vdots \\ \dfrac{\partial}{\partial x_n} \\ \end{array} \right] \left[ \begin{array}{c} a_{11}x_1 + \cdots + a_{1j}x_j + \cdots + a_{1n}x_n \\ \vdots \\ a_{i1}x_1 + \cdots + a_{ij}x_j + \cdots + a_{in}x_n \\ \vdots \\ a_{m1}x_1 + \cdots + a_{mj}x_j + \cdots + a_{mn}x_n \\ \end{array} \right]^T \\ &= \left[ \begin{array}{ccccc} a_{11} & \cdots & a_{i1} & \cdots & a_{m1} \\ \vdots && \vdots && \vdots \\ a_{1j} & \cdots & a_{ij} & \cdots & a_{mj} \\ \vdots && \vdots && \vdots \\ a_{1n} & \cdots & a_{in} & \cdots & a_{mn} \\ \end{array} \right] = \boldsymbol{A}^T \end{align} \end{equation*}$

x^2の形式

(22) $\begin{equation*} \frac{d(\boldsymbol{x}^T \boldsymbol{x})}{d\boldsymbol{x}} = 2 \boldsymbol{x} \end{equation*}$

［証明］

(23) $\begin{equation*} \left[ \begin{array}{c} \dfrac{\partial}{\partial x_1} \\ \vdots \\ \dfrac{\partial}{\partial x_n} \end{array} \right] [ x_1^2 + \cdots + x_n^2 ] = \left[ \begin{array}{c} 2x_1 \\ \vdots \\ 2x_n \end{array} \right] \end{equation*}$

x^TAxの形式

この場合、 $\boldsymbol{A}$ は正方行列で、 $\boldsymbol{x}$ と同じ次数でなければならない。

(24) $\begin{equation*} \frac{d}{d\boldsymbol{x}} \left(\boldsymbol{x}^T \boldsmbol{A} \boldsymbol{x} \right)= \left( \boldsymbol{A} + \boldsymbol{A}^T \right) \boldsymbol{x} \end{equation*}$

［証明］

(25) $\begin{align*} &\frac{d}{d \boldsymbol{x}} \left( [x_1 \; \cdots \; x_n] \left[ \begin{array}{ccc} a_{11} & \cdots & a_{1n} \\ \vdots & & \vdots \\ a_{n1} & \cdots & a_{nn} \\ \end{array} \right] \left[ \begin{array}{c} x_1 \\ \vdots \\ x_n \end{array} \right] \right)\\ &=\frac{d}{d \boldsymbol{x}} \left( [x_1 \; \cdots \; x_n] \left[ \begin{array}{c} a_{11} x_1 + \cdots + a_{1n} x_n \\ \vdots \\ a_{n1} x_1 + \cdots + a_{nn} x_n \end{array} \right] \right)\\ &= \frac{d}{d \boldsymbol{x}} \left( \left( a_{11} {x_1}^2 + \cdots + a_{1n} x_1 x_n \right) + \cdots + \left( a_{n1} x_n x_1 + \cdots + a_{nn} {x_n}^2 \right) \right) \end{array}\\ &=\left[ \begin{array}{c} \left( 2 a_{11} x_1 + \cdots + a_{1n} x_n \right) + a_{21} x_2 + \cdots + a_{n1} x_n \\ \vdots \\ a_{1n} x_1 + \cdots + a_{1n-1} x_{n-1} + \left( a_{n1} x_1 + \cdots + 2a_{nn} x_n \right) \end{array} \right] \\ &=\left[ \begin{array}{c} \left( a_{11} x_1 + \cdots + a_{1n} x_1 \right) + \left( a_{11} x_1 + \cdots + a_{n1} x_n \right) \\ \vdots \\ \left( a_{n1} x_1 + \cdots + a_{nn} x_n \right) + \left( a_{1n} x_1 + \cdots + a_{nn} x_n \right) \end{array} \right] \end{align*}$

転置行列

2020-04-18 / tau / コメントする

定義

(1) $\begin{equation*} {\boldsymbol{A}^T}_{ij} = {\boldsymbol{A}}_{ji} \end{equation*}$

性質

単独の行列

転置の転置

(2) $\begin{equation*} \left(\boldsymbol{A}^T\right)^T = \boldsymbol{A} \end{equation*}$

逆行列

(3) $\begin{equation*} \left( \boldsymbol{A}^T \right)^{-1} = \left( \boldsymbol{A}^{-1} \right)^T \end{equation*}$

［証明］

(4) $\begin{align*} &\boldsymbol{AA}^{-1} = \boldsymbol{I} \; \Leftrightarrow \; \left( \boldsymbol{AA}^_{-1} \right)^T = \boldsymbol{I}^T \; \Leftrightarrow \; \left( \boldsymbol{A}^{-1} \right)^T \boldsymbol{A}^T = \boldsymbol{I} \\ &\boldsymbol{A}^{-1} \boldsymbol{A} = \boldsymbol{I} \; \Leftrightarrow \quad \left( \boldsymbol{A}^{-1} \boldsymbol{A} \right)^T = \boldsymbol{I}^T \; \Leftrightarrow \; \boldsymbol{A}^T \left( \boldsymbol{A}^{-1} \right)^T = \boldsymbol{I} \end{align*}$

行列式

(5) $\begin{equation*} \left| \boldsymbol{A}^T \right| = \left| \boldsymbol{A} \right| \end{equation*}$

行列演算

線形性

(6) $\begin{equation*} (\alpha \boldsymbol{A})^T = \alpha \boldsymbol{A}^T \end{equation*}$

(7) $\begin{equation*} \left(\boldsymbol{A} + \boldsymbol{B}\right)^T = \boldsymbol{A}^T + \boldsymbol{B}^T \end{equation*}$

積

交換法則は成り立たない。

(8) $\begin{equation*} (\boldsymbol{AB})^T = \boldsymbol{B}^T \boldsymbol{A}^T \end{equation*}$

［証明］

(9) $\begin{align*} \left( [ \boldsymbol{AB} ]_{ij} \right)^T = \left( \sum_k \boldsymbol{A}_{ik} \boldsymbol{B}_{kj} \right)^T = \sum_k \boldsymbol{B}_{jk} \boldsymbol{A}_{ki} =\boldsymbol{B}^T \boldsymbol{A}^T \end{align*}$

行列とベクトル

行列とベクトルの積

(10) $\begin{equation*} (\boldsymbol{Ax})^T = \boldsymbol{x}^T \boldsymbol{A}^T \end{equation*}$

すなわち

(11) $\begin{equation*} \left[ \begin{array}{ccc} a_{11} & \cdots & a_{1n} \\ \vdots & & \vdots \\ a_{m1} & \cdots & a_{mn} \\ \end{array} \right] \left[ \begin{array}{c} x_1 \\ \vdots \\ x_n \end{array} \right] = [x_1 \; \cdots \; x_n] \left[ \begin{array}{ccc} a_{11} & \cdots & a_{m1} \\ \vdots & & \vdots \\ a_{1n} & \cdots & a_{mn} \\ \end{array} \right] \end{equation*}$

内積

(12) $\begin{equation*} \boldsymbol{x}^T \boldsymbol{x} = \langle \boldsymbol{x} , \boldsymbol{x} \rangle \end{equation*}$

すなわち

(13) $\begin{equation*} [x_1 \; \cdots \; x_n] \left[ \begin{array}{c} x_1 \\ \vdots \\ x_n \end{array} \right] = x_1^2 + \cdots + x_n^2 \end{equation*}$

Gradient～勾配ベクトル

2020-02-27 / tau / コメントする

定義

多変数の関数のグラジエント(gradient)は勾配とも呼ばれ、以下で定義される。

(1) $\begin{equation*} \nabla f( x_1, \ldots, x_n ) = \left( \frac{\partial f}{\partial x_1} , \ldots , \frac{\partial f}{\partial x_n} \right) = \left( f_{x_1} , \ldots , f_{x_n} \right) \end{equation*}$

2変数の場合

(2) $\begin{equation*} \nabla f( x, y ) = \left( \frac{\partial f}{\partial x} , \ldots , \frac{\partial f}{\partial y} \right) = \left( f_x , \ldots , f_y \right) \end{equation*}$

準備～1変数関数の微分係数

微分係数の符号だけを見た場合

gradientの前に、1変数関数の微分係数についてもう一度考えてみる。

たとえば $y = (x - 1)^2 + 1$ は $x = 1$ で極値を持ち、その前後で微分係数の符号が変わる。

(3) $\begin{equation*} y' = 2(x - 1) \left\{ \begin{align} < 0 \quad (x < 1) \\ = 0 \quad (x = 1) \\ > 0 \quad (x > 1) \end{align} \right. \end{equation*}$

教科書的には、極値をとる点の前後で減少から増加に変わるということになる。

import numpy as np
import matplotlib.pyplot as plt


x_min = -1
x_max = 3

x = np.linspace(x_min, x_max, 50)
y = (x - 1)**2 + 1

n_points = 21
u0 = np.linspace(x_min, x_max, n_points)
u = 2 * (u0 - 1)
zeros = np.zeros(n_points)

fig, ax = plt.subplots()

ax.plot(x, y)
ax.quiver(u0, zeros, zeros, u)

ax.set_xlim(x_min, x_max)
ax.set_ylim(x_min, x_max)

ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')

ax.set_aspect('equal')

plt.show()

import numpy as np

import matplotlib.pyplot as plt

x_min = -1

x_max = 3

x = np.linspace(x_min, x_max, 50)

y = (x - 1)**2 + 1

n_points = 21

u0 = np.linspace(x_min, x_max, n_points)

u = 2 * (u0 - 1)

zeros = np.zeros(n_points)

fig, ax = plt.subplots()

ax.plot(x, y)

ax.quiver(u0, zeros, zeros, u)

ax.set_xlim(x_min, x_max)

ax.set_ylim(x_min, x_max)

ax.spines['top'].set_visible(False)

ax.spines['right'].set_visible(False)

ax.spines['bottom'].set_position('zero')

ax.spines['left'].set_position('zero')

ax.set_aspect('equal')

plt.show()

微分係数の符号を方向としてみた場合

同じ微分係数の値を、それが得られる値に対して、符号を考慮して矢印で表してみる。

ただしここでは、矢印の重なりを避けるために、少しずつ縦方向にずらして描いている。

こう描いてみると分かるが、微分係数の正負が方向を表すと考えると、「微分係数は、その方向に関数がどれだけ増加の傾きを持っているか」という意味を持つことがわかる。

import numpy as np
import matplotlib.pyplot as plt


x_min = -1
x_max = 3

x = np.linspace(x_min, x_max, 50)
y = (x - 1)**2 + 1

n_points = 21
u0 = np.linspace(x_min, x_max, n_points)
u = 2 * (u0 - 1)
v0 = np.linspace(-0.5, 0.5, n_points)
zeros = np.zeros(n_points)

fig, ax = plt.subplots()

ax.plot(x, y)
ax.quiver(u0, v0, u, zeros)

ax.set_xlim(x_min, x_max)
ax.set_ylim(x_min, x_max)

ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')

ax.set_aspect('equal')

plt.show()

import numpy as np

import matplotlib.pyplot as plt

x_min = -1

x_max = 3

x = np.linspace(x_min, x_max, 50)

y = (x - 1)**2 + 1

n_points = 21

u0 = np.linspace(x_min, x_max, n_points)

u = 2 * (u0 - 1)

v0 = np.linspace(-0.5, 0.5, n_points)

zeros = np.zeros(n_points)

fig, ax = plt.subplots()

ax.plot(x, y)

ax.quiver(u0, v0, u, zeros)

ax.set_xlim(x_min, x_max)

ax.set_ylim(x_min, x_max)

ax.spines['top'].set_visible(False)

ax.spines['right'].set_visible(False)

ax.spines['bottom'].set_position('zero')

ax.spines['left'].set_position('zero')

ax.set_aspect('equal')

plt.show()

微分係数の曲線上での分布

関数の曲線上で、微分係数の方向と大きさを考慮してベクトルを描いてみる。ベクトルの水平方向の成分を微分係数と同じ符号で長さ1、垂直方向の成分を微分係数の絶対値としている。

こうしてみると、微分係数は各位置での関数の増加方向とその増分を表していることがわかる。

import numpy as np
import matplotlib.pyplot as plt


x_min = -1
x_max = 3

x = np.linspace(x_min, x_max, 50)
y = (x - 1)**2 + 1

n_points = 21
u0 = np.linspace(x_min, x_max, n_points)
v0 = (u0 - 1)**2 + 1
v = 2 * (u0 - 1)
u = np.sign(v)
v = np.abs(v)

fig, ax = plt.subplots()

ax.plot(x, y, color='royalblue', alpha=0.5)
ax.quiver(u0, v0, u, v)

ax.set_xlim(x_min, x_max)
ax.set_ylim(x_min, x_max)

ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')

ax.set_aspect('equal')

plt.show()

import numpy as np

import matplotlib.pyplot as plt

x_min = -1

x_max = 3

x = np.linspace(x_min, x_max, 50)

y = (x - 1)**2 + 1

n_points = 21

u0 = np.linspace(x_min, x_max, n_points)

v0 = (u0 - 1)**2 + 1

v = 2 * (u0 - 1)

u = np.sign(v)

v = np.abs(v)

fig, ax = plt.subplots()

ax.plot(x, y, color='royalblue', alpha=0.5)

ax.quiver(u0, v0, u, v)

ax.set_xlim(x_min, x_max)

ax.set_ylim(x_min, x_max)

ax.spines['top'].set_visible(False)

ax.spines['right'].set_visible(False)

ax.spines['bottom'].set_position('zero')

ax.spines['left'].set_position('zero')

ax.set_aspect('equal')

plt.show()

2変数関数のgradient

勾配ベクトルの分布

2変数の関数について、gradientの分布を描いてみる。

関数として $z = x^2 + y^2$ を考えると、そのgradientは $( 2x, 2y)$ であり、ベクトル場は以下のようになる。

中心から外側に向かって増加しており、その傾きは外側ほど大きくなっている(ベクトルの長さが長くなり、コンター間隔は小さくなっている)。

import numpy as np
import matplotlib.pyplot as plt

lim_min = -1
lim_max = 1

x = np.linspace(lim_min, lim_max, 51)
y = np.linspace(lim_min, lim_max, 51)
x, y = np.meshgrid(x, y)
z = x * x + y* y

u0 = np.linspace(lim_min, lim_max, 21)
v0 = np.linspace(lim_min, lim_max, 21)
u0, v0 = np.meshgrid(u0, v0)
u = 2 * u0
v = 2 * v0

fig, ax = plt.subplots()

ax.contour(x, y, z)
ax.quiver(u0, v0, u, v)

ax.set_aspect('equal')

plt.show()

import numpy as np

import matplotlib.pyplot as plt

lim_min = -1

lim_max = 1

x = np.linspace(lim_min, lim_max, 51)

y = np.linspace(lim_min, lim_max, 51)

x, y = np.meshgrid(x, y)

z = x * x + y* y

u0 = np.linspace(lim_min, lim_max, 21)

v0 = np.linspace(lim_min, lim_max, 21)

u0, v0 = np.meshgrid(u0, v0)

u = 2 * u0

v = 2 * v0

fig, ax = plt.subplots()

ax.contour(x, y, z)

ax.quiver(u0, v0, u, v)

ax.set_aspect('equal')

plt.show()

勾配ベクトルの曲面上での分布

微分係数の時と同じように、曲面上にgradientoの分布を描いてみたのが以下の図。

2次元の曲線の時と同じく、勾配ベクトルが曲面に沿った傾きを表していることがわかる。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

lim_min = -1
lim_max = 1
n_points=15

x = np.linspace(lim_min, lim_max, n_points)
y = np.linspace(lim_min, lim_max, n_points)
x, y = np.meshgrid(x, y)
z = x * x + y* y

u0 = x
v0 = y
w0 = u0 * u0 + v0 * v0
u = 2 * u0
v = 2 * v0
w = np.sqrt(u * u + v * v)
(u, v) = (u, v) / w

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.plot_wireframe(x, y, z, color='cyan')
ax.quiver(u0, v0, w0, u, v, w, length=0.2, normalize=True)

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

lim_min = -1

lim_max = 1

n_points=15

x = np.linspace(lim_min, lim_max, n_points)

y = np.linspace(lim_min, lim_max, n_points)

x, y = np.meshgrid(x, y)

z = x * x + y* y

u0 = x

v0 = y

w0 = u0 * u0 + v0 * v0

u = 2 * u0

v = 2 * v0

w = np.sqrt(u * u + v * v)

(u, v) = (u, v) / w

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

ax.plot_wireframe(x, y, z, color='cyan')

ax.quiver(u0, v0, w0, u, v, w, length=0.2, normalize=True)

plt.show()

なお上記で(u, v) = (u, v) / wとしているが、 $\nabla f$ 方向で長さ1のベクトルを意味していて、これに対して $\nabla f$ をz方向の成分とすると曲面に沿ったベクトルとなる。

gradientの意味

2変数の場合

(4) $\begin{equation*} df(x, y) = \frac{\partial f}{\partial x} dx + \frac{\partial f}{\partial y} dy = f_x dx + f_y dy \end{equation*}$

gradientの方向の意味

先の微分係数の考え方から、勾配の各成分は関数の各点においてその成分方向に増加する。

したがって、gradientの方向は、関数が増加する方向を指している。

gradientの大きさの意味

コーシー・シュワルツの不等式から、

(5) $\begin{equation*} | df(x, y) | ^2 = \left( f_x dx + f_y dy \right) ^2 \le ( f_x^2 + f_y^2 ) ( dx^2 + dy^2) \end{equation*}$

ここで右辺の最初の項がgradientの絶対値だから、

(6) $\begin{align*} | df(x, y) | ^2 \le | \nabla f |^2 ( dx^2 + dy^2) \end{align*}$

${\boldsymbol d} = (dx, dy)$ 、 $| {\boldsymbol d} | = C$ とすると、 $f(x, y)$ の ${\boldsymbol d}$ 方向の変化量の大きさは $C | \nabla f |$ 以下である。

また等号が成立するためには、

(7) $\begin{gather*} \left( f_x dx + f_y dy \right) ^2 - ( f_x^2 + f_y^2 ) ( dx^2 + dy^2) = 0 \\ 2 f_x dx f_y dy - f_x^2 dy^2 - f_y^2 dx^2 = 0 \\ (f_x dy - f_y dx)^2 = 0 \\ (f_x , f_y) \cdot (dy, -dx) = 0 \end{gather*}$

これは、 $\nabla f$ と $(dx, dy)$ の2つのベクトルが平行であることを意味している。

すなわちgradientの方向は、関数の変化量の大きさが最も大きくなる方向を指している。

gradientの意味

上記のことを併せると、gradientは以下のように最も勾配が大きな方向を意味することがわかる。

gradientの方向は、関数がその点で最も大きく増加する方向
gradientと逆の方向は、関数がその点で最も大きく減少する方向

n変数の場合

gradientの方向と大きさの意味は、2変数の場合と同じ。

不等式で等号が成立する場合を考える。

(8) $\begin{equation*} (f_{x_1} dx_1 + \cdots + f_{x_n} dx_n)^2 - (f_{x_1} ^2 + \cdots + f_{x_n} ^2) (dx_1 ^2 + \cdots + dx_n ^2) = 0 \end{equation*}$

左辺第1項は、

(9) $\begin{equation*} \begin{array}{cccccccc} f_{x_1} ^2 dx_1 ^2 &+& f_{x_1} dx_1 f_{x_2} dx_2 &+& \cdots &+& f_{x_1} dx_1 f_{x_n} dx_n &+\\ f_{x_2} dx_2 f_{x_1} dx_1 &+& f_{x_2} ^2 dx_2 ^2 &+& \cdots &+& f_{x_2} dx_2 f_{x_n} dx_n &+\\ &&&& \cdots &&& \\ f_{x_n} dx_n f_{x_1} dx_1 &+& f_{x_n} dx_n f_{x_1} dx_1 &+& \cdots &+& f_{x_n} ^2 dx_n ^2 \end{array} \end{equation*}$

また左辺第2項は、

(10) $\begin{equation*} \begin{array}{cccccccc} f_{x_1} ^2 dx_1 ^2 &+& f_{x_1} ^2 dx_2 ^2 &+& \cdots &+& f_{x_1} ^2 dx_n ^2 &+\\ f_{x_2} ^2 dx_1 ^2 &+& f_{x_2} ^2 dx_2 ^2 &+& \cdots &+& f_{x_2} ^2 dx_n ^2 &+\\&&&& \cdots &&& \\ f_{x_n} ^2 dx_1 ^2 &+& f_{x_n} ^2 dx_2 ^2 &+& \cdots &+& f_{x_n} ^2 dx_n ^2 &+\\\end{array} \end{equation*}$

それぞれ各項が行列でいえば対角になっていることに留意しながら、(10)から(9)を差し引いて以下を得る。

(11) $\begin{equation*} \sum_{i=1}^{n} \sum_{j=i}^{n} \left( f_{x_i} dx_j - f_{x_j} dx_i \right) ^2 = 0 \end{equation*}$

これは以下のようにも表せる。

(12) $\begin{equation*} f_{x_i} dx_j = f_{x_j} dx_i \quad (i \ne j) \end{equation*}$

これは以下を意味する。

(13) $\begin{equation*} f_{x_i} : dx_i = f_{x_j} : dx_j \quad (i \ne j) \end{equation*}$

すなわちn変数の場合でも、 $\nabla f$ と ${\boldsymbol d}$ の方向が一致するときに $f$ の変化量が最も大きい。

コーシー・シュワルツの不等式

2020-02-26 / tau / コメントする

公式

Cauchy-Schwaltz inequality

(1) $\begin{equation*} \left( \sum_{i=1}^n a_i ^2 \right) \left( \sum_{i=1}^n b_i ^2 \right) \ge \left( \sum_{i=1}^n a_i b_i \right) ^2 \end{equation*}$

証明

n=2の場合

(2) $\begin{equation*} \left( a_1^2 + a_2^2 \right) \left(b_1^2 + b_2^2 \right) \ge \left( a_1 b_1 + a_2 b_2 \right) ^2 \end{equation*}$

(3) $\begin{align*} &\left( a_1^2 + a_2^2 \right) \left(b_1^2 + b_2^2 \right) - \left( a_1 b_1 + a_2 b_2 \right) ^2 \\ &= a_1 ^2 b_1^2 + a_1^2 b_2^2 + a_2^2 b_1^2 + a_2^2 b_2^2 - a_1^2 b_1^2 - 2 a_1 b_1 a_2 b_2 - a_2^2 b_2^2 \\ &= a_1^2 b_2^2 + a_2^2 b_1^2 - 2 a_1 b_1 a_2 b_2 \\ &= \left( a_1 b_2 - a_2 b_1 \right) ^2 \ge 0 \end{align*}$

nが任意の場合

2次方程式の判別式による方法

以下の2次方程式を考える。

(4) $\begin{equation*} \sum_{i=1}^n \left( a_i x + b_i \right)^2 = 0 \end{equation*}$

ここで関数 $f(x) = \sum_{i=1}^2 (a_i x + b_i)^2 \ge 0$ であり、上記の2次方程式の数は0個または1個である。

この方程式は以下のように変形できる。

(5) $\begin{equation*} \left( \sum a_i^2 \right) x^2 + 2 \left( \sum a_i b_i \right) x + \left( \sum b_i^2 \right) = 0 \end{equation*}$

もとの方程式の解の個数が0 or 1なので、上記の方程式の判別式から

(6) $\begin{gather*} 4 \left( \sum a_i b_i \right)^2 - 4 \left( \sum a_i^2 \right) \left( \sum b_i^2 \right) \le 0 \\ \therefore \left( \sum a_i^2 \right) \left( \sum b_i^2 \right) \ge \left( \sum a_i b_i \right)^2 \end{gather*}$

イメージ

${\boldsymbol a} = (a_1, \ldots , a_n)$ 、 ${\boldsymbol b} = (b_1, \ldots , b_n)$ とすると、ベクトルの内積となす角の関係から

(7) $\begin{gather*} \left( {\boldsymbol a}{\boldsymbol b} \right)^2 = \left( \sum a_i b_i \right)^2 = | {\boldsymbol a} |^2 | {\boldsymbol b} |^2 \cos^2 \theta \le | {\boldsymbol a} |^2 | {\boldsymbol b} |^2 = \left( \sum a_i^2 \right) \left( \sum a_i^2 \right) \end{gather*}$

ラグランジュの未定乗数法～等式制約条件

2020-02-23 / tau / コメントする

準備

長方形の面積最大化

たとえば2変数の問題として、長方形の周囲長 $L$ を一定として、その面積が最大となる長方形の形状と面積はどのようになるかを考える。この場合、長方形の辺の長さを $x, y$ とすると、問題は以下のように表せる。

(1) $\begin{equation*} \mathrm{maximize} \quad S(x, y) = xy \quad \mathrm{subject~to} \quad x+y = \frac{L}{2} \end{equation*}$

これは以下のように代数的に簡単に解けて、答えは正方形とわかる。

(2) $\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*}$

ただし変数の数が増えたり、目的関数や制約条件が複雑になると、解析的に解くのが面倒になる。

Lgrangeの未定乗数法による解

解法から先に示す。Lagrangeの未定乗数法では、目的関数 $L(x, y)$ に対して以下の問題となる。

(3) $\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*}$

$L(x, y, \lambda)$ を最大化するために、 $x, y, \lambda$ で偏微分した以下の方程式を設定する。

(4) $\begin{gather*} \frac{\partial L}{\partial x} = 0, \quad \frac{\partial L}{\partial y} = 0, \quad \frac{\partial L}{\partial \lambda} = 0 \end{gather*}$

これを計算すると

(5) $\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*}$

Lagrangeの未定乗数法の一般形

一般には、変数 $\boldsymbol{x} = (x_1, \ldots, x_n)$ について、目的関数 $f(\boldsymbol{x})$ を制約条件 $g(\boldsymbol{x})=0$ の下で最大化／最小化する問題として与えられる。

(6) $\begin{align*} & \mathrm{maximize/minimize} \quad f(\boldsymbol{x}) \\ & \mathrm{subject~to} \quad g(\boldsymbol{x}) = 0 \end{align*}$

この等式制約条件付き最大化／最小化問題は、以下のように $L(\boldsymbol{x}, \lambda)$ を導入して、連立方程式として表現される。

(7) $\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*}$

例題

例題1：平面と円

平面 $x + y - 1$ について、制約条件 $f(x, y) = x^2 + y^2$ の下での極値を求める。

(8) $\begin{align*} & \mathrm{minimize} \quad f(x, y) = x + y - 1 \\ & \mathrm{subject~to} \quad x^2 + y^2 = 1 \end{align*}$

lagrangeの未定乗数を導入して問題を定式化すると以下のようになる。

(9) $\begin{equation*} L(x, y, \lambda) = x + y - 1 - \lambda (x^2 + y^2 - 1) \end{equation*}$

(10) $\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*}$

この連立方程式を解くと以下のようになり、解として2つの極値を得るが、それらは最大値と最小値に相当する。

(11) $\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*}$

これを目的関数のコンターと制約条件の線で表すと以下の通り。

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()

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：凸関数と直線

下に凸な関数 $f(x, y) = x^2 + y^2$ について、直線 $x + y = 1$ の制約条件下での最小値を求める。

(12) $\begin{align*} & \mathrm{minimize} \quad f(x, y) = x^2 + y^2 \\ & \mathrm {subject~to} \quad x + y - 1 = 0 \end{align*}$

ここでlagrangeの未定乗数を導入して問題を定式化すると以下のようになる。

(13) $\begin{equation*} L(x, y, \lambda) = x^2 + y^2 - \lambda (x + y - 1) \end{equation*}$

(14) $\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*}$

この連立方程式を解くと以下のようになり、解は最小値1つとなる。

(15) $\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*}$

これを目的関数のコンターと制約条件の線で表すと以下の通り。

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()

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次元で表示すると以下のようになる。青い曲面が目的関数で、赤い直線が制約条件となる。最適化問題は、制約条件を満たす曲面上の点(図中、赤い放物線)の最小値を求めることになる。

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()

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) $\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*}$

さらにgradientで表すと

(17) $\begin{gather*} \nabla f = \lambda \nabla g \\ g(x_1, \ldots, x_n) = 0 \end{gather*}$

すなわちこの式の解 $(x_1, \ldots, x_n)$ は、制約条件である $g(x_1, \ldots, x_n)=0$ を満足し、その曲線上にある。さらに解の点において $f(x_1, \ldots, x_n)$ の勾配ベクトルとゼロ平面上における $g(x_1, \ldots, x_n)$ の勾配ベクトルが平行になる。これはゼロ平面上の解の点において制約条件の曲線と $f(x_1, \ldots, x_n)$ のコンターのうち特定の曲線が接するのと同義であり、この点は停留点(stationary point)である。

つまりこのような停留点を発見する手順は、制約条件を満たす(制約条件の線上にある)点のうち、その点において目的関数のgradientと制約条件の関数のgradientが平行となる点を求めるということになる。再度これを式で表すと、

(18) $\begin{gather*} \nabla f(\boldsymbol{x}) = \lambda \nabla g(\boldsymbol{x}) \\ g(\boldsymbol{x}) = 0 \end{gather*}$

となるが、これをLangrange関数 $L=f - \lambda g$ と定義したうえで各変数で偏微分したものをゼロと置いた方程式を解くと表現している。未定乗数λは停留点における目的関数のgradientと制約条件の関数のgradientの比を表している。

λの符号

λの符号に意味があるかどうか。

たとえば、以下の制約条件付き最適化問題を考える。

(19) $\begin{align*} & \mathrm{minimize} \quad f(x, y) = x^2 + y^2 \\ & \mathrm {subject~to} \quad x + y - 2 = 0 \end{align*}$

(20) $\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*}$

この問題の制約条件を変更して以下のようにした場合。

(21) $\begin{align*} & \mathrm{minimize} \quad f(x, y) = x^2 + y^2 \\ & \mathrm {subject~to} \quad - x - y + 2 = 0 \end{align*}$

(22) $\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*}$

このように、制約条件の正負を反転するとλの符号が逆になるが解は変わらない。これを表したのが以下の図。

等式制約条件の場合、制約条件の線上でgradientが平行になりさえすればよいので、λ符号(制約条件式の正負)には拘らなくてよい。

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()

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) $\begin{align*} & \mathrm{maxmize} \quad f(x, y) = x^3 + y^3 \\ & \mathrm{subject~to} \quad g(x, y) = x - y \end{align*}$

(24) $\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*}$

第1式から第2式を引き、第3式を適用して、

(25) $\begin{gather*} \3x^2 - 3y^2 - 2\lambda = 0 \\ 3(x + y)(x - y) = 2\lambda = 0 \\ \therefore x = y = \lambda = 0 \end{gather*}$

ここで $z=f(x, y)$ の曲面と制約条件 $g(x, y)=0$ に対する曲面上の軌跡を描くと下図左のようになる。下図の右は $t = x = y$ として曲線を表したもので、 $x = y = 0$ で勾配は水平になっており、停留点ではあるが極大／極小となっていない。

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()

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()

参考サイト

本記事をまとめるにあたって、下記サイトが大変参考になったことに感謝したい。

その２極値を求めるラグランジュの未定乗数法～○×（まるぺけ）つくろーどっとコム

定義

L1正則化の意味

準備

場合分け

劣微分の導入

参考サイト

概要

例題

停留点が制約条件の範囲内の場合

停留点が制約条件の範囲外の場合

λの符号について

有界でない例

最小化問題でgradientが逆向き(極値が制約条件外)

最小化問題でgradientが同じ向き(極値が制約条件内)

最大化問題でgradientが同じ向き(極値が制約条件範囲外)

最大化問題でgradientが逆向き(極値が制約条件範囲内)

λの符号について

KKT条件

参考サイト

定義式

計算例

次数2の場合

次数3の場合

記号の定義

ベクトル・行列をスカラーで微分

スカラーをベクトルで微分

スカラーを行列で微分

ベクトルをベクトルで微分

公式

一般形

単位行列

合成関数

積の微分

行列の積のスカラーによる微分

一次～二次形式の微分

Axの形式

x^2の形式

xTAxの形式

定義

性質

単独の行列

転置の転置

逆行列

行列式

行列演算

線形性

積

行列とベクトル

行列とベクトルの積

内積

定義

準備～1変数関数の微分係数

微分係数の符号だけを見た場合

微分係数の符号を方向としてみた場合

微分係数の曲線上での分布

2変数関数のgradient

勾配ベクトルの分布

勾配ベクトルの曲面上での分布

gradientの意味

2変数の場合

gradientの方向の意味

gradientの大きさの意味

gradientの意味

n変数の場合

公式

証明

n=2の場合

nが任意の場合

2次方程式の判別式による方法

イメージ

準備

長方形の面積最大化

Lgrangeの未定乗数法による解

Lagrangeの未定乗数法の一般形

例題

例題1：平面と円

例題2：凸関数と直線

幾何学的説明

λの符号

停留点が極致とならない例

x^TAxの形式