Python3 – 配列要素の重複を除く

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

リスト要素の重複を除く

リストをset()関数の引数にすると、重複する要素がなくなり全ての要素がユニークになる。

import random

a = [random.randint(0, 2) for n in range(20)]

print(a)
print(set(a))

# [2, 2, 2, 2, 1, 0, 1, 0, 0, 0, 0, 1, 1, 2, 1, 1, 0, 2, 2, 0]
# {0, 1, 2}

import random

a = [random.randint(0, 2) for n in range(20)]

print(a)

print(set(a))

# [2, 2, 2, 2, 1, 0, 1, 0, 0, 0, 0, 1, 1, 2, 1, 1, 0, 2, 2, 0]

# {0, 1, 2}

ただし結果は集合なので、これをリスト化するにはlist()関数を使う。

set()関数の結果、要素は昇順に並んでいるが、hashの計算方法によって必ずしも昇順になると決まってはいない。そこでsorted()関数でリストをソートしておく(生成されたリストをほかで再利用しないなら、sort()メソッドを使ってもよい)。

import random

a = [random.randint(0, 2) for n in range(20)]

print(a)
print(sorted(list(set(a))))

# [2, 2, 2, 2, 1, 0, 1, 0, 0, 0, 0, 1, 1, 2, 1, 1, 0, 2, 2, 0]
# [0, 1, 2]

import random

a = [random.randint(0, 2) for n in range(20)]

print(a)

print(sorted(list(set(a))))

# [2, 2, 2, 2, 1, 0, 1, 0, 0, 0, 0, 1, 1, 2, 1, 1, 0, 2, 2, 0]

# [0, 1, 2]

ndarrayの要素の重複を除く

元の配列がndarrayで与えられた場合でも、set()関数を適用すると結果は集合となるが、ndarrayの生成時に引数を集合とすると配列として生成されずに集合のまま。

import numpy as np
import numpy.random as rnd

a = rnd.randint(0, 3, 20)

print(a)
print(set(a))
print(np.array(set(a)))

# [2 2 1 2 2 0 1 1 0 1 0 1 2 1 0 1 0 2 2 0]
# {0, 1, 2}
# {0, 1, 2}

import numpy as np

import numpy.random as rnd

a = rnd.randint(0, 3, 20)

print(a)

print(set(a))

print(np.array(set(a)))

# [2 2 1 2 2 0 1 1 0 1 0 1 2 1 0 1 0 2 2 0]

# {0, 1, 2}

そこで、list()関数でいったん集合をリスト化してからndarrayにする必要がある。

また、要素の昇順を保証するためにnp.array()関数でソートしておく。生成した配列を再利用しないのであれば、ndarrayのsort()メソッドを使ってもよい。

import numpy as np
import numpy.random as rnd

a = rnd.randint(0, 3, 20)

print(a)
print(np.array(list(set(a))))

# [2 1 2 0 0 0 2 2 2 1 0 1 0 2 0 0 2 0 2 2]
# [0 1 2]

import numpy as np

import numpy.random as rnd

a = rnd.randint(0, 3, 20)

print(a)

print(np.array(list(set(a))))

# [2 1 2 0 0 0 2 2 2 1 0 1 0 2 0 0 2 0 2 2]

# [0 1 2]

利用例～クラス値を持つデータの分類

たとえば多数のデータの特性値とクラス区分が配列で与えられた場合、クラスごとにマーカーの形や色を変えてプロットするなど、クラスごとに元のデータを分けて処理したい場合。

以下のようにset()関数で重複を除いてループ処理できる。

import numpy as np

X = np.array(['swallow', 'carp', 'dog', 'horse', 'hawk', 'bonito'])
y = np.array([1, 0, 3, 3, 1, 0])

for cls in set(y):
    print(X[y==cls])

# ['carp' 'bonito']
# ['swallow' 'hawk']
# ['dog' 'horse']

import numpy as np

X = np.array(['swallow', 'carp', 'dog', 'horse', 'hawk', 'bonito'])

y = np.array([1, 0, 3, 3, 1, 0])

for cls in set(y):

print(X[y==cls])

# ['carp' 'bonito']

# ['swallow' 'hawk']

# ['dog' 'horse']

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

参考サイト

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

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

matplotlib.pyplot.quiver – ベクトル場

2020-02-23 / tau / 3件のコメント

概要

matplotlib.pyplot.quiver()はベクトル場を可視化する。基本的なパラメーターは以下の通り。

