行列の固有値分解（スペクトル分解）の写像を可視化

行列によるベクトルの写像は，行列の固有値分解で3段階で写像されている．それを可視化してみる．

固有値分解

固有値分解の定義は

$\mathbf{A} \vec{a} = \lambda \vec{a}$

固有方程式は

$|\mathbf{A} - \lambda \mathbf{I}| = 0$

行列Aをp次正方行列として固有値$\lambda_{1}, \cdots , \lambda_{p}$と固有ベクトル$\vec{a}_{1}, \cdots, \vec{a}_{p}$として，行列を分解する．

$\begin{eqnarray} \mathbf{A} \begin{pmatrix} \vec{a}_{1} &\cdots & \vec{a}_{p} \end{pmatrix} &=& \begin{pmatrix} \vec{a}_{1} & \cdots & \vec{a}_{p} \end{pmatrix} \begin{pmatrix} \lambda_{1} & 0 & \cdots & 0\\ 0 &\lambda_{2} & \cdots & 0\\ 0 & 0 &\cdots & \lambda_{p} \end{pmatrix}\\\\ \mathbf{A}\mathbf{P} &=& \mathbf{P} \mathbf{\Lambda} \\\\ \mathbf{A} &=& \mathbf{P} \mathbf{\Lambda} \mathbf{P}^{-1} \\\\ \end{eqnarray}$

対称行列の固有値分解の性質

特に実対称行列の固有値分解では，以下2つの性質がある．

固有値が必ず実数
異なる固有値に対応する固有ベクトルは直交する（同じ固有値でもグラムシュミットの直交基底化で直行するベクトルが作れる）（対称行列でない場合，直交しないので斜交になる．）

対称行列の固有値と固有ベクトルの性質の証明 | 高校数学の美しい物語

応用例：

固有ベクトルを並べた行列は直交行列になるので，ベクトルを写像すると直交座標系の変換になる．
特に統計では共分散・相関係数行列が対称行列になり，それを固有値分解するとそのデータをよく表現する基底ベクトルを作る直交座標系の変換になり主成分分析(PCA), 因子分析(FA)に繋がる．

固有値分解で見える行列の写像

$\begin{eqnarray} \mathbf{A} = \mathbf{P} \mathbf{\Lambda} \mathbf{P}^{-1} \end{eqnarray}$

より行列は3段階に分けられている． $\vec{x}$を行列$\mathbf{A}$で写像してこれを順番に見ていく．

まず行列を決める．

$\begin{eqnarray} \mathbf{A} = \begin{pmatrix} 8 & -2 \\ 3 & 1 \end{pmatrix} \end{eqnarray}$

ベクトル$\vec{x} = (0, 0)^{T}$から$(3, 3)^{T}$までの写像を描画してみる．

import scipy as sp
from scipy import linalg
import matplotlib.pyplot as plt

A = sp.matrix([[8, -2],
               [3, 1]])

for x in sp.linspace(0, 3, 10):
    for y in sp.linspace(0, 3, 10):
        v = sp.array([x, y])
        map_v = sp.asarray(A.dot(v)).reshape(-1)
        plt.scatter(x, y, color="k", s=1)
        plt.quiver(0, 0, *map_v, color="r",angles='xy', scale_units='xy', scale=1)


plt.quiver(0, 0, *(16*sp.array([0.31622777, 0.9486833])), angles='xy', scale_units='xy', scale=1)
plt.quiver(0, 0, *(16*sp.array([0.89442719, 0.4472136])), angles='xy', scale_units='xy', scale=1)
plt.grid()
plt.xlim(-10, 25)
plt.ylim(-5, 20)
plt.show()

固有ベクトル上のベクトルであれば同じ向きに，他のベクトルは違う向きに写像されていることがわかる．

次にベクトル1つ$\vec{x}=(1, 1)^{T}$に対してどのような変換が行われていくか段階ごとに見る．

1段階: $P^{-1}\vec{x}$

標準基底から固有ベクトルへの基底の取替え（固有ベクトル基底を恒等変換して標準基底に関する表現行列）を考える．座標系変換後の(1,1)の座標を$(\alpha_{1},\alpha_{2} )^{T}$とする．

