Kei Minagawa's Blog

皆川圭(@keimina)のブログ、Pythonで試したことを書いていきます

点を移動させる行列と、その固有ベクトルについて可視化して理解を深める

1. はじめに

2x2の行列を点(x, y)に適応すると、新しい点に移動します。本記事は、この行列による点の変換を可視化し行列の理解を深める記録です。本題ではありませんが、行列の固有ベクトルも表示し、固有ベクトルについても理解を深めました。

2. 固有ベクトルについて

固有ベクトルは本記事の本題ではありませんが、固有ベクトルとは何なのかというと、行列を適応しても定数倍の効果しか及ぼさないベクトルをその行列の固有ベクトルと呼ぶようです。詳細は以下を参照してください。

ja.wikipedia.org

3. 仕様

行列による点の変換を可視化しするために考えた仕様は以下になります。

  1. 行列適用前の格子状の点を表示する(グレー)
  2. 行列適用後の点を表示する(赤)
  3. 固有ベクトルを表示する(青)
  4. 固有ベクトルが存在しないとき(固有ベクトルの要素に複素数が含まれるとき)は固有ベクトルを表示しない
  5. 点の移動量を矢印で表示する(グレー)
  6. 変換行列をスライダで調整できるようにする

4. ソースコード

仕様に従い実装しました。ソースコードは以下の通りです。

# xy座標上の点 A(x1, y1) を点B(x2, y2)に移動(写像)する行列 M を考えます。
# 点B は行列 M と 点A を用いて以下のように表せます。
# (a, b は点A と点B を表す縦ベクトルとします)
# b = M・a
import numpy as np
from sympy import Matrix, MatrixSymbol

# 描画関連
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button
from matplotlib.patches import Arrow

# デバッグ用
# plt.ion()

# 行列適用前の点
xs, ys = np.mgrid[-10:11:2, -10:11:2] / 10

# 行列演算の定義
N = 2
M = MatrixSymbol('M', N, N)
a = MatrixSymbol('a', N, 1)
b = M*a

# 行列適用後の点
def calc_point(mat, point):
    return b.subs({M:mat, a:point}).doit()

# スライダーの axis を作成
fig, (ax, ax_button, *ax_sldrs) = plt.subplots(6, 1, gridspec_kw = {'height_ratios':[28, 1, 1, 1, 1, 1]})

# スライダーの名前を定義
snames = ['p11', 'p12', 'p21', 'p22']

# axis にボタンを作成
button = Button(ax_button, 'show/hide vector')

# axis にスライダーを作成
sldrs = [Slider(ax_sldr, sname, -1.0, 1.0, valinit=0.0) for ax_sldr, sname in zip(ax_sldrs, snames)]

# 行列 M は
# M = [[sldrs[0].val, sldrs[1].val],
#      [sldrs[2].val, sldrs[3].val]]
# とする

# 矢印の表示の有無
patches_show_state = True

def init_sldrs():
    # スライダの初期値を設定する
    # 単位行列
    sldrs[0].set_val(1.0)
    sldrs[1].set_val(0.0)
    sldrs[2].set_val(0.0)
    sldrs[3].set_val(1.0)

def calc_points():
    points = []
    # スライダから値を取得
    p11, p12 = sldrs[0].val, sldrs[1].val
    p21, p22 = sldrs[2].val, sldrs[3].val
    for x, y in zip(xs.flatten(), ys.flatten()):
        mat = Matrix([[p11, p12], [p21, p22]])
        # 行列適用前の点の座標を取得
        point1 = Matrix([[x], [y]])
        # 行列を適応
        point2 = calc_point(mat, point1)
        # 行列適用後の座標をリストに保存
        points.append(point2)
    # 保存したリストを x座標のリスト、 y座標のリストに変換する
    x, y = list(zip(*points))
    return x, y

def init_ax():
    # 行列適用後の点を描画する
    x, y = calc_points()
    ax_line, = ax.plot(x, y, color='r', marker='.', linestyle='None')
    # オリジナルの点も描画する
    ax.plot(xs, ys, color='k', alpha=0.1, marker='.', linestyle='None')
    # 矢印も描画する
    for x1, y1, x2, y2 in zip(xs.flatten(), ys.flatten(), x, y):
        # Arrow(x, y, dx, dy) で dx と dy を float にしないとエラーになる
        ax.add_patch(Arrow(x1, y1, float(x2-x1), float(y2-y1)))
    # 虚数が含まれない場合、固有ベクトルを表示する
    p11, p12 = sldrs[0].val, sldrs[1].val
    p21, p22 = sldrs[2].val, sldrs[3].val
    w, v = np.linalg.eig([[p11, p12],[p21, p22]])
    if np.all(np.isreal(v)):
        ax.add_patch(Arrow(0, 0, *(2*v[:,0]),  linestyle='solid', color='blue', width=0.02))
        ax.add_patch(Arrow(0, 0, *(2*v[:,1]),  linestyle='solid', color='blue', width=0.02))

    return ax_line