quiver(X, Y, U, V, [C]): X, Yはベクトルの開始点、U, Vはベクトルの成分、Cはベクトルの大きさに応じたカラーマップ上の色をつけるための配列。

単一のベクトルの描画例

以下の例では、始点の位置と成分を1つずつ指定してベクトルを描画している。デフォルトではベクトルのスケールは描画領域に対して自動的に調節されるが、ここでは描画領域のスケールと同じになるようパラメーターを設定している。

matplotlibのドキュメントでは、scale_unitsのところに以下のように書かれている。

“To plot vectors in the x-y plane, with u and v having the same units as x and y, use angles='xy', scale_units='xy', scale=1”

import numpy as np
import matplotlib.pyplot as plt

fig, ax = plt.subplots()

ax.quiver(1, 2, color='blue', angles='xy', scale_units='xy', scale=1)
ax.quiver(1, 2, 2, -1, color='red', angles='xy', scale_units='xy', scale=1)

ax.set_xlim(0, 3)
ax.set_ylim(0, 3)

ax.set_xticks(np.arange(0, 4, 1))
ax.set_yticks(np.arange(0, 4, 1))

ax.grid()
ax.set_aspect('equal')

plt.show()

import numpy as np

import matplotlib.pyplot as plt

fig, ax = plt.subplots()

ax.quiver(1, 2, color='blue', angles='xy', scale_units='xy', scale=1)

ax.quiver(1, 2, 2, -1, color='red', angles='xy', scale_units='xy', scale=1)

ax.set_xlim(0, 3)

ax.set_ylim(0, 3)

ax.set_xticks(np.arange(0, 4, 1))

ax.set_yticks(np.arange(0, 4, 1))

ax.grid()

ax.set_aspect('equal')

plt.show()

ベクトル場の描画例

以下の例では、xy平面上の位置に応じた成分を持つベクトルを描画している。関数のgradientのイメージ。

描画にあたって、開始点の座標とベクトルの成分をmeshgridで生成している。

1つ目の図は単に開始点と成分を与えただけで、単一の色で、ベクトルのスケールは自動調節されている。

2つ目の図はスケールと色付けのための配列を指定し、ベクトルの大きさに応じてcolormapで色を付けている。

import numpy as np
import matplotlib.pyplot as plt

x = y = np.arange(-5, 6)
u = 2 * x
v = 3 * y

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

fig, axes = plt.subplots(1, 2, figsize=(12, 4.8))

lim = 8
for ax in axes:
    ax.set_xlim(-lim, lim)
    ax.set_ylim(-lim, lim)

    ax.set_xticks(np.arange(-lim, lim, 1))
    ax.set_yticks(np.arange(-lim, lim, 1))

    ax.grid()
    ax.set_aspect('equal')

C = np.sqrt(u * u + v * v)
axes[0].quiver(x, y, u, v)
axes[1].quiver(x, y, u, v, C, scale=100, cmap='Blues')

plt.show()

import numpy as np

import matplotlib.pyplot as plt

x = y = np.arange(-5, 6)

u = 2 * x

v = 3 * y

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

u, v = np.meshgrid(u, v)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.8))

lim = 8

for ax in axes:

ax.set_xlim(-lim, lim)

ax.set_ylim(-lim, lim)

ax.set_xticks(np.arange(-lim, lim, 1))

ax.set_yticks(np.arange(-lim, lim, 1))

ax.grid()

ax.set_aspect('equal')

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

axes[0].quiver(x, y, u, v)

axes[1].quiver(x, y, u, v, C, scale=100, cmap='Blues')

plt.show()

不等式の公式

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

相加相乗平均の不等式

一般系

(1) $\begin{equation*} a_i \ge 0 \quad \Rightarrow \quad \sum_{i=1}^{n} a_i \ge n \sqrt[n]{\prod_{i=1}^{n} a_i} \end{equation*}$

2項の場合

(2) $\begin{equation*} a, b \ge 0 \quad \Rightarrow \quad a + b \ge 2 \sqrt{ab} \end{equation*}$

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

Cauchy-Shwaltzの不等式一般系

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

2項の場合

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

Axes.spines～軸の設定

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

概要

グラフのx軸、y軸の位置や表示の有無については、Axesオブジェクトのspinesプロパティーで制御する。

spinesは辞書型でbottom、top、left、rightのキーで対象を指定し、表示位置はset_position()メソッド、表示の有無はset_visible()で操作する。

