numpyで数理最適化 メモ

• サポートページ : https://www.sbcr.jp/product/4797393965/
• 第1版 正誤表 : https://s3-ap-northeast-1.amazonaws.com/sbcr-dl-tec/2018-12-12-93965-%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92%E3%81%AE%E3%82%A8%E3%83%83%E3%82%BB%E3%83%B3%E3%82%B9/errata.pdf

from IPython.display import HTML

HTML('''
<script>  
code_show=true; 
function code_toggle() {
  if (code_show){
    $(\'div.input\').hide();
  } else {
    $(\'div.input\').show();
  }
  code_show = !code_show
}  
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()">
    <input type="submit" value="Toggle code">
</form>
''')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
plt.rcParams['font.size']=15

def plt_legend_out():
    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0)

ブロードキャスト

import time

a = np.arange(10000000)

time_s1 = time.time()
b = np.array([x**2 for x in a])
time_e1 = time.time() - time_s1

b = np.empty(a.shape[0])
time_s2 = time.time()
for i,x in enumerate(a):
    b[i] = x**2
time_e2 = time.time() - time_s2

df_plot = pd.DataFrame({'value':[time_e1,time_e2],'name':['ブロードキャスト','非ブロードキャスト']})
sns.barplot(data=df_plot,x='name',y='value')
plt.ylabel('time [s]')
plt.xlabel('')
plt.show()

線形代数

$Ax=b$

$
\left(
\begin{array}{ccc}
3 & 1 & 1 \
1 & 2 & 1 \
0 & -1 & 1
\end{array}
\right)
\left(
\begin{array}{ccc}
x \
y \
z
\end{array}
\right)=
\left(
\begin{array}{ccc}
1 \
2 \
3
\end{array}
\right)
$

a = np.array([[3,1,1],[1,2,1],[0,-1,1]])
a

array([[ 3,  1,  1],
       [ 1,  2,  1],
       [ 0, -1,  1]])

b = np.array([1,2,3])
b

array([1, 2, 3])

np.linalg.solve(a,b)

array([-0.57142857, -0.14285714,  2.85714286])

LU分解。逆行列を計算する方法よりも高速らしい。

$A=PLU$

$PLUx=b$

乱数

シード有無の比較

def throw_dice1(n):
    return np.random.randint(1, 7, size=n).sum()

def throw_dice2(n, random_seed=10):
    np.random.seed(random_seed)
    return np.random.randint(1, 7, size=n).sum()

print('シード指定なし：'+str([throw_dice1(10) for x in range(10)]))
print('シード指定あり：'+str([throw_dice2(10) for x in range(10)]))

シード指定なし：[30, 42, 37, 28, 27, 26, 46, 35, 26, 38]
シード指定あり：[34, 34, 34, 34, 34, 34, 34, 34, 34, 34]

関数から呼び出す

良くない例

# よくない実装の例
class Dice:
    def __init__(self):
        np.random.seed(0)
        self.sum_ = 0

    def throw(self):
        self.sum_ += np.random.randint(1, 7)

    def get_sum(self):
        return self.sum_

良い例

# コンストラクタの引数に種を与えると結果に再現性がある例
class Dice:
    def __init__(self, random_seed=None):
        self.random_state_ = np.random.RandomState(random_seed)
        self.sum_ = 0

    def throw(self):
        self.sum_ += self.random_state_.randint(1, 7)

    def get_sum(self):
        return self.sum_

import dice2a
import dice2b

print('良くない例')

l = []
for __ in range(10):
    d1 = dice2a.Dice()
    for _ in range(10):
        d1.throw()
    l.append(d1.get_sum())

print('v1 : '+str(l))

l = []
for __ in range(10):
    d2 = dice2a.Dice()
    for _ in range(10):
        d2.throw()
    l.append(d1.get_sum())

print('v2 : '+str(l))

l1 = []
l2 = []
for __ in range(10):
    d1 = dice2a.Dice()
    d2 = dice2a.Dice()
    for _ in range(10):
        d1.throw()
        d2.throw()
    l1.append(d1.get_sum())
    l2.append(d2.get_sum())

print('v3 : '+str(l1))
print('v4 : '+str(l2))

良くない例
v1 : [39, 39, 39, 39, 39, 39, 39, 39, 39, 39]
v2 : [39, 39, 39, 39, 39, 39, 39, 39, 39, 39]
v3 : [34, 34, 34, 34, 34, 34, 34, 34, 34, 34]
v4 : [33, 33, 33, 33, 33, 33, 33, 33, 33, 33]

print('良い例')