def update_ax(ax_line):
    # 行列適用後の点を描画する
    x, y = calc_points()
    ax_line.set_xdata(x)
    ax_line.set_ydata(y)
    update_patches()

def update_patches():
    x, y = calc_points()
    # 矢印を更新
    # https://stackoverflow.com/questions/49791848/matplotlib-remove-all-patches-from-figure?noredirect=1&lq=1
    for patch in reversed(ax.patches):
        patch.remove()
    if patches_show_state:
        for x1, y1, x2, y2 in zip(xs.flatten(), ys.flatten(), x, y):
            # Arrow(x, y, dx, dy) で dx と dy を float にしないとエラーになる
            ax.add_patch(Arrow(x1, y1, float(x2-x1), float(y2-y1), linestyle='solid', color='grey', width=0.02))
        # 固有ベクトルを求める
        p11, p12 = sldrs[0].val, sldrs[1].val
        p21, p22 = sldrs[2].val, sldrs[3].val
        w, v = np.linalg.eig([[p11, p12],[p21, p22]])
        # 虚数が含まれない場合、固有ベクトルを表示する
        if np.all(np.isreal(v)):
            ax.add_patch(Arrow(0, 0, *(2*v[:,0]),  linestyle='solid', color='blue', width=0.02))
            ax.add_patch(Arrow(0, 0, *(2*v[:,1]),  linestyle='solid', color='blue', width=0.02))
    fig.canvas.draw()

def button_clicked(x):
    global patches_show_state
    patches_show_state = not patches_show_state
    update_patches()

# 描画範囲を固定する
ax.set_xlim(-2.0, 2.0)
ax.set_ylim(-2.0, 2.0)

# x軸のラベルを上に持ってくる
# https://stackoverflow.com/questions/14406214/moving-x-axis-to-the-top-of-a-plot-in-matplotlib
ax.xaxis.tick_top()

# アスペクト比を 1:1 にする
ax.set_aspect("equal")

# スライダを初期化
init_sldrs()

# グラフを初期化
ax_line = init_ax()

# スライダが変わった時に呼ばれる関数を登録
for sldr in sldrs:
    sldr.on_changed(lambda t: update_ax(ax_line))

# ボタンがクリックした時に呼ばれる関数を登録
button.on_clicked(button_clicked)

plt.show()

5. ソースコードの説明

calc_point(mat, point) 関数により、与えられた行列を点(ベクトル)の計算を行っています。行列とベクトルの計算の定義には sympy.MatrixSymbol を使用しています。calc_point 関数にすべての格子状の点に適応し、格子状の点を移動させています。スライダを動かすと update_ax(ax_line) 関数がよばれ、これらの計算処理と矢印の描画などが行われ、表示が更新されます。


固有ベクトルを計算するのに np.linalg.eigh を使用しました。
docs.scipy.org

矢印の表示には matplotlib.patches.Arrow を使用しました。矢印の再描画ではまったのですが、矢印を描画する前に矢印をすべて消去(remove)して再描画しています。矢印の消去には以下を参考にしました。

stackoverflow.com

6. 出力結果

出力されるグラフは以下になります。スライダで2x2行列の各要素の値を変更し変換行列の作用を見ることができます。

f:id:keimina:20190429125528g:plain

点や矢印の意味は以下の通りです。

No. 点や矢印 意味
1 赤の点 移動後の点
2 グレーの点 変換前の点
3 グレーの矢印 変換前と変換後の点を結んだ矢印
4 青の矢印 行列の固有ベクトル

7. 考察

行列の初期状態は単位行列ですので、この状態では変換前と変換後で点の位置はかわりません。行列は[[p11, p12], [p21, p22]] で p11, p12, p21, p22 のパラメータからなります。 p11 を動かすと x 方向に変形します。 p22 を動かすと y 方向に変形します。 p12 を動かすと平行四辺形のように x 平方向に変形します。 p21 を動かすと平行四辺形のように y 方向に変形します。p12 と p21 の符号が違うとき 点は回転するように変形します。そしてこのとき固有ベクトルの要素が複素数になります(そのため固有ベクトルの矢印は表示していません)。点の移動量をみると、あたかも固有ベクトルの方向を目指すように点が移動しているように見えます。このように固有ベクトルは、行列の作用を説明してくれているようにも見えます。

8. まとめ

2x2 行列による点の変換と、固有ベクトルについて可視化しました。スライダーを動かして変換行列と、それに対応する固有ベクトルの性質を見ることができました。

以上です。