数学とPythonのメモ帳

僕のメモ帳です。同じ疑問を感じた方の参考になれば良いです。

サポートベクターマシンのパーセプトロンアルゴリズムをPythonで

2クラス分類を解く一番簡単なアルゴリズムです。

#適当な10点のx1,x2成分
x1 = np.array([0,2,4,6,8,10,12,14,16,18])
x2 = np.array([1,2,1,4,3,9,7,6,7,11])

#それぞれのラベル
y = np.array([-1,-1,-1,-1,-1,1,1,1,1,1])

#点を全てプロット
plt.scatter(px, py)

#重みは0に近い値で初期化
w = np.array([0.001,0.001,0.001])

#失敗の数を表す
k = 0

#サンプルの中での、ノルムの最大値
R = np.linalg.norm(np.array([x1[9], x2[9]]))

#学習率 適当に小さく
eta = 0.001

#繰り返し継続判定のフラグ
miss = True

#パーセプトロンの学習
while(miss):
    miss = False
    for i in range(10):
        if(y[i] * (w[0]*x1[i] + w[1]*x2[i] + w[2]) <= 0): #分類ミス
            miss = True
            w[0] = w[0] + eta * y[i] * x1[i]
            w[1] = w[1] + eta * y[i] * x2[i]
            w[2] = w[2] + eta * y[i] * R * R
            k += 1

#-3から20まで0.1刻みでnumpy配列で用意
x = np.arange(-3, 20, 0.1)

#学習して得た直線の式
y = -(w[0] / w[1])*x - (w[2] / w[1])

plt.plot(x, y)

plt.show()

f:id:keitainoue157:20181214010343p:plain

いい感じに分類してくれていますね!

最小二乗法をPythonで実装

最小二乗法では、重みはこうなります。

\begin{align}
w\bf = (X^TX)^{-1}X^Ty\bf
\end{align}

#適当な6点のx,y座標
px = np.array([0,1,2,3,4,5])
py = np.array([1,5,6,8,9,9])

e = np.array([1,1,1,1,1,1]) #成分が全て1のベクトル
X = np.array([px, e]).T #6点のx成分と1を縦に並べた行列(最小二乗法に必要)

#最小二乗法による重みの解
w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(py)

#点を全てプロット
plt.scatter(px, py)

#-3から7まで0.1刻みでnumpy配列で用意
x = np.arange(-3, 7, 0.1)

#得られた重みを用いて、直線の式を作る
y = w[0]*x + w[1]

#得られた直線を引く
plt.plot(x, y)
plt.show()

f:id:keitainoue157:20181214002132p:plain

ちなみに、Windrow-Hoffアルゴリズムという最小二乗法を解く手法も紹介する

#適当な6点のx,y座標
px = np.array([0,1,2,3,4,5])
py = np.array([1,5,6,8,9,9])

#重みは0で初期化
w = np.array([0, 0])

#学習率 適当に小さく
eta = 0.001

#Widrow-Hoffアルゴリズムで最小二乗法を解く
for _ in range(1000):
    for i in range(6):
        w = w - eta*(w[0]*px[i] + w[1] - py[i])*np.array([px[i], 1])

#点を全てプロット
plt.scatter(px, py)

#-3から7まで0.1刻みでnumpy配列で用意
x = np.arange(-3, 7, 0.1)

#得られた重みを用いて、直線の式を作る
y = w[0]*x + w[1]

#得られた直線を引く
plt.plot(x, y)
plt.show()

f:id:keitainoue157:20181214003948p:plain

点に沿ったいい感じの直線が引けてて気持ちいいね!

線形計画問題


線形計画問題の解法の一つである「シンプレックス法」をPythonで実装してみました。