$\begin{eqnarray} \begin{pmatrix} \vec{a}_{1} & \vec{a}_{2} \end{pmatrix} \begin{pmatrix} \alpha_{1} \\ \alpha_{2} \end{pmatrix} &=& \begin{pmatrix} \vec{e}_{1} & \vec{e}_{2} \end{pmatrix}\vec{x} \\\\ \begin{pmatrix} \alpha_{1} \\ \alpha_{2} \end{pmatrix}_{\left< \vec{a}_{1}, \vec{a}_{2} \right> } &=& \begin{pmatrix} \vec{a}_{1} & \vec{a}_{2} \end{pmatrix}^{-1}\mathbf{I} \vec{x} \\\\ &=& \mathbf{P}^{-1}\vec{x} \end{eqnarray}$

つまり，固有値分解の右側は，標準基底から固有ベクトルを新たな基底とする基底の取替（座標系の変換）を行っている．元のベクトル自体は変わっていない．

（斜交座標系になるので直角でなく平行四辺形になる．）

2段階：$\mathbf{\Lambda} \mathbf{P}^{-1}\vec{x} $

固有値を並べきた行列$\Lambda$は対角行列より，固有ベクトルに沿って（固有ベクトル基底の各成分を）固有値でスケール倍することになる．ここで線形写像される．

$\begin{eqnarray} \begin{pmatrix} \alpha_{1}' \\ \alpha_{2}' \end{pmatrix}_{\left< \vec{a}_{1}, \vec{a}_{2} \right> } &=& \mathbf{\Lambda} \mathbf{P}^{-1}\vec{x} \\\\ &=& \mathbf{\Lambda} \begin{pmatrix} \alpha_{1} \\ \alpha_{2} \end{pmatrix}_{\left< \vec{a}_{1}, \vec{a}_{2} \right> } \\\\ &=& \begin{pmatrix} \lambda_{1} \alpha_{1} \\ \lambda_{2} \alpha_{2} \end{pmatrix}_{\left< \vec{a}_{1}, \vec{a}_{2} \right> } \\\\ \end{eqnarray}$

3段階：$\mathbf{P} \mathbf{\Lambda} \mathbf{P}^{-1}\vec{x} $

最後に，固有ベクトル基底から標準基底に戻している．

$\begin{eqnarray} \vec{x}' &=& \mathbf{P} \begin{pmatrix} \alpha_{1}' \\ \alpha_{2}' \end{pmatrix}_{\left< \vec{a}_{1}, \vec{a}_{2} \right> } &=& \mathbf{P} \mathbf{\Lambda} \mathbf{P}^{-1}\vec{x} \\\\ \end{eqnarray}$

最後に図をまとめる

import scipy as sp
from scipy import linalg
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

%matplotlib inline

fig = plt.figure(figsize=(5*3, 5))
gs = gridspec.GridSpec(nrows=1, ncols=3)

# 標準基底
def view_basic_vec(ax):
    e_vec_prop = {
        "linestyle": "--",
        "length_includes_head": True,
        "head_width": .05,
        "width": .008,
        "color": "k",
    }
    e1 = sp.array([1, 0])
    e2 = sp.array([0, 1])
    ax.arrow(0, 0, *e1, **e_vec_prop)
    ax.arrow(0, 0, *e2, **e_vec_prop)
    ax.text(*e1,"$\\vec{e}_{1}$")
    ax.text(*e2,"$\\vec{e}_{2}$")
    
# 固有ベクトル
def view_eig_vec(ax):
    vec_a1 = sp.asarray(eig_vecs[:,0]).reshape(-1)
    vec_a2 = sp.asarray(eig_vecs[:,1]).reshape(-1)
    ax.text(*vec_a1,"$\\vec{a}_{1}$", color="r")
    ax.text(*vec_a2,"$\\vec{a}_{2}$", color="r")
    ax.quiver(0, 0, *vec_a1, angles='xy', scale_units='xy', scale=1, color="r")
    ax.quiver(0, 0, *vec_a2, angles='xy', scale_units='xy', scale=1, color="r")

def subtitle(ax, txt):
    plt.annotate(txt, xytext=(.5, -.15), xy=(.5, -.15), textcoords=ax, ha="center", fontsize=16)



A = sp.matrix([[8, -2],
               [3, 1]])
# 1 (1, 1)
ax = fig.add_subplot(gs[0, 0])
ax.axis("equal")
# basic
view_basic_vec(ax)