軸の指定

spines['bottom']とspines['left']は下と左の軸で、軸の値が表示される。

spines['top']とspines['right']は上と右の軸で、ただ線が引かれるだけ。

各軸に対して、set_positon()、set_visible()の各メソッドを実行して、位置や可視／不可視を設定する。

軸の表示・非表示

set_visible(False)で軸を非表示にする。

以下の例では、上の軸と右の軸を非表示にしている。

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 1, 50)

fig , axes = plt.subplots(1, 2)

for ax in axes:
    ax.plot(x, 1 - x)

    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1)

    ax.grid()
    ax.set_aspect('equal')

axes[1].spines['right'].set_visible(False)
axes[1].spines['top'].set_visible(False)

plt.show()

import numpy as np

import matplotlib.pyplot as plt

x = np.linspace(0, 1, 50)

fig , axes = plt.subplots(1, 2)

for ax in axes:

ax.plot(x, 1 - x)

ax.set_xlim(0, 1)

ax.set_ylim(0, 1)

ax.grid()

ax.set_aspect('equal')

axes[1].spines['right'].set_visible(False)

axes[1].spines['top'].set_visible(False)

plt.show()

ゼロ位置／中央に軸位置を設定

set_position(‘zero’)でゼロの位置に、set_position(‘center’)で描画位置の中央に軸をセットできる。set_visible()と組み合わせて使うケースが多そう。

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-1, 2, 50)

fig , axes = plt.subplots(1, 2)

for ax in axes:
    ax.plot(x, 1 - x)

    ax.set_xlim(-1, 2)
    ax.set_ylim(-1, 2)

    ax.set_xticks(np.arange(-1, 2.5, 0.5))
    ax.set_yticks(np.arange(-1, 2.5, 0.5))

    ax.grid()
    ax.set_aspect('equal')

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

axes[0].spines['bottom'].set_position('zero')
axes[0].spines['left'].set_position('zero')

axes[1].spines['bottom'].set_position('center')
axes[1].spines['left'].set_position('center')

plt.show()

import numpy as np

import matplotlib.pyplot as plt

x = np.linspace(-1, 2, 50)

fig , axes = plt.subplots(1, 2)

for ax in axes:

ax.plot(x, 1 - x)

ax.set_xlim(-1, 2)

ax.set_ylim(-1, 2)

ax.set_xticks(np.arange(-1, 2.5, 0.5))

ax.set_yticks(np.arange(-1, 2.5, 0.5))

ax.grid()

ax.set_aspect('equal')

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

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

axes[0].spines['bottom'].set_position('zero')

axes[0].spines['left'].set_position('zero')

axes[1].spines['bottom'].set_position('center')

axes[1].spines['left'].set_position('center')

plt.show()

軸の位置の数値指定

set_position()の引数として、タプルで('指定方法', 値)の形で与える。

指定方法	値
data	各軸を配置するx、yの値。
outward	単位はポイントで、正なら描画領域の内側、負なら外側に配置。
axes	描画領域の高さ・幅に対する割合。

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 10, 50)

fig , axes = plt.subplots(2, 2)

for row in axes:
    for ax in row:
        ax.plot(x, 10 - x)

        ax.set_xlim(0, 10)
        ax.set_ylim(0, 10)

        ax.set_xticks(np.arange(0, 12, 2))
        ax.set_yticks(np.arange(0, 12, 2))

        ax.grid()
        ax.set_aspect('equal')

ax = axes[0, 1]
ax.set_title('data')
ax.spines['bottom'].set_position(('data', 1))
ax.spines['left'].set_position(('data', 2))
ax.spines['top'].set_position(('data', 7))
ax.spines['right'].set_position(('data', 8))

ax = axes[1, 0]
ax.set_title('outward')
ax.spines['bottom'].set_position(('outward', 5))
ax.spines['left'].set_position(('outward', -10))
ax.spines['top'].set_position(('outward', 15))
ax.spines['right'].set_position(('outward', -20))

ax = axes[1, 1]
ax.set_title('axes')
ax.spines['bottom'].set_position(('axes', 0.1))
ax.spines['left'].set_position(('axes', 0.2))
ax.spines['top'].set_position(('axes', 0.7))
ax.spines['right'].set_position(('axes', 0.8))

fig.subplots_adjust(hspace=0.5)
plt.show()

import numpy as np

import matplotlib.pyplot as plt

x = np.linspace(0, 10, 50)

fig , axes = plt.subplots(2, 2)