def LP(c, A, b, xb_idx):
    xn_idx = np.arange(A.shape[1])
    xn_idx = np.delete(xn_idx, xb_idx)

    while True:
        B = A[:, xb_idx]
        N = A[:, xn_idx]
        B_inv = np.linalg.inv(B) # Bの逆行列
        b_bar = np.dot(B_inv, b) # bバー。(xB, xN) = (bバー, 0)
        if np.any(b_bar < 0) :
            print('initial basic solution isn\'t solution!')
            return 0
        cB = c[xb_idx]
        cN = c[xn_idx]
        pi = np.dot(B_inv.T, cB) # シンプレックス乗数
        rcf = cN - np.dot(N.T, pi) # 相対コスト係数
        
        if np.all(rcf >= 0) :
            solution = np.zeros_like(c, dtype=float)
            solution[xb_idx] = b_bar
            return solution
        
        for i in range(len(rcf)):
            if rcf[i] < 0 :
                choose = xn_idx[i]
                choose_idx = i
                break
            
        y = np.dot(B_inv, A[:, choose])
        y = y.reshape(len(y), 1)
        
        if np.all(y <= 0) :
            print('Unbounded')
            return 0
        
        theta_opt = b_bar / y
        
        theta = np.min(theta_opt)
        dropout = np.argmin(theta_opt)
        
        temp = xn_idx[choose_idx]
        xn_idx[choose_idx] = xb_idx[dropout]
        xb_idx[dropout] = temp
        xb_idx = sorted(xb_idx)
        xn_idx = sorted(xn_idx)
具体的にはこう使います。

\begin{align}
&minimize &-2x_1-3x_2 \\
&\text{subject to} &4x_1+3x_2+x_3&=120\\
&&x_1+2x_2+x_4&=60\\
&&x_i \geq 0 &&\forall{i = 1, 2, 3, 4}
\end{align}

こんな時は

A = np.array([[4,3,1,0],[1,2,0,1]])
b = np.array([[120,60]]).T
c = np.array([[-2,-3,0,0]]).T
xb_idx = [2,3]
x = LP(c, A, b, xb_idx)
print('Optimum solution')
print(x)
print('objective function -> ',np.round(np.dot(c.T, x)[0][0]))
print('constraints -> ', np.all(np.round(np.dot(A, x) - b, 7) == 0))
出力は


Optimum solution
[[12.]
[24.]
[ 0.]
[ 0.]]
objective function -> -96.0
constraints -> True

matplotlibで陰関数のグラフを描画してみた

陰関数とは、 f(x, y) = 0というものです。

Pythonのmatplotlibを使って、 \boldsymbol{L_\inftyノルム}を、陰関数に直して図示してみました。

まずは  \boldsymbol{L_pノルム}の定義から。

 \boldsymbol{R}ⁿ の元  x = (x_1, x_2, \ldots, x_n) に対して、その \boldsymbol{Lpノルム} \|x\|_p
\begin{align}
\|x\|_p = \sqrt[p]{|x_1|^p + |x_2|^p + \cdots + |x_n|^p}
\end{align}と定義する。

さて、定義から

\begin{align}
\|x\|_1 = |x_1| + |x_2| + \cdots + |x_n|
\end{align}

\begin{align}
\|x\|_2 = \sqrt{|x_1|² + |x_2|² + \cdots + |x_n|²}
\end{align}

ですね。

 x = (x_1, x_2), \|x\|_p = 1とし、 p = 1, 2, \infty(200)の場合で図示してみましょう。

この場合、

import matplotlib.pyplot as plt
import numpy as np

#Lpノルムの p
p = 200

x1range = np.arange(-2, 2, 0.025)
x2range = np.arange(-2, 2, 0.025)
X1, X2 = np.meshgrid(x1range, x2range)

plt.axis([-2, 2, -2, 2])
plt.gca().set_aspect('equal', adjustable='box')

Z = abs(X1)**p + abs(X2)**p - 1
plt.contour(X1, X2, Z, [0])
plt.show()

f:id:keitainoue157:20180905035617j:plain
 L_1ノルム
f:id:keitainoue157:20180905035754j:plain
 L_2ノルム
f:id:keitainoue157:20180905040030j:plain
 L_{200}ノルム

 p = 200 で、正方形っぽくなりましたね。

実は、 p \to \inftyで完全な正方形になります。

\begin{align}
\|x\|_\infty = max \left\{|x_1|, |x_2|, \cdots \right\}
\end{align}

 \|x\|_\infty = 1 の時、
 |x_i| < 0 なら、  |x_i|^\infty = 0
 |x_i| > 0 なら、  |x_i|^\infty = \infty
となるため、
\begin{align}
\|x\|_\infty = 1 \longrightarrow |x_i| \leq 1 (i = 1, 2, \cdots, n)
\end{align}
となり、このようなグラフになるわけですね。

ではまた。