l = []
for __ in range(10):
    d1 = dice2b.Dice()
    for _ in range(10):
        d1.throw()
    l.append(d1.get_sum())

print('v1 : '+str(l))

l = []
for __ in range(10):
    d2 = dice2b.Dice()
    for _ in range(10):
        d2.throw()
    l.append(d1.get_sum())

print('v2 : '+str(l))

l1 = []
l2 = []
for __ in range(10):
    d1 = dice2b.Dice()
    d2 = dice2b.Dice()
    for _ in range(10):
        d1.throw()
        d2.throw()
    l1.append(d1.get_sum())
    l2.append(d2.get_sum())

print('v3 : '+str(l1))
print('v4 : '+str(l2))

良い例
v1 : [35, 30, 33, 35, 33, 36, 41, 35, 35, 36]
v2 : [36, 36, 36, 36, 36, 36, 36, 36, 36, 36]
v3 : [29, 24, 32, 30, 38, 34, 30, 32, 33, 32]
v4 : [38, 31, 37, 37, 43, 34, 30, 35, 50, 32]

数理最適化

線形計画問題

利益を最大化するには、XとYをそれぞれどれだけ製造すれば良いか。

	原料A	原料B	原料C	利益
製品X	1kg	2kg	2kg	3ドル
製品Y	4kg	3kg	1kg	4ドル
在庫	1700kg	1400kg	1000kg	-

以下のように書き換えられる。目的関数も制約式も1次式なので線形計画問題と言われる。

\begin{eqnarray}
&\text{Maximize}\ \ & 3x+4y\\
&\text{Subject to}\ \ & x+4y\leq 1700\\
&& 2x+3y\leq 1400\\
&& 2x+y\leq 1000\\
&& x\geq 0\\
&& y\geq 0
\end{eqnarray}

x = np.linspace(-1,500,100)
y = np.linspace(-1,500,100)

f1 = (1700-x)/4
f2 = (1400-2*x)/3
f3 = (1000-2*x)
plt.plot(x,f1,label='$y=(1700-x)\ /\ 4$',color='red')
plt.plot(x,f2,label='$y=(1400-2x)\ /\ 3$',color='blue')
plt.plot(x,f3,label='$y=1000-2x$',color='purple')

X,Y = np.meshgrid(x,y)
Z = 3*X + 4*Y
plt.contour(X,Y,Z,50)

upper = pd.DataFrame([f1,f2,f3]).min()
bottom = np.full(100,-1)
plt.fill_between(x,upper,bottom,where=upper>bottom,color='gray',alpha=0.3)

plt.colorbar()
plt.legend()
#plt_legend_out()
plt.ylim(0,500)
plt.xlabel('x')
plt.ylabel('y')
plt.scatter(400,200,color='red')
plt.show()

一般化すると以下のようになる。

\begin{eqnarray}
&\text{Minimize}\ & c^Tx\\
&\text{Subject to}\ & Gx\leq h\ (\text{or}\ Ax=b)
\end{eqnarray}

今回のケースを上式に当てはめると、以下のようになる。目的変数の符号を変えている。

\begin{eqnarray}
&\text{Minimize}\ &
  \left(
    \begin{array}{ccc}
      -3 \\
      -4
    \end{array}
  \right)^T
  \left(
    \begin{array}{ccc}
      x  \\
      y 
    \end{array}
  \right)\\
&\text{Subject to}\ &
  \left(
    \begin{array}{ccc}
      1 & 4\\
      2 & 3\\
      2 & 1\\
      -1 & 0\\
      0 & -1
    \end{array}
  \right)
  \left(
    \begin{array}{ccc}
      x\\
      y
    \end{array}
  \right)\leq
    \left(
    \begin{array}{ccc}
      1700\\
      1400\\
      1000\\
      0\\
      0
    \end{array}
  \right)
\end{eqnarray}

import numpy as np
from scipy import optimize

c = np.array([-3, -4], dtype=np.float64)
G = np.array([[1, 4], [2, 3], [2, 1]], dtype=np.float64)
h = np.array([1700, 1400, 1000], np.float64)
sol = optimize.linprog(c, A_ub=G, b_ub=h, bounds=(0, None))

print('(x, y) = '+str(np.round(sol.x,0)))
print('obj.   = '+str(np.round(sol.fun,0)))

(x, y) = [400. 200.]
obj.   = -2000.0

2次計画法

制約なし

以下を解く。

$f(x,y)=x^2+xy+y^2+2x+4y$

cvxoptライブラリでは、以下を標準形とする。