for row in axes:

for ax in row:

ax.plot(x, 10 - x)

ax.set_xlim(0, 10)

ax.set_ylim(0, 10)

ax.set_xticks(np.arange(0, 12, 2))

ax.set_yticks(np.arange(0, 12, 2))

ax.grid()

ax.set_aspect('equal')

ax = axes[0, 1]

ax.set_title('data')

ax.spines['bottom'].set_position(('data', 1))

ax.spines['left'].set_position(('data', 2))

ax.spines['top'].set_position(('data', 7))

ax.spines['right'].set_position(('data', 8))

ax = axes[1, 0]

ax.set_title('outward')

ax.spines['bottom'].set_position(('outward', 5))

ax.spines['left'].set_position(('outward', -10))

ax.spines['top'].set_position(('outward', 15))

ax.spines['right'].set_position(('outward', -20))

ax = axes[1, 1]

ax.set_title('axes')

ax.spines['bottom'].set_position(('axes', 0.1))

ax.spines['left'].set_position(('axes', 0.2))

ax.spines['top'].set_position(('axes', 0.7))

ax.spines['right'].set_position(('axes', 0.8))

fig.subplots_adjust(hspace=0.5)

plt.show()

ヘッセの標準形～点と平面の距離

2020-02-16 / tau / 2件のコメント

2次元の場合

直線 $l$ 、点 $Q$ の距離を考える。直線と各点の記号、座標を以下のように定義する。

(1) $\begin{equation*} l:ax + by + c = 0 \end{equation*}$

Step-1：直線に直交するベクトル

まず、ベクトル $(a, b)$ が直線 $l$ に直交することを示す。直線は以下のように媒介変数表示できて、ベクトル $(u_x , u_y)$ は直線に平行なベクトル。

(2) $\begin{equation*} \left\{ \begin{array}{lll} x &=& u_x t + x_0 \\ y &=& u_y t+ y_0 \end{array} \right. \end{equation*}$

これを直線の式に代入して、

(3) $\begin{gather*} a (u_x t + x_0) + b (u_y t + y_0) + c = 0 \\ (a u_x + b u_y) t + (a x_0 + b y_0 + c) = 0 \end{gather*}$

ここで任意の $t$ に対して上式が成り立つことから、 $a u_x + b u_y = 0$ となり、ベクトル $(a, b)$ は直線に垂直であることが示された。

別解：点と直線の距離をパラメーター（媒介変数）によって愚直に求める方法

Step-2：法線ベクトルとの平行条件による導出

与えられた点 $P$ から直線 $l$ への垂線の足を $H$ とすると、 $\overrightarrow{PH} \parallel (a, b)$ なので、以下が成り立つ。

(4) $\begin{equation*} \left| \overrightarrow{PH} \cdot (a, b) \right| = \left| \overrightarrow{PH} \right| \cdot | (a, b) | \end{equation*}$

ここで $(x_h, y_h)$ は直線 $l$ 上にあることを考慮し、上式の左辺を以下のように変形できる。

(5) $\begin{equation*} \begin{array}{lll} \left| \overrightarrow{PH} \cdot (a, b) \right| &=& | (x_p - x_h ) a + (y_p - y_h) b | \\ &=& | a x_p + b y_p - (a x_h + b y_h) | \\ &=& | a x_p + b y_p + c | \end{equation*}$

これより

(6) $\begin{equation*} \left| \overrightarrow{PH} \right| = \frac{\left| a x_p + b y_p + c \right|}{\sqrt{a^2 + b^2}} \end{equation*}$

3次元の場合

3次元平面の式

3次元空間内の平面は、たとえば以下のように表すことができる。

(7) $\begin{equation*} \pi : w_x x + w_y y + w_z z + w_0 = 0 \end{equation*}$

一方、3次元平面上の点と法線ベクトル ${\boldsymbol n} = (x_x, n_y, n_z)$ との直行条件から、以下のようにも表現できる。

(8) $\begin{gather*} {\boldsymbol n} ({\boldsymbol x} - {\boldsymbol x_0} ) = 0 \\ (n_x, n_y, n_z) \cdot (x - x_0, y - y_0, z - z_0) = 0 \\ n_x x + n_y y + n_z z + (-n_x x_0 -n_y y_0 -n_z z_0) = 0 \end{gather*}$

上記2つの式より、ベクトル ${\boldsymbol w} = (w_x, w_y, w_z)$ は平面に対する法線ベクトルであることがわかる。

法線ベクトルとの平行条件

この法線ベクトルがベクトル $\left| \overrightarrow{PH} \right|$ と平行であることから、

(9) $\begin{equation*} \left| \overrightarrow{PH} \cdot {\boldsymbol w} \right| = \left| \overrightarrow{PH} \right| \cdot \left| {\boldsymbol w} \right| \end{equation*}$

上式の左辺は以下のように変形できる。

(10) $\begin{equation*} \begin{array}{lll} \left| \overrightarrow{PH} \cdot {\boldsymbol w} \right| &=& \left| (x_p - x_h, y_p - y_h, z_p - z_h) \cdot (w_x, w_y, w_z) \right| \\ &=& \left| w_x x_p + w_y y_p + w_z z_p - (w_x x_h +w_y y_h + w_z z_h) \right| \\ &=& \left| w_x x_p + w_y y_p + w_z z_p + w_0 \right| \end{array} \end{equation*}$

以上のことから、点 ${\boldsymbol x}_p$ から三次元平面 $\pi$ への距離については、以下で表される。

(11) $\begin{equation*} \left| \overrightarrow{PH} \right| = \frac{\left| w_x x_p + w_y y_p + w_z z_p + w_0 \right|}{| {\boldsymbol w} |} \end{equation*}$

多次元の場合

n次元の超平面を以下の式で与える。

(12) $} \begin{equation*} {\boldsymbol w} \cdot {\boldsymbol x} + w_0 = 0 \; \Leftrightarrow \; w_0 + w_1 x_1 + \cdots + w_n x_n = 0 \end{equation*}$

このとき、これまでと同様の考え方により、点 ${\boldsymbol x}_p (x_{p1}, \ldots , x_{pn})$ と上記の超平面との距離は以下で表される。

(13) $\begin{equation*} \left| \overrightarrow{PH} \right| = \frac{ {\boldsymbol w} \cdot {\boldsymbol x} + w_0 }{ \| {\boldsymbol w} \|} \end{equation*}$

Python3 – コレクションのソート

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

リストのソート

`sort()`メソッドは破壊的処理

sort()はリストのメソッドで、元のリストの内容を変更する(破壊的処理)。メソッドの実行結果はNone。

降順にソートしたいときは、引数reverseをTrueで指定。

lst = [3, 2, 1, 5, 4]

print(lst.sort())
print(lst)

lst.sort(reverse=True)
print(lst)

# None
# [1, 2, 3, 4, 5]
# [5, 4, 3, 2, 1]

lst = [3, 2, 1, 5, 4]

print(lst.sort())

print(lst)

lst.sort(reverse=True)

print(lst)

# None

# [1, 2, 3, 4, 5]

# [5, 4, 3, 2, 1]

`sorted()`関数は非破壊的処理

sorted()関数は引数のリストのソート結果を返す。元のリストの内容は変更されない(非破壊的処理)。

降順ソートの指定はsort()メソッドと同じ。

lst = [3, 2, 1, 5, 4]

print(sorted(lst))
print(lst)

print(sorted(lst, reverse=True))

print()

# [1, 2, 3, 4, 5]
# [3, 2, 1, 5, 4]
# [5, 4, 3, 2, 1]

lst = [3, 2, 1, 5, 4]

print(sorted(lst))

print(lst)

print(sorted(lst, reverse=True))

print()

# [1, 2, 3, 4, 5]

# [3, 2, 1, 5, 4]

# [5, 4, 3, 2, 1]

文字列は辞書順でソートされる

lst = ['ca', 'ba', 'aa', 'bb', 'ab']

print(sorted(lst))

print()

# ['aa', 'ab', 'ba', 'bb', 'ca']

lst = ['ca', 'ba', 'aa', 'bb', 'ab']

print(sorted(lst))

print()

# ['aa', 'ab', 'ba', 'bb', 'ca']

`ndarray`の場合の注意

`sorted()`はそのままでは`ndarray`にならない

ndarrayをsorted()関数の引数にすると、エラーにはならないが結果はリストで返されるため、配列への変換が必要。

import numpy as np

a = np.array([3, 2, 1, 5, 4])
print(sorted(a))
print(np.array(sorted(a)))

print()

# [1, 2, 3, 4, 5]
# [1 2 3 4 5]

import numpy as np

a = np.array([3, 2, 1, 5, 4])

print(sorted(a))

print(np.array(sorted(a)))

print()

# [1, 2, 3, 4, 5]

# [1 2 3 4 5]

`numpy.sort()`は非破壊的に`ndarray`をソートできる

numpy.sort()関数は、引数のndarrayのソート結果を返し、元のndarrayは変更しない。リストの場合のsorted()関数と同じ動作。

a = np.array([3, 2, 1, 5, 4])

print(np.sort(a))
print(a)

# [1 2 3 4 5]
# [3 2 1 5 4]

a = np.array([3, 2, 1, 5, 4])

print(np.sort(a))

print(a)

# [1 2 3 4 5]

# [3 2 1 5 4]

`ndarray`の`sort()`メソッドは破壊的

ndarrayのsort()メソッドは、元の配列の内容を書き換える。リストのsort()メソッドと同じ挙動で、実行結果の戻り値はNone。

a = np.array([3, 2, 1, 5, 4])
print(a.sort())
print(a)

# None
# [1 2 3 4 5]

a = np.array([3, 2, 1, 5, 4])

print(a.sort())

print(a)

# None

# [1 2 3 4 5]

辞書のソート

今後

Python3 – zipによる複数リストの並行ループ

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

2つのリストの要素を並行して取得しつつ処理したい場合、zip()関数を用いる。

names = ['Jane', 'Bill', 'Lucy', 'Amanda']
ages = [34, 18, 25, 44]

for name, age in zip(names, ages):
    print("{} is {} years old.".format(name, age))

# Jane is 34 years old.
# Bill is 18 years old.
# Lucy is 25 years old.
# Amanda is 44 years old.

names = ['Jane', 'Bill', 'Lucy', 'Amanda']

ages = [34, 18, 25, 44]

for name, age in zip(names, ages):

print("{} is {} years old.".format(name, age))

# Jane is 34 years old.

# Bill is 18 years old.

# Lucy is 25 years old.

# Amanda is 44 years old.

zip()関数は、引数で与えた複数のコレクションの要素が対になったタプルのイテレーターを返す。各コレクションの長さが異なる場合、イテレーターの長さは最も短いコレクションの長さとなり、それ以降の各コレクションの要素は無視される。

lst1 = ['A', 'B', 'C']
lst3 = ['zero', 'one', 'two', 'three']
lst2 = [0, 1, 2, 3, 4]

z = zip(lst1, lst2, lst3)

for t in z:
    print(t)

# ('A', 0, 'zero')
# ('B', 1, 'one')
# ('C', 2, 'two')

lst1 = ['A', 'B', 'C']

lst3 = ['zero', 'one', 'two', 'three']

lst2 = [0, 1, 2, 3, 4]

z = zip(lst1, lst2, lst3)

for t in z:

print(t)

# ('A', 0, 'zero')

# ('B', 1, 'one')

# ('C', 2, 'two')

リスト要素の重複を除く

ndarrayの要素の重複を除く

利用例～クラス値を持つデータの分類

定義

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

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

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

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

2変数関数のgradient

勾配ベクトルの分布

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

gradientの意味

2変数の場合

gradientの方向の意味

gradientの大きさの意味

gradientの意味

n変数の場合

公式

証明

n=2の場合

nが任意の場合

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

イメージ

準備

長方形の面積最大化

Lgrangeの未定乗数法による解

Lagrangeの未定乗数法の一般形

例題

例題1：平面と円

例題2：凸関数と直線

幾何学的説明

λの符号

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

参考サイト

概要

単一のベクトルの描画例

ベクトル場の描画例

相加相乗平均の不等式

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

概要

軸の指定

軸の表示・非表示

ゼロ位置／中央に軸位置を設定

軸の位置の数値指定

2次元の場合

Step-1：直線に直交するベクトル

Step-2：法線ベクトルとの平行条件による導出

3次元の場合

3次元平面の式

法線ベクトルとの平行条件

多次元の場合

リストのソート

sort()メソッドは破壊的処理

sorted()関数は非破壊的処理

文字列は辞書順でソートされる

ndarrayの場合の注意

sorted()はそのままではndarrayにならない

numpy.sort()は非破壊的にndarrayをソートできる

ndarrayのsort()メソッドは破壊的

辞書のソート

`sort()`メソッドは破壊的処理

`sorted()`関数は非破壊的処理

`ndarray`の場合の注意

`sorted()`はそのままでは`ndarray`にならない

`numpy.sort()`は非破壊的に`ndarray`をソートできる

`ndarray`の`sort()`メソッドは破壊的