xy_coord = sp.array([1, 1])
ax.scatter(*xy_coord, color="k")
ax.text(*xy_coord, "$({},{})".format(*xy_coord)+"_{xy}$")

# Ax
Ax = sp.asarray(A.dot(xy_coord)).reshape(-1)
# print("Ax", Ax)
ax.scatter(*Ax)

for x in sp.linspace(0, .1, 5):
    for y in sp.linspace(0, .1, 5):
        v = sp.array([x, y])
        map_v = sp.asarray(A.dot(v)).reshape(-1)
#         print(map_v)
        ax.scatter(x, y, color="k", s=1)
#         plt.scatter(*map_v, color="r")
#         ax.arrow(0, 0, *map_v,**e_vec_prop)

eig_vs, eig_vecs = linalg.eig(A)
# print(eig_vecs)
view_eig_vec(ax)

P = sp.matrix(eig_vecs)
P_inv = linalg.inv(P)
coord_xy2eigvec = P_inv.dot(xy_coord)
# print(coord_xy2eigvec)
proj_a1_coord_xy2eigvec_xy = eig_vecs.dot(sp.array([coord_xy2eigvec[0], 0]))
proj_a2_coord_xy2eigvec_xy = eig_vecs.dot(sp.array([0, coord_xy2eigvec[1]]))
# print("on a1 2xy", proj_a1_coord_xy2eigvec_xy)
# print("on a2 2xy", proj_a2_coord_xy2eigvec_xy)
ax.quiver(0, 0, *proj_a1_coord_xy2eigvec_xy, color="b", angles='xy', scale_units='xy', scale=1)
ax.quiver(0, 0, *proj_a2_coord_xy2eigvec_xy, color="b", angles='xy', scale_units='xy', scale=1)
ax.arrow(*proj_a1_coord_xy2eigvec_xy, *(xy_coord - proj_a1_coord_xy2eigvec_xy), color="b", linestyle="--")
ax.arrow(*proj_a2_coord_xy2eigvec_xy, *(xy_coord - proj_a2_coord_xy2eigvec_xy), color="b", linestyle="--")
ax.text(*coord_xy2eigvec, "$({:.3f},{:.3f})".format(*coord_xy2eigvec) + "_{\\vec{a}_{1}\\vec{a}_{2}}$", color="b", fontsize=12)
ax.grid(True)
subtitle(ax, "$(\\alpha_{1}, \\alpha_{2})^{T}=\mathbf{P}^{-1}(x, y)^{T}$")
ax.set_title("(1) change $\\vec{e}_{1}, \\vec{e}_{2}(xy)$ to eigen basis $\\vec{a}_{1}, \\vec{a}_{2}$")


#2 scale eig
ax = fig.add_subplot(gs[0, 1])
ax.axis("equal")

view_eig_vec(ax)
ax.scatter(*Ax)

ax.quiver(0, 0, *proj_a1_coord_xy2eigvec_xy, color="b", angles='xy', scale_units='xy', scale=1)
ax.quiver(0, 0, *proj_a2_coord_xy2eigvec_xy, color="b", angles='xy', scale_units='xy', scale=1)

scale_proj_a1_coord_xy2eigvec = eig_vs[0].real * coord_xy2eigvec[0]
scale_proj_a1_coord_xy2eigvec_xy = eig_vecs.dot(sp.array([scale_proj_a1_coord_xy2eigvec, 0]))
scale_proj_a2_coord_xy2eigvec = eig_vs[1].real * coord_xy2eigvec[1]
scale_proj_a2_coord_xy2eigvec_xy = eig_vecs.dot(sp.array([0, scale_proj_a2_coord_xy2eigvec]))

# print("2:")
# print("scale_a1",scale_proj_a1_coord_xy2eigvec)
# print("scale_a2",scale_proj_a2_coord_xy2eigvec)
# print("scale_a1_xy",scale_proj_a1_coord_xy2eigvec_xy)
# print("scale_a2_xy",scale_proj_a2_coord_xy2eigvec_xy)
ax.quiver(0, 0, *scale_proj_a1_coord_xy2eigvec_xy, color="g", angles='xy', scale_units='xy', scale=1)
ax.quiver(0, 0, *scale_proj_a2_coord_xy2eigvec_xy, color="g", angles='xy', scale_units='xy', scale=1)
ax.text(*scale_proj_a1_coord_xy2eigvec_xy/2 - .5, "${:.1f}".format(eig_vs[0].real) + "\\vec{a}_{1} \\alpha_{1}$", color="g")
ax.text(*scale_proj_a2_coord_xy2eigvec_xy/2 + .5, "${:.1f}".format(eig_vs[1].real) + "\\vec{a}_{2} \\alpha_{2}$", color="g")

scale_a1a2_xy = sp.asarray(eig_vecs.dot(sp.array([scale_proj_a1_coord_xy2eigvec, scale_proj_a2_coord_xy2eigvec]))).reshape(-1)
ax.text(scale_a1a2_xy[0]-2,scale_a1a2_xy[1]-2.3, "$({:.3f},{:.3f})".format(scale_proj_a1_coord_xy2eigvec,scale_proj_a2_coord_xy2eigvec) + "_{\\vec{a}_{1}\\vec{a}_{2}}$" , color="g", fontsize=12)
ax.arrow(*scale_proj_a1_coord_xy2eigvec_xy, *(scale_a1a2_xy - scale_proj_a1_coord_xy2eigvec_xy), color="g", linestyle="--")
ax.arrow(*scale_proj_a2_coord_xy2eigvec_xy, *(scale_a1a2_xy - scale_proj_a2_coord_xy2eigvec_xy), color="g", linestyle="--")

ax.grid(True)
ax.set_title("(2) scale to eigenV $\lambda_{1} \\alpha_{1}, \lambda_{2}\\alpha_{2}$ along eigen vec axis")
subtitle(ax, "$(\\alpha_{1}', \\alpha_{2}')^{T}=(\\lambda_{1} \\alpha_{1}, \\lambda_{2}\\alpha_{2})^{T}$")
    
# 3. P 
ax = fig.add_subplot(gs[0, 2])
ax.axis("equal")
# basic
view_basic_vec(ax)

vec_xy = sp.asarray(P.dot(sp.array([scale_proj_a1_coord_xy2eigvec, scale_proj_a2_coord_xy2eigvec]))).reshape(-1)
# print(vec_xy)
ax.scatter(*Ax)
ax.quiver(0, 0, *vec_xy, color="k", angles='xy', scale_units='xy', scale=1)
ax.text(vec_xy[0]-1.2,vec_xy[1]+.2, "$({:.1f},{:.1f})".format(*vec_xy) + "_{xy}$" , color="k", fontsize=12)
plt.vlines(vec_xy[0], 0, vec_xy[1], linestyle="--")
plt.hlines(vec_xy[1], 0, vec_xy[0], linestyle="--")
ax.grid(True)
ax.set_title("(3) reverse to basis $\\vec{e}_{1},\\vec{e}_{2}(xy)$")
subtitle(ax, "$(x', y')^{T}=\mathbf{P}(\\alpha_{1}', \\alpha_{2}')^{T}$")


plt.tight_layout()
plt.savefig("s.png")
plt.show()

応用例

行列の累乗が容易に計算できるようになる．

$\begin{eqnarray} \mathbf{A}^{n} &=& (\mathbf{P}\mathbf{\Lambda}\mathbf{P}^{-1}) (\mathbf{P}\mathbf{\Lambda}\mathbf{P}^{-1}) \cdots (\mathbf{P}\mathbf{\Lambda}\mathbf{P}^{-1})\\\\ &=& \mathbf{P}\mathbf{\Lambda}(\mathbf{P}^{-1}\mathbf{P}) \mathbf{\Lambda}\mathbf{P}^{-1} \cdots (\mathbf{P}\mathbf{\Lambda}\mathbf{P}^{-1})\\\\ &=& (\mathbf{P}\mathbf{\Lambda}^{n}\mathbf{P}^{-1}) \end{eqnarray}$

マルコフ遷移の遷移行列Tの固有値$|\lambda| \leq 1$をみれば発散するか収束（極限）があるかがわかる．
制御では，微分方程式のシステム行列の固有値からシステムの発散，安定がわかる．常微分方程式の解から実時間空間に戻すと$Ce^{\mathbf{A}t}$になるので，固有値の実部がマイナスであれば時間遷移していっても安定するとわかる．