$\frac{1}{2}x^TPx+q^Tx$

以下のように書き換える。

f(x,y)=\frac{1}{2}(x\ y)
\left(
\begin{array}{ccc}
  2 & 1 \\
  1 & 2
\end{array}
\right)
\left(
\begin{array}{ccc}
  x \\
  y
\end{array}
\right)+
(2\ 4)
\left(
\begin{array}{ccc}
  x \\
  y
\end{array}
\right)

import cvxopt

P = cvxopt.matrix(np.array([[2, 1], [1, 2]], dtype=np.float64))
q = cvxopt.matrix(np.array([2, 4], dtype=np.float64))

sol = cvxopt.solvers.qp(P, q)

print('(x, y) = '+str(np.round(np.array(sol["x"]).reshape(-1),1)))
print('obj.   = '+str(np.round(np.array(sol["primal objective"]),1)))

(x, y) = [ 0. -2.]
obj.   = -4.0

制約あり

\begin{eqnarray}
&\text{Minimize}\ & f(x,y)=x^2+xy+y^2+2x+4y\\
&\text{Subject to}\ & x+y=0
\end{eqnarray}

P = cvxopt.matrix(np.array([[2, 1], [1, 2]], dtype=np.float64))
q = cvxopt.matrix(np.array([2, 4], dtype=np.float64))
A = cvxopt.matrix(np.array([[1, 1]], dtype=np.float64))
b = cvxopt.matrix(np.array([0], dtype=np.float64))

sol = cvxopt.solvers.qp(P, q, A=A, b=b)

print('(x,y) = ',end='')
print(np.round(np.array(sol["x"]).reshape(-1),1))
print('obj.  = ',end='')
print(np.round(np.array(sol["primal objective"]),1))

(x,y) = [ 1. -1.]
obj.  = -1.0

#from matplotlib import ticker, cm

x = np.linspace(-5,5,20)
y = np.linspace(-5,5,20)
below = np.full(20,-5)
X, Y = np.meshgrid(x, y)
f1 = X**2 + X*Y + Y**2 + 2*X + 4*Y
f2 = -x

#plt.contourf(X,Y,f1,cmap=cm.PuBu_r)
plt.contour(X,Y,f1,20)
plt.plot(x,f2,color='k')
plt.colorbar()
#plt.fill_between(x, f2, below, where=f2>below, facecolor='gray', alpha=0.3)
#plt.axhline(y=0,color='gray',lw=1)
plt.axvline(x=1,color='k',ls='dashed')
plt.axhline(y=-1,color='k',ls='dashed')
plt.scatter(1,-1,color='red',s=100)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

\begin{eqnarray}
&\text{Minimize}\ & f(x,y)=x^2+xy+y^2+2x+4y\\
&\text{Subject to}\ & 2x+3y\leq 3
\end{eqnarray}

P = cvxopt.matrix(np.array([[2, 1], [1, 2]], dtype=np.float64))
q = cvxopt.matrix(np.array([2, 4], dtype=np.float64))
G = cvxopt.matrix(np.array([[2, 3]], dtype=np.float64))
h = cvxopt.matrix(np.array([3], dtype=np.float64))

sol = cvxopt.solvers.qp(P, q, G=G, h=h)

print()
print('(x,y) = ',end='')
print(np.round(np.array(sol["x"]).reshape(-1),1))
print('obj.  = ',end='')
print(np.round(np.array(sol["primal objective"]),1))

     pcost       dcost       gap    pres   dres
 0:  1.8858e+00  2.9758e-01  2e+00  4e-17  2e+00
 1: -2.1066e+00 -2.1546e+00  5e-02  3e-16  7e-01
 2: -3.9999e+00 -4.0665e+00  7e-02  6e-16  8e-17
 3: -4.0000e+00 -4.0007e+00  7e-04  3e-16  1e-16
 4: -4.0000e+00 -4.0000e+00  7e-06  3e-16  6e-17
 5: -4.0000e+00 -4.0000e+00  7e-08  3e-16  2e-16
Optimal solution found.

(x,y) = [-0. -2.]
obj.  = -4.0

from matplotlib import ticker, cm

x = np.linspace(-5,5,20)
y = np.linspace(-5,5,20)
below = np.full(20,-5)
X, Y = np.meshgrid(x, y)
f1 = X**2 + X*Y + Y**2 + 2*X + 4*Y
f2 = (3 - 2*x) / 3

#plt.contourf(X,Y,f1,cmap=cm.PuBu_r)
plt.contour(X,Y,f1,20)
plt.plot(x,f2,color='k')
plt.colorbar()
plt.fill_between(x, f2, below, where=f2>below, facecolor='gray', alpha=0.3)
plt.axvline(x=0,color='k',ls='dashed')
plt.axhline(y=-2,color='k',ls='dashed')
plt.scatter(0,-2,color='red',s=100)
plt.xlabel('x')
plt.ylabel('y')
plt.show()

勾配降下法

$\text{Minimize}\ \ 5x^2-6xy+3y^2+6x-6y$

アルゴリズム

1. パラメータ$\alpha$,$\epsilon$、初期点$x_k$を定める。
2. $k$を$0$から増やしながら以下を繰り返す。
   1. $||\nabla f(x_k)||\leq \epsilon$ であれば終了する。
   2. $x_{k+1}=x_k-\alpha \nabla f(x_k)$を計算する。

x = np.linspace(-10,10,100)
y = np.linspace(-10,10,100)
X,Y = np.meshgrid(x,y)
Z = 5*X**2 -6*X*Y + 3*Y**2 + 6*X -6*Y
plt.contour(X,Y,Z,50)
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar()
plt.show()

class GradientDescent:
    def __init__(self, f, df, alpha=0.01, eps=1e-6):
        self.f = f
        self.df = df
        self.alpha = alpha
        self.eps = eps
        self.path = None

    def solve(self, init):
        x = init
        path = []
        grad = self.df(x)
        path.append(x)
        while (grad**2).sum() > self.eps**2:
            x = x - self.alpha * grad
            grad = self.df(x)
            path.append(x)
        self.path_ = np.array(path)
        self.x_ = x
        self.opt_ = self.f(x)

def f(xx):
    x = xx[0]
    y = xx[1]
    return 5 * x**2 - 6 * x * y + 3 * y**2 + 6 * x - 6 * y

def df(xx):
    x = xx[0]
    y = xx[1]
    return np.array([10 * x - 6 * y + 6, -6 * x + 6 * y - 6])

algo = GradientDescent(f, df)
initial = np.array([1, 1])
algo.solve(initial)
plt.scatter(initial[0], initial[1], color="red", marker="o",s=50)
plt.plot(algo.path_[:, 0], algo.path_[:, 1], color="k", linewidth=1.5)
plt.scatter(algo.x_[0],algo.x_[1],color='red',s=50)

algo = GradientDescent(f, df)
initial = np.array([1, -1])
algo.solve(initial)
plt.scatter(initial[0], initial[1], color="red", marker="o",s=50)
plt.plot(algo.path_[:, 0], algo.path_[:, 1], color="k", linewidth=1.5)
plt.scatter(algo.x_[0],algo.x_[1],color='red',s=50)

algo = GradientDescent(f, df)
initial = np.array([-1, -1.5])
algo.solve(initial)
plt.scatter(initial[0], initial[1], color="red", marker="o",s=50)
plt.plot(algo.path_[:, 0], algo.path_[:, 1], color="k", linewidth=1.5)
plt.scatter(algo.x_[0],algo.x_[1],color='red',s=50)

algo = GradientDescent(f, df)
initial = np.array([0, -1.5])
algo.solve(initial)
plt.scatter(initial[0], initial[1], color="red", marker="o",s=50)
plt.plot(algo.path_[:, 0], algo.path_[:, 1], color="k", linewidth=1.5)
plt.scatter(algo.x_[0],algo.x_[1],color='red',s=50)

algo = GradientDescent(f, df)
initial = np.array([-1.5, -1.2])
algo.solve(initial)
plt.scatter(initial[0], initial[1], color="red", marker="o",s=50)
plt.plot(algo.path_[:, 0], algo.path_[:, 1], color="k", linewidth=1.5)
plt.scatter(algo.x_[0],algo.x_[1],color='red',s=50)

xs = np.linspace(-2, 2, 300)
ys = np.linspace(-2, 2.300)
xmesh, ymesh = np.meshgrid(xs, ys)
xx = np.r_[xmesh.reshape(1, -1), ymesh.reshape(1, -1)]
plt.contour(xs, ys, f(xx).reshape(xmesh.shape),50)
plt.axvline(x=algo.x_[0],color='k',ls='dashed')
plt.axhline(y=algo.x_[1],color='k',ls='dashed')
plt.colorbar()
plt.xlabel('x')
plt.ylabel('y')
plt.show()

plt.figure(figsize=(8,3))

plt.subplot(121)
algo = GradientDescent(f, df)
initial = np.array([1, 1.3])
algo.solve(initial)
plt.scatter(initial[0], initial[1], color="red", marker="o",s=50)
plt.plot(algo.path_[:, 0], algo.path_[:, 1], color="k", linewidth=1.5)
plt.scatter(algo.x_[0],algo.x_[1],color='red',s=50)

xs = np.linspace(-2, 2, 300)
ys = np.linspace(-2, 2.300)
xmesh, ymesh = np.meshgrid(xs, ys)
xx = np.r_[xmesh.reshape(1, -1), ymesh.reshape(1, -1)]
plt.contour(xs, ys, f(xx).reshape(xmesh.shape),50)
plt.axvline(x=algo.x_[0],color='k',ls='dashed')
plt.axhline(y=algo.x_[1],color='k',ls='dashed')
#plt.colorbar()
plt.xlabel('x')
plt.ylabel('y')
plt.title(r'$\alpha=0.01$')

plt.subplot(122)

algo = GradientDescent(f, df, alpha=0.1)
initial = np.array([1, 1.3])
algo.solve(initial)
plt.scatter(initial[0], initial[1], color="red", marker="o",s=50)
plt.plot(algo.path_[:, 0], algo.path_[:, 1], color="k", linewidth=1.5)
plt.scatter(algo.x_[0],algo.x_[1],color='red',s=50)

xs = np.linspace(-2, 2, 300)
ys = np.linspace(-2, 2.300)
xmesh, ymesh = np.meshgrid(xs, ys)
xx = np.r_[xmesh.reshape(1, -1), ymesh.reshape(1, -1)]
plt.contour(xs, ys, f(xx).reshape(xmesh.shape),50)
plt.axvline(x=algo.x_[0],color='k',ls='dashed')
plt.axhline(y=algo.x_[1],color='k',ls='dashed')
#plt.colorbar()
plt.xlabel('x')
#plt.ylabel('y')
plt.tick_params(left=False, labelleft=False)
plt.title(r'$\alpha=0.1$')

plt.show()

plt.figure(figsize=(8,3))

plt.subplot(121)
algo1 = GradientDescent(f, df)
initial = np.array([1, 1.3])
algo1.solve(initial)
plt.scatter(initial[0], initial[1], color="red", marker="o",s=50)
plt.plot(algo1.path_[:, 0], algo1.path_[:, 1], color="red", linewidth=1.5)
plt.scatter(algo1.x_[0],algo1.x_[1],color='red',s=50)

xs = np.linspace(-2, 2, 300)
ys = np.linspace(-2, 2.300)
xmesh, ymesh = np.meshgrid(xs, ys)
xx = np.r_[xmesh.reshape(1, -1), ymesh.reshape(1, -1)]
plt.contour(xs, ys, f(xx).reshape(xmesh.shape),50)
plt.axvline(x=algo1.x_[0],color='red',ls='dashed')
plt.axhline(y=algo1.x_[1],color='red',ls='dashed')
#plt.colorbar()
plt.xlabel('x')
plt.ylabel('y')
plt.title(r'$\epsilon=10^{-6}$')
plt.xlim(-2,2)
plt.ylim(-2,2)

plt.subplot(122)

algo = GradientDescent(f, df, eps=2)
initial = np.array([1, 1.3])
algo.solve(initial)
plt.scatter(initial[0], initial[1], color="red", marker="o",s=50)
plt.plot(algo.path_[:, 0], algo.path_[:, 1], color="red", linewidth=1.5)
plt.scatter(algo.x_[0],algo.x_[1],color='red',s=50)

xs = np.linspace(-2, 2, 300)
ys = np.linspace(-2, 2.300)
xmesh, ymesh = np.meshgrid(xs, ys)
xx = np.r_[xmesh.reshape(1, -1), ymesh.reshape(1, -1)]
plt.contour(xs, ys, f(xx).reshape(xmesh.shape),50)
plt.axvline(x=algo.x_[0],color='red',ls='dashed')
plt.axhline(y=algo.x_[1],color='red',ls='dashed')
plt.axvline(x=algo1.x_[0],color='red',ls='dashed')
plt.axhline(y=algo1.x_[1],color='red',ls='dashed')
plt.scatter(algo1.x_[0],algo1.x_[1],color='red',s=50)
plt.xlim(-2,2)
plt.ylim(-2,2)
#plt.colorbar()
plt.xlabel('x')
#plt.ylabel('y')
plt.tick_params(left=False, labelleft=False)
plt.title(r'$\epsilon=2$')

plt.show()

ニュートン法

1次元

$f(x)$（$x\in\mathbb{R}$）について、$f(x)=0$を求めたいとする。

アルゴリズム

1. 初期値$x_0$、パラメータ$\epsilon$を定める。
2. 以下を繰り返す。
   1. $y=f(x)$の$x=x_k$における接戦とx軸の交点を$x_{k+1}$とする。$x_{k+1}$を接線の方程式$y=f'(x_k)(x_{k+1}-x_k)+f(x_k)$で$y=0$として求める。
   2. $|x_{k+1}-x_k|\leq \epsilon$となったら終了。

$f(x)=x^3-5x+1=0$

def newton1dim(f, df, x0, eps=1e-10, max_iter=1000):
    l = []
    x = x0
    iter = 0
    while True:
        x_new = x - f(x)/df(x)
        l.append(x)
        if abs(x-x_new) < eps:
            break
        x = x_new
        iter += 1
        if iter == max_iter:
            break
#    return x_new
    return l

def f(x):
    return x**3-5*x+1


def df(x):
    return 3*x**2-5

plt.figure(figsize=(8,3))
plt.subplot(121)

x = np.linspace(-10,10,100)
y = x**3 - 5*x + 1
plt.plot(x,y,color='k',alpha=0.5)
plt.xlabel('x')
plt.ylabel('y')
plt.ylim(-10,10)
plt.xlim(-3,3)
plt.axhline(y=0,color='k')

x_plot = newton1dim(f, df, 10)
y_plot = [f(i) for i in x_plot]
plt.scatter(x_plot,y_plot,color='red',alpha=0.5)
plt.plot(x_plot,y_plot,color='red',alpha=0.5)
#for i in range(len(x_plot)):
#    plt.axvline(x=x_plot[i],color='red',lw=1,alpha=0.5)

x_plot = newton1dim(f, df, -1)
y_plot = [f(i) for i in x_plot]
plt.scatter(x_plot,y_plot,color='orange',alpha=0.5)
plt.plot(x_plot,y_plot,color='orange',alpha=0.5)
#for i in range(len(x_plot)):
#    plt.axvline(x=x_plot[i],color='orange',lw=1,alpha=0.5)

x_plot = newton1dim(f, df, -0.5)
y_plot = [f(i) for i in x_plot]
plt.scatter(x_plot,y_plot,color='blue',alpha=0.5)
plt.plot(x_plot,y_plot,color='blue',alpha=0.5)

#for i in range(len(x_plot)):
#    plt.axvline(x=x_plot[i],color='blue',lw=1,alpha=0.5)

x_plot = newton1dim(f, df, -10)
y_plot = [f(i) for i in x_plot]
plt.scatter(x_plot,y_plot,color='green',alpha=0.5)
plt.plot(x_plot,y_plot,color='green',alpha=0.5)
#for i in range(len(x_plot)):
#    plt.axvline(x=x_plot[i],color='green',lw=1,alpha=0.5)
plt.title('$\epsilon =10^{-10}$')    
    
plt.subplot(122)    
x = np.linspace(-10,10,100)
y = x**3 - 5*x + 1
plt.plot(x,y,color='k',alpha=0.5)
plt.xlabel('x')
plt.ylim(-10,10)
plt.xlim(-3,3)
plt.axhline(y=0,color='k')
plt.tick_params(left=False,labelleft=False)
x_plot = newton1dim(f, df, 10, eps=0.1)
y_plot = [f(i) for i in x_plot]
plt.scatter(x_plot,y_plot,color='red',alpha=0.5)
plt.plot(x_plot,y_plot,color='red',alpha=0.5)
#for i in range(len(x_plot)):
#    plt.axvline(x=x_plot[i],color='red',lw=1,alpha=0.5)

x_plot = newton1dim(f, df, -0.5, eps=0.1)
y_plot = [f(i) for i in x_plot]
plt.scatter(x_plot,y_plot,color='blue',alpha=0.5)
plt.plot(x_plot,y_plot,color='blue',alpha=0.5)
#for i in range(len(x_plot)):
#    plt.axvline(x=x_plot[i],color='blue',lw=1,alpha=0.5)

x_plot = newton1dim(f, df, -1, eps=0.1)
y_plot = [f(i) for i in x_plot]
plt.scatter(x_plot,y_plot,color='orange',alpha=0.5)
plt.plot(x_plot,y_plot,color='orange',alpha=0.5)

x_plot = newton1dim(f, df, -10, eps=0.1)
y_plot = [f(i) for i in x_plot]
plt.scatter(x_plot,y_plot,color='green',alpha=0.5)
plt.plot(x_plot,y_plot,color='green',alpha=0.5)
#for i in range(len(x_plot)):
#    plt.axvline(x=x_plot[i],color='green',lw=1,alpha=0.5)

plt.title('$\epsilon =0.1$')
plt.show()

2次元

解きたい問題

f(x)=
\left(
\begin{array}{ccc}
  f_1(x) \\
  f_2(x)
\end{array}
\right)=
\left(
\begin{array}{ccc}
  x^3-2y \\
  x^2+y^2-1
\end{array}
\right)=
\left(
\begin{array}{ccc}
  0 \\
  0
\end{array}
\right)

更新式。考え方は1次元の場合と同じ。

$\ \ \ \ x_{k+1}=x_k-\color{red}{J_f(x_k)}^{-1}f(x_k)$

収束条件

$\ \ \ \ ||x_{k+1}-x_k||\leq \epsilon$

ヤコビアン行列

\ \ \ \ \color{red}{J_f(x)}=
\left(
\begin{array}{ccc}
  \frac{\partial f_1}{\partial x}(x) & \frac{\partial f_1}{\partial y}(x)\\
  \frac{\partial f_2}{\partial x}(x) & \frac{\partial f_2}{\partial y}(x)
\end{array}
\right)=
\left(
\begin{array}{ccc}
  x^3-2y \\
  x^2+y^2-1
\end{array}
\right)=
\left(
\begin{array}{ccc}
  0 \\
  0
\end{array}
\right)

from numpy import linalg

class Newton:
    def __init__(self, f, df, eps=1e-10, max_iter=1000):
        self.f = f
        self.df = df
        self.eps = eps
        self.max_iter = max_iter

    def solve(self, x0):
        x = x0
        iter = 0
        self.path_ = x0.reshape(1, -1)
        while True:
            x_new = x - np.dot(linalg.inv(self.df(x)), self.f(x))
            self.path_ = np.r_[self.path_, x_new.reshape(1, -1)]
            if ((x-x_new)**2).sum() < self.eps*self.eps:
                break
            x = x_new
            iter += 1
            if iter == self.max_iter:
                break
        return x_new

def f1(x, y):
    return x**3-2*y

def f2(x, y):
    return x**2+y**2-1

def f(xx):
    x = xx[0]
    y = xx[1]
    return np.array([f1(x, y), f2(x, y)])

def df(xx):
    x = xx[0]
    y = xx[1]
    return np.array([[3*x**2, -2], [2*x, 2*y]])

plt.figure(figsize=(8,3))
xmin, xmax, ymin, ymax = -3, 3, -3, 3

plt.subplot(121)
x = np.linspace(xmin, xmax, 200)
y = np.linspace(ymin, ymax, 200)
xmesh, ymesh = np.meshgrid(x, y)
z1 = f1(xmesh, ymesh)
z2 = f2(xmesh, ymesh)
plt.contour(xmesh, ymesh, z1, colors="k", levels=[0])
plt.contour(xmesh, ymesh, z2, colors="K", levels=[0])
#solver = newton.Newton(f, df)
solver = Newton(f, df)

initials = [np.array([1, 1]),
            np.array([-2, 2]),
            np.array([-1, -1]),
            np.array([1, -1])]
#markers = ["+", "*", "x"]

c = sns.color_palette('husl',len(initials))
i = 0
#for x0, m in zip(initials, markers):
for x0 in initials:
    sol = solver.solve(x0)
#    plt.scatter(solver.path_[:, 0], solver.path_[:, 1], color="k", marker=m)
    plt.scatter(solver.path_[:, 0], solver.path_[:, 1], color=c[i], alpha=0.8)
    plt.plot(solver.path_[:, 0], solver.path_[:, 1], color=c[i], alpha=0.8)
    i += 1
#    print(sol)
plt.xlabel('x')
plt.ylabel('y')
xmin, xmax, ymin, ymax = -3, 3, -3, 3
plt.xlim(xmin, xmax)
plt.ylim(ymin, ymax)
plt.title('$\epsilon =10^{-10}$')

##########

plt.subplot(122)
x = np.linspace(xmin, xmax, 200)
y = np.linspace(ymin, ymax, 200)
xmesh, ymesh = np.meshgrid(x, y)
z1 = f1(xmesh, ymesh)
z2 = f2(xmesh, ymesh)
plt.contour(xmesh, ymesh, z1, colors="k", levels=[0])
plt.contour(xmesh, ymesh, z2, colors="K", levels=[0])
#solver = newton.Newton(f, df)
solver = Newton(f, df, eps=2)

initials = [np.array([1, 1]),
            np.array([-2, 2]),
            np.array([-1, -1]),
            np.array([1, -1])]
#markers = ["+", "*", "x"]

c = sns.color_palette('husl',len(initials))
i = 0
#for x0, m in zip(initials, markers):
for x0 in initials:
    sol = solver.solve(x0)
#    plt.scatter(solver.path_[:, 0], solver.path_[:, 1], color="k", marker=m)
    plt.scatter(solver.path_[:, 0], solver.path_[:, 1], color=c[i], alpha=0.8)
    plt.plot(solver.path_[:, 0], solver.path_[:, 1], color=c[i], alpha=0.8)
    i += 1
#    print(sol)
plt.xlabel('x')
plt.tick_params(left=False,labelleft=False)
#plt.ylabel('y')
plt.xlim(xmin, xmax)
plt.ylim(ymin, ymax)
plt.title('$\epsilon =2$')
plt.show()

/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:11: MatplotlibDeprecationWarning: Support for uppercase single-letter colors is deprecated since Matplotlib 3.1 and will be removed in 3.3; please use lowercase instead.
  # This is added back by InteractiveShellApp.init_path()

ラグランジュ未定乗数法

\begin{eqnarray}
&\text{Minimize}\ \ & f(x,y)=5x^2+6xy+5y^2-26x-26y\\
&\text{Subject to}\ \ & g(x,y)=x^2+y^2-4=0
\end{eqnarray}

以下は$f$が等高線（破線）、$g$が制約条件（実線）。$f$は最小値の方向に進んでいる。

（省略）

考え方

$\nabla f=-\lambda\nabla g$となる$\lambda$が存在する。これを式変形すると$\nabla f+\lambda\nabla g=0$となる。

新しく$\nabla L$を定義。$L$を最大化する？

\begin{eqnarray}
\nabla L &=& \nabla f+\lambda \nabla g = 0
\end{eqnarray}

統計

# 入手元：気象庁ホームページ
# （https://www.data.jma.go.jp/gmd/risk/obsdl/index.php）
# 2018年4月の東京の最高気温（日別）
x = np.array([21.9, 24.5, 23.4, 26.2, 15.3, 22.4, 21.8, 16.8,
              19.9, 19.1, 21.9, 25.9, 20.9, 18.8, 22.1, 20.0,
              15.0, 16.0, 22.2, 26.4, 26.0, 28.3, 18.7, 21.3,
              22.5, 25.0, 22.0, 26.1, 25.6, 25.7])

m = x.sum() / len(x)
s = np.sqrt(((x - m)**2).sum() / len(x))
print("平均値:{:.4f}".format(m))
print("標準偏差:{:.4f}".format(s))

平均値:22.0567
標準偏差:3.4908

# 入手元：気象庁ホームページ（https://www.data.jma.go.jp/gmd/risk/obsdl/index.php）
# 2018年4月の東京の最高気温（日別）
x = np.array([21.9, 24.5, 23.4, 26.2, 15.3, 22.4, 21.8, 16.8,
              19.9, 19.1, 21.9, 25.9, 20.9, 18.8, 22.1, 20.0,
              15.0, 16.0, 22.2, 26.4, 26.0, 28.3, 18.7, 21.3,
              22.5, 25.0, 22.0, 26.1, 25.6, 25.7])
# 2018年4月の札幌の最高気温（日別）
y = np.array([8.3, 13.0, 8.4, 7.9, 7.0, 3.7, 6.1, 8.5, 8.6,
              11.9, 12.1, 14.4, 7.0, 10.5, 6.6, 10.6, 16.6,
              19.1, 20.1, 19.8, 24.5, 12.6, 16.4, 13.0, 13.3,
              14.1, 14.4, 17.0, 21.3, 24.5])


mx = x.sum() / len(x)
my = y.sum() / len(y)
sx = np.sqrt(((x - mx)**2).sum() / len(x))
sy = np.sqrt(((y - my)**2).sum() / len(y))
sxy = ((x - mx) * (y - my)).sum() / len(x)
print("東京の最高気温の標準偏差:{:4f}".format(sx))
print("札幌の最高気温の標準偏差:{:4f}".format(sy))
print("共分散:{:4f}".format(sxy))
print("相関係数:{:4f}".format(sxy / (sx * sy)))

東京の最高気温の標準偏差:3.490815
札幌の最高気温の標準偏差:5.425414
共分散:5.487211
相関係数:0.289729

plt.scatter(x,y)
plt.xlabel('東京の気温')
plt.ylabel('札幌の気温')
plt.show()