統計検定準一級の２２章は主成分分析ですが、内容はかなり詳細で理解に苦しみました。そこでまとめてみました。sympyのプログラミング、論理の整合性のチェックなどにChatGPT４(open AI)を使っています。いつも通り間違いがあると思いますので、ご指摘ただけると助かります。

関連するsympyのシートは
https://github.com/innovation1005/qiita_innovation1005

また、日曜の朝9:30と火曜の午後7:00から勉強会を開いています。
https://study-data-analysis.connpass.com/
参考文献は統計学実践ワークブックです。今回の投稿は２２章に関するものです。

要約

平均
$n$個の個体があり、それは$p$個の変数で構成されているとします。

これらを $n \times p$行列としてまとめたものを$X$とします。

$X=(x_{i,j})_{1\le i \le n,1\le j \le p}$

X=
\begin{pmatrix}
x_{11} & \cdots & x_{1j} & \cdots & x_{1p}\\
\vdots & \ddots &        &        & \vdots \\
x_{i1} &        & x_{ij} &        & x_{ip} \\
\vdots &        &        & \ddots & \vdots \\
x_{n1} & \cdots & x_{nj} & \cdots & x_{np}
\end{pmatrix}

行列$X$の行と列の和を明確にするために、ドット記号を使って

$x_{\cdot ,j}=\sum_{i=1}^n x_{i,j} (j=1,\cdots,p)$

$x_{i,\cdot}=\sum_{j=1}^p x_{i,j} (i=1,\cdots,n)$

$x_{\cdot,\cdot}=\sum_{i=1}\sum_{j=1}^p x_{i,j} (i,j=1,\cdots,n)$

と書きます。変数$j$に対する$n$個の平均$\bar{x}_{\cdot,j}$ は

\bar{x}_{\cdot,j} =x_{\cdot,j} /n

となります。

内積

ベクトル$\textbf{x}=(x_1,\cdots,x_p)$,$\textbf{y}=(y_1,\cdots,y_p) \in \mathbb{R}^p$に

内積を

$$\langle\textbf{x},\textbf{y}\rangle=\textbf{x}^T\textbf{y}=\sum_{i=1}^p x_iy_i$$

で定義します。

直交行列と固有値

直交行列とは$U^TU=UU^T=I_p$となる$p$次正方行列です。$p$次行列$A$を実対称行列とするとき、適当な直交行列$U=(u_1,\cdots,u_p)$をとれば

$U^T A U = diag(\lambda_1,\cdots,\lambda_p)$

この式に左から$U$を掛けると$AU=Udiag(\lambda_1,\cdots,\lambda_p)$となり、固有値問題

A\textbf{u}_j=\lambda_j\textbf{u}_j,\langle\textbf{u}_j\textbf{u}_k\rangle=\delta_{j,k}

$\delta_{j,k}$はクロネッカーのデルタを表します。

標本分散共分散行列と固有値

観測データ$B$を$n \times p$行列とし、$A=B^TB$となっている場合は、$A$は分散共分散行列$S$や標本相関行列$R$であることを前提としています。

主成分

$j$番目の固有値$\lambda_j$に対応する固有ベクトル$\textbf{u}_j$と元の変数$\textbf{x}=(x_1,\cdots,x_p)^T$との内積を変数$y_j$で表し、第$j$主成分とします。

y_j=\langle\textbf{u}_j,\textbf{x}\rangle=x_iu_{1,j}+\cdots+x_pu_{p,j}

$y_j$ は $j$ 列の $n$ 個について、データ点(個体)を第 $j$ 主成分へ射影した射影値(主成分得点)です。

主成分得点

つぎに$\textbf{x}_i$を定義します。これは元データと中心化されたデータのベクトルの両方を意味します。

\textbf{x}_i=(x_{i,1}, x_{i,2}, \ldots, x_{i,p})^T

または中心化された

\textbf{x}_i=(x_{i,1}-\bar{x}_{\cdot,1},\cdots,x_{i,p}-\bar{x}_{\cdot,p})^T

です。

このベクトルを主成分に代入した
$${y_{i,j}=\langle\textbf{x}_i,\textbf{u}_j\rangle|i=1,\cdots,n,j=1,\cdots,p}$$
を主成分得点といいます。

中心化された$\textbf{x}_i,i=1,\cdots,n$を利用して標本の大きさ$n$の第$j$主成分$y_j$を

\textbf{y}_i=(y_{1,j},\cdots,y_{n,j})^T=(\langle\textbf{x}_1,\textbf{u}_j\rangle,\cdots,\langle\textbf{x}_n,\textbf{u}_j\rangle)^T

とおきます。

yjの分散と固有値

$\textrm{V}[y_j]=\frac{1}{n-1}||\textbf{y}_j||^2=\textbf{u}_j^TS\textbf{u}_j=\lambda_j$

yjとykの共分散と固有値

$\textrm{Cov}[y_j,y_k]=\sqrt{\lambda_j\lambda_k}\delta_{j,k}$

ここで、一般的に中心化されたSから求められた２つのi.j主成分は直交していますので$i\ne j$ではゼロです。したがって、

$\textrm{Cov}[y_j,y_k]=\lambda_j\delta_{j,k}$

と書くのが普通です。主成分がどのように求められたのかについて知っている必要があります。

主成分の持つ意味

\lambda_j\textbf{u}_j=S\textbf{u}_j=\frac{1}{n-1}X_C^T(X_C\textbf{u}_j)=\frac{1}{n-1}X_C^T\textbf{y}_j=
\begin{pmatrix}
\textrm{Cov}[x_1,y_j]\\
\vdots\\
\textrm{Cov}[x_p,y_j]
\end{pmatrix}

$k$行目について
$\mathrm{Cov}[x_k,y_j]=\lambda_ju_{k,j}$
となります。

主成分負荷量

元の変数と主成分との相関係数を主成分負荷量という。

$r_{y_j,x_k}=\frac{\textrm{Cov}[x_k,y_j]}{\sqrt{V[y_j]s_{k,k}}}=\frac{\sqrt{\lambda_j}u_{k,j}}{\sqrt{s_{k,k}}}$

$s_{k,k}$は$k$次の分散を表します。

sympyによる確認

観測値をx,yで表現しています。

import sympy as sp
sp.init_printing(use_latex='mathjax')

sp.var('x11, x12, x13, x21, x22, x23, x31, x32, x33') 
sp.var('y11, y12, y13, y21, y22, y23, y31, y32, y33') 
X = sp.Matrix([
[x11,x12,x13],
[x21,x22,x23],
[x31,x32,x33]
])
Y = sp.Matrix([
[y11,y12,y13],
[y21,y22,y23],
[y31,y32,y33]
])
display(X)
display(Y)

内積

ベクトル$\textbf{x}=(x_1,\cdots,x_p)$,$\textbf{y}=(y_1,\cdots,y_p) \in \mathbb{R}^p$に
内積$\langle \textbf{x},\textbf{y}\rangle$を
$\langle \textbf{x},\textbf{y} \rangle =\textbf{x}^T\textbf{y}=\sum_{i=1}^p x_iy_i$で定義する。

bf_x=X.col(0)
bf_y=Y.col(0)
display(bf_x)
display(bf_y)
display(sp.transpose(bf_x)*bf_y)

bf_x=X.row(0)
bf_y=Y.row(0)
display(bf_x)
display(bf_y)
display(sp.transpose(bf_x)*bf_y)

通常の Euclidean 内積の場合、$\langle \textbf{a}, \textbf{b}\rangle$ では $\textbf{a}$ と $\textbf{b}$ は列ベクトル（縦ベクトル）として扱われます。一方で、内積の定義が異なる場合や特定の文脈では、異なる方向のベクトルを取ることもありますが、その際には適切なベクトルの形を代入するように注意する必要があります。

直交行列

直交行列とは$U^TU=UU^T=I_p$となる$p$次正方行列です。$p$次行列$A$を実対称行列とするとき、適当な直交行列$U=(u_1,\cdots,u_p)$をとれば

$U^T A U = diag(\lambda_1,\cdots,\lambda_p)$

この式に左から$U$を掛けると$AU=Udiag(\lambda_1,\cdots,\lambda_p)$となり、固有値問題

A \textbf{u}_j=\lambda_j \textbf{u}_j, \langle\textbf{u}_j\textbf{u}_k\rangle= \delta_{j,k}

直交行列を得る方法にはいろいろあります。直交の条件を与えて最適化する。回転行列、置換行列、反射行列、QR分解など

# 変数の定義
u11, u12, u21, u22 = sp.symbols('u11 u12 u21 u22')

# 直交行列の要素に関する関係
relation1 = sp.Eq(u11**2 + u21**2, 1)
relation2 = sp.Eq(u12**2 + u22**2, 1)
relation3 = sp.Eq(u11*u12 + u21*u22, 0)

# u11, u12, u21, u22に関する関係を解く
solutions = sp.solve([relation1, relation2, relation3], (u11, u12, u21, u22))

# 解を表示
print("Solution:")
for solution in solutions:
    display(solution)

# 解を代入して直交行列を作成
U = sp.Matrix([
    [solutions[0][0], solutions[0][1]],
    [solutions[0][2], solutions[0][3]]
])

# 直交行列を表示
print("Orthogonal Matrix U:")
display(U)
display(sp.simplify(U.T*U))

この方法では実際に二次元までしか解けません。もっとも簡単に直交行列を意識できるのは三角関数を利用したものです。

#三角関数を利用した直交行列の構築

# 変数の定義
theta = sp.symbols('theta')

# 直交行列の定義
U = sp.Matrix([
    [sp.cos(theta), -sp.sin(theta)],
    [sp.sin(theta), sp.cos(theta)]
])

# Uを表示
display(U,sp.simplify(U.T*U))

# 実対称行列の固有値分解
sp.var('a11, a12, a21, a22')
# 行列の定義
A = sp.Matrix([
    [a11, a12],
    [a12, a22]
])

U, D = A.diagonalize()
display(U)
display(D)
display(sp.simplify(U * D * U.inv()))

Aが標本分散共分散行列Sや相関行列R

Bを$n \times p$行列として$A=B^TB$とすると、$A$の固有値は非負となります。一般に$B$は中心されていいます。

A\textbf{u}_j=\lambda_j \textbf{u}_j, \langle \textbf{u}_j,\textbf{u}_k \rangle=\delta_{j,k}

の両辺に$\textbf{u}_j^T$を掛けると

\textbf{u}_j^TA\textbf{u}_j=\textbf{u}_j^TB^TB\textbf{u}_j=(B\textbf{u}_j)^TB\textbf{u}_j=\textbf{u}_j^T\lambda_j \textbf{u}_j=\lambda_j

$\lambda_j =||Bu_j||^2\ge0$

sp.var('b11, b12, b21, b22, b31, b32') 
B = sp.Matrix([
[b11,b12],
[b21,b22],
[b31,b32]
])
display("B",B)
A=B.T*B
display("A",A)
U, D = A.diagonalize()
display("A from U,D",sp.simplify(U * D * U.inv()))

主成分

ｊ番目の固有値$\lambda_j$に対応する固有ベクトル$\textbf{u}_j$と元の変数 $\textbf{x}=(x_1,\cdots,x_p)^T$との内積を変数$y_j$で表し、第$j$主成分という。

y_j=\langle \textbf{u}_j,\textbf{x} \rangle= x_iu_{1,j}+\cdots+x_p u_{p,j}

X = sp.Matrix([
[x11,x12,x13],
[x21,x22,x23],
[x31,x32,x33]
])

U = sp.Matrix([
[u11,u12,u13],
[u21,u22,u23],
[u31,u32,u33]
])

bf_x=X.col(0)#\textbf{x}=(x_1,\cdots,x_p)^Tが縦ベクトル。(x_1,\cdots,x_p)は横ベクトル
j=0
bf_u_j=U.col(j)
display(bf_x)     
display(bf_u_j)
y_j=sp.transpose(bf_u_j)*bf_x
display(y_j)

主成分得点

\textbf{x}_i=(x_{i,1},x_{i,2},\ldots,x_{i,p})^T

、または中心化された

\textbf{x}_i=(x_{i,1}-\bar{x}_{\cdot,1},\cdots,x_{i,p}-\bar{x}_{\cdot,p})^T

について主成分に代入した

\{y_{i,j}=\langle\textbf{x}_i,\textbf{u}_j\rangle|i=1,\cdots,n,j=1,\cdots,p\}

を主成分得点という。$y_{i,j}$は個体１つの第j主成分分析ですが、${\cdot}$で囲われているので、$n \times p$個の主成分得点を表します。

i=0
bf_x_i=sp.transpose(X.row(i))
j=0
bf_u_j=U.col(j)
y_ij=sp.transpose(bf_x_i)*bf_u_j
display(bf_x_i)
display(bf_u_j)
display(y_ij)

yjの分散と固有値

$\textrm{V}[y_j]=\frac{1}{n-1}||\textbf{y}_j||^2=\textbf{u}_j^TS\textbf{u}_j=\lambda_j$

# 元データ行列Xの定義
X = sp.Matrix([
    [1, 2],
    [3, 4],
    [4, 6]
])

# Xの列数（特徴数）を取得
p = X.shape[1]
n = X.shape[0]

# Xの各列の平均を計算

# 偏差の計算
Xc=(sp.eye(n)-sp.ones(n,n)/n)*X
#display(sp.simplify(Xc))

# 共分散行列の計算
S = (Xc.T @ Xc) / (n-1 )
# 固有値と固有ベクトルの計算
P,D=S.diagonalize()
adjs=[(sp.transpose(P.col(i))*P.col(i)) for i in range(p)]
U=[sp.simplify(P.col(i)) for i in range(p)]

# y_jの計算と分散の計算
variance=[]
yjs=[]
for u in U:
    display(u)
    v=0
    yj=[]
    for i in range(n):
        y=Xc.row(i) * u
        v+=y[0]**2/(n-1)
        yj.append(y)
    variance.append(v)
    yjs.append(yj)
adj1=adjs[0]
adj2=adjs[1]
display(sp.simplify(variance[0]/adj1[0]))
display(sp.simplify(variance[1]/adj2[0]))
display(D)

sympyのeigenvectorは$a^2+b^2=1$になっていないので注意。

yjとykの共分散と固有値

$\textrm{Cov}[y_j,y_k]=\sqrt{\lambda_j\lambda_k}\delta_{j,k}$

display(S)
yj=sp.Matrix(yjs[0])
yk=sp.Matrix(yjs[1])
sp.simplify(yj.T*yk)

主成分の持つ意味

\lambda_j\textbf{u}_j=
\begin{pmatrix}
\lambda_j\textbf{u}_{1,j}\\
\vdots\\
\lambda_j\textbf{u}_{p,j} 
\end{pmatrix}=
\begin{pmatrix}
S\textbf{u}_{1,j}\\
\vdots\\
S\textbf{u}_{p,j} 
\end{pmatrix}=
\frac{1}{n-1}\begin{pmatrix}
\textbf{x}_1^T\textbf{x}_1\textbf{u}_{1,j}\\
\vdots\\
\textbf{x}_p^T\textbf{x}_p\textbf{u}_{p,j} 
\end{pmatrix}

=
\frac{1}{n-1}\begin{pmatrix}
\textbf{x}_1^T\textbf{y}_j\\
\vdots\\
\textbf{x}_p^T\textbf{y}_j 
\end{pmatrix}=
\begin{pmatrix}
\textrm{Cov}[x_1,y_j]\\
\vdots\\
\textrm{Cov}[x_p,y_j]
\end{pmatrix}

$k$行目とは$\mathrm{Cov}$に代入した$X_C^T$の$k$行目のことなので、$X_C$の$k$列目です。
$\mathrm{Cov}[x_k,y_j]=\lambda_ju_{k,j}$
となります。

L=D*sp.Matrix([1,1])
for u,la in zip(U,L):
    display(sp.simplify(u*la))

for u in U:
    display(S*u)

#共分散による主成分の理解　右辺
for j in range(p):
    x=Xc.col(j)
    for k in range(p):
        print("j=",j,"k=",k)
        y=yjs[j][:]
        cov=0
        for i in range(n):
            cov+=y[i][0]*x[i]
        display(cov)

主成分負荷量

元の変数と主成分との相関係数を主成分負荷量という。

$r_{y_j,x_k}=\frac{\textrm{Cov}[x_k,y_j]}{\sqrt{V[y_j]s_{k,k}}}=\frac{\sqrt{\lambda_j}u_{k,j}}{\sqrt{s_{k,k}}}$

$x_k$とは$X_C$の$k$列目です。
$\mathrm{Cov}[x_k,y_j]=\lambda_j u_{k,j}$
となります。

$y_j$の分散と固有値の関係は$\textrm{V}[y_j]=\lambda_j$です。よって、

$\frac{\textrm{Cov}[x_k,y_j]}{\sqrt{V[y_j]s_{k,k}}}=\frac{\lambda_j u_{k,j}}{\sqrt{\lambda_js_{k,k}}}=\frac{\sqrt{\lambda_j}u_{k,j}}{\sqrt{s_{k,k}}}$

分散共分散行列ではなく、相関行列をつかって固有値問題を扱った場合、

\lambda_ju_j=
\begin{pmatrix}
\lambda_ju_{1,j}\\
\vdots\\
\lambda_ju_{p,j} 
\end{pmatrix}=
\begin{pmatrix}
\textrm{R}u_{1,j}\\
\vdots\\
\textrm{R}u_{p,j} 
\end{pmatrix}=
\frac{1}{n-1}\begin{pmatrix}
\textbf{x}_{Z,1}^T\textbf{x}_{Z,1}u_{1,j}\\
\vdots\\
\textbf{x}_{Z,p}^T\textbf{x}_{Z,p}u_{p,j} 
\end{pmatrix}

=
\frac{1}{n-1}\begin{pmatrix}
\textbf{x}_{Z,1}^T\textbf{y}_j\\
\vdots\\
\textbf{x}_{Z,p}^T\textbf{y}_j 
\end{pmatrix}=
\begin{pmatrix}
\textrm{Cov_z}[x_1,y_j]\\
\vdots\\
\textrm{Cov_z}[x_p,y_j]
\end{pmatrix}

ここで、$\textbf{x}_{Z,p}$は標準化した$x$を表し、$\textrm{Cov_Z}$は標準化した変数を代入した分散共分散行列を表します。したがって、$V[y_j]=1$とは限りません。ゆえに

$r_{y_j,x_k}=\frac{\textrm{Cov_Z}[x_k,y_j]}{\sqrt{V[y_j]s_{k,k}}}=\frac{\textrm{Cov_Z}[x_k,y_j]}{\sqrt{V[y_j]}}=\frac{\lambda_j}{\sqrt{V[y_j]}}=\sqrt{\lambda_j}\textbf{u}_{k,j}$

多くの文献で説明されないけれども最も重要な点

主成分分析は元の多変量データを線形結合した「主成分」と呼ばれる新しい変数を生成します。つまり、標本データに基づいています。

線形性の仮定: 主成分分析は線形性の仮定が基礎にあります。これは、非線形な関係がある場合には限定的な解釈しかできないという点で注意が必要です。

解釈可能性: 主成分分析の結果は、元の変数に対する新しい「主成分」の線形結合として表されます。これは、新しい主成分が元のデータの何を表しているのか解釈する際に注意が必要です。

適用範囲: 主成分分析は探索的な分析手法ともされ、データのパターンを理解するための一手段です。そのため、標本に基づいているという制限は、その目的と照らし合わせて考慮されるべきです。

標本サイズと品質: 実際に主成分分析を適用する場合、標本サイズやデータの品質（例：外れ値、欠損値等）も考慮する必要があります。

2023/9/5 1947

疑問点1：共分散行列を使うべきか相関行列を使うべきか？

100点満点の複数の試験のようにデータの単位がそろっている場合には、分散共分散行列と相関行列のどちらを使って、主成分分析をしてもよい。しかし、単位がそろっていない場合は分析結果が単位に依存するため、相関行列を使って無単位化を行った主成分分析が望ましい。というのは本当でしょうか？

主成分分析の目的は、多次元データセットの主要な「方向」を見つけ出すことです。その「方向」は元のデータのスケールや単位に依存します。

単位が揃っている場合：たとえば、全ての試験が100点満点の場合、各試験のスケール（単位）は既に同じです。そのため、分散共分散行列を用いた主成分分析は問題ありません。

単位が揃っていない場合：異なる単位を持つ変数を分析する場合、スケールの大きい変数が結果に大きな影響を与える可能性があります。このような場合、相関行列に基づいて主成分分析を行い、データを無単位化することが一般的です。

分散共分散行列に基づく主成分分析は、スケールの大きな変数に結果が偏る可能性があります。相関行列を使用すると、この問題を緩和できます。しかし、スケールの違いがデータセットの特徴を表しているかもしれません。したがって、分析の前にどのようなデータを扱っているのか、その特性を理解することが重要です。それに基づいて、分散共分散行列を使うか、相関行列を使うかを決定すると良いでしょう。

主成分分析への入力は分散共分散行列や相関行列だけではありません。つぎの条件を満たす行列であれば問題ありません。

PCAの入力として分散共分散行列か相関行列のどちらを用いるべきかについては、具体的な分析の目的やデータの性質に依存します。こちらに一概に適用できる「正解」は存在しないため、注意が必要です。

行列の性質について

正方行列: これは固有値と固有ベクトルを求めるための基本的な条件です。
対称行列: 対称行列であれば、その固有ベクトルも直交するという便利な性質があります。
正定値または半正定値: 正定値または半正定値でないと、固有値が負になる可能性があります。

各行列のこれらの性質は、数学的な計算が正確かつ意味のあるものになるために重要です。どの行列を選ぶかは最終的には解析の目的に依存します。特定の目的に最も適した行列を選ぶことが重要です。

対称行列

対称行列は数学と工学の多くの分野で重要な役割を果たします。その基本的な性質はつぎのとおりです：

対称性: 対称行列はその転置と等しい（$A = A^T$）。
正方形: 対称行列は必ず正方行列です。
固有値と固有ベクトル: 対称行列の固有値は全て実数です。また、異なる固有値に対応する固有ベクトルは直交します。
対角化可能: 対称行列は常に対角化可能です。すなわち、$A = PDP^{-1}$ と分解できるとき、$P$ はその固有ベクトルで形成され、$D$ は対角行列でその対角成分が固有値です。さらに、$P$ は直交行列であることが多い（$P^T = P^{-1}$）。
分解: 対称行列は特定の方法で分解できます。例えば、コレスキー分解や固有値分解があります。
エネルギーの保存: 物理学の文脈では、対称行列はしばしばエネルギーの保存則と関連しています。
計算効率: 対称行列の性質により、行列計算が簡単になる場合があります。
正定値、半正定値: 対称行列はしばしば正定値または半正定値です。
トレースと行列式: 対称行列のトレース（対角成分の和）は固有値の和と等しく、行列式は固有値の積と等しい。

画像解析

共分散行列と相関行列の選択は重要です。かつ、その重要性は具体的な解析目的や画像データの性質に依存します。

共分散行列: 通常、各ピクセルの明るさ（または色）が特徴量となります。すべてのピクセルが同じ単位（例えば、0-255の整数値）で表されている場合、共分散行列を用いることが多いです。
相関行列: 一方で、もし異なる種類の画像特徴（例えば、色とテクスチャ）を組み合わせる必要がある場合や、スケールの違いが解析に影響を与える可能性がある場合は、相関行列を用いてスケーリングの影響を排除することが推奨されます。

エネルギー解析

共分散行列と相関行列の選択は非常に重要です。

共分散行列: エネルギー解析で通常使用される電力消費、温度、湿度などではすべて同じ単位で表されています。または単位の違いが解析結果に重要でないと考えられる場合、共分散行列が使用されます。
相関行列: 一方、解析に用いる各種指標が異なる単位で表されている、またはスケールの違いが解析結果に影響を与える可能性がある場合は、相関行列を使用することが推奨されます。
解析目的: 共分散行列と相関行列の選択は、解析の目的にも強く依存します。たとえば、ある特定のエネルギー変数（例：電力消費）に最も影響を与える要因を特定することが目的であれば、その変数と他の変数との関係性に注目し、適切な行列を選択する必要があります。

共分散行列を主に使う分野

データの単位が揃っている、または単位の違いが分析において重要でないとされるケースが多いです。たとえば、

金融のポートフォリオ最適化: 金融業界では、異なる資産間のリスクとリターンの関係を評価するために共分散行列がしばしば使用されます。
信号処理: 同じ種類の信号を扱う場合、共分散行列が用いられることが多いです。
品質管理: 製造プロセスにおいて、同じ単位で計測される複数の品質指標を同時に監視する必要がある場合、共分散行列が用いられます。
気象学: 同じ単位で表される気象データ(気温、湿度など)の関連性を調べる場合にも共分散行列が使用されます。
機械学習: 特に教師なし学習において、データの単位が同じであれば、共分散行列を用いてデータの構造を理解することがあります。

疑問点2：複数の主成分ベクトルの表現

固有値分解から得られる特定の固有値に対する固有ベクトルは無数にあります。しかし、主成分分析の場合、元データの分散を維持するために、固有ベクトルの要素の２乗の和が１になるという制限が設けられています。したがって、ランク$r$の固有ベクトルは$2r$の$\pm$の組み合わせに制限されます。もともと固有ベクトルは方向性を示す尺度で、$r$個ある固有値に対してそれぞれについて2種類存在する固有ベクトルの表現は同じ方向を指しています。

たとえば、２次元の場合、固有ベクトルは原点を通る直線になります。これに要素の２乗和が１という制限を加えると解は2つに限定されます。２つの解の違いは±の符号です。つまりどちらも同じ方向を向いています。この際に、違う符号の解釈と表現には注意が必要です。誤解を招かないようにしなければなりません。

例１

import numpy as np
np.set_printoptions(precision=2)
import pandas as pd
pd.set_option('display.float_format', lambda x: '%.4f' % x)
import matplotlib.pyplot as plt
import seaborn as sns

# あるクラスの6人（No.1〜6）の生徒の国語(x1)、数学(x2)、理科(x3)、社会(x4)
# の小テストの結果

X = {'x1':[2,9,8,7,2,5],
     'x2':[2,8,3,1,9,4],
     'x3':[3,10,2,3,8,5],
     'x4':[1,9,7,8,2,5]}
index = ['No.1','No.2','No.3','No.4','No.5','No.6']

X = pd.DataFrame(X, index=index)
pd.set_option('display.float_format', lambda x: '%.2f' % x)
print(X)
XC=(np.identity(len(X))-1/len(X)*np.ones((len(X),len(X))))@X
print(XC)
# 標本分散共分散行列
S = 1/(len(XC)-1)*XC.T@XC
S

R = np.corrcoef(X.T)
print('標本相関行列')
display(pd.DataFrame(R,index=X.columns, columns=X.columns))
# 標本分散共分散行列から固有値・固有ベクトルを求める
eigvals, eigvecs = np.linalg.eigh(S)

eigvals = eigvals
eigvecs = eigvecs

print(f'固有値:\n{eigvals[::-1]}') # 大きい順にソート

print(f'固有ベクトル:\n{eigvecs}') # 第二主成分以降の符号が逆になっている

df_eigvecs = pd.DataFrame(np.c_[eigvecs[:, -1], eigvecs[:, -2], eigvecs[:, -3]\
                                , eigvecs[:, -4]], index=X.columns\
                          ,columns=['PC1','PC2','PC3','PC4']).T
print(df_eigvecs)

結果の解釈：

第一主成分：x1,x2,x3,x4はすべて正なので総合評価を表しています。
第二主成分：x1,x4が同じ方向を示し、x2,x3がそれとは逆の方向を示しています。
第三主成分：x1,x2が同じ方向を示し、x3,x4はそれとは逆の方向を示しています。
第四主成分：x1,x3が同じ方向を示し、x2,x4はそれとは逆の方向を示しています。

特定の固有値に対する固有ベクトルは無数に存在する可能性があります。そして、その方向は一意です。さらに、主成分分析の文脈では通常、固有ベクトルはベクトルの要素の二乗和が1になるように正規化されたベクトルとして扱われます。

$n$次元であっても固有ベクトルは直線です。一般に、固有ベクトルはある線形変換が作用するときにその方向が変わらないベクトルです。つまり、この線形変換を適用しても、ベクトルは同じ「方向」にあるというわけです。

固有ベクトルはその「方向」を表すものであり、スカラー倍（この場合は-1倍を含む）したベクトルもまた同じ固有値に対する固有ベクトルです。つまり、ある固有値に対する固有ベクトルが$ \vec{v} $であれば、$-\vec{v} $もまたその固有値に対する固有ベクトルです。

具体的には、ある行列$ A $と固有値$ \lambda $に対して、次の式が成り立つ場合、

$$ A\vec{v} = \lambda \vec{v} $$

この式の両辺に-1をかけても、式は成立します。

$$ -A\vec{v} = -\lambda \vec{v} $$
$$ A(-\vec{v}) = \lambda (-\vec{v}) $$

これが意味するのは、$-\vec{v} $もまた$ \lambda $という固有値に対する固有ベクトルである、ということです。

# 主成分得点
pca_score = XC @ df_eigvecs.iloc[:2, :].T
pca_score.index=index
pca_score

df_eigvecs.loc["PC2"][0]

plt.figure(figsize=(3, 2))
plt.scatter(pca_score.iloc[:,0], pca_score.iloc[:,1])
labels = index
for i, label in enumerate(labels):
    plt.text(pca_score.iloc[:,0][i], pca_score.iloc[:,1][i],label)

arrow_adjust=8
text_adjust=1.1
feature_names = df_eigvecs.index
for i in range(len(feature_names)):
        plt.arrow(0, 0,
                  df_eigvecs.loc['PC1'][i]*arrow_adjust,
                  df_eigvecs.loc['PC2'][i]*arrow_adjust,
                  color='r', head_width=0.15, head_length=0.2)
        plt.text(df_eigvecs.loc['PC1'][i]*arrow_adjust*text_adjust,
                 df_eigvecs.loc['PC2'][i]*arrow_adjust*text_adjust,
                 X.columns[i],
                 color='r')
plt.axhline(0, ls="--", color="y")
plt.axvline(0, ls="--", color="y")
plt.xlabel('PC1')
plt.ylabel('PC2')

これは主成分分析の結果を図で表しています。赤字はそれぞれの要素が＋ーのどちらの側にいるのかを示しています。横軸を見ると全ての要素が＋側にいます。一方で縦軸を見るとx1とx4がマイナス側にいて、x2とx3がプラスの側にいるのが分かります。

主成分ベクトルの正負を逆にしても同じです。

df_eigvecs_pmrev=df_eigvecs*np.transpose([[1,-1,-1,-1],
                             [1,-1,-1,-1],
                             [1,-1,-1,-1],
                             [1,-1,-1,-1]])
pca_score_pmrev = XC @ df_eigvecs_pmrev.iloc[:2, :].T
pca_score_pmrev.index=index
pca_score_pmrev,df_eigvecs_pmrev

plt.figure(figsize=(3, 2))
plt.scatter(pca_score_pmrev.iloc[:,0], pca_score_pmrev.iloc[:,1])
labels = index
for i, label in enumerate(labels):
    plt.text(pca_score_pmrev.iloc[:,0][i], pca_score_pmrev.iloc[:,1][i],label)

arrow_adjust=8
text_adjust=1.1
feature_names = df_eigvecs.index
for i in range(len(feature_names)):
        plt.arrow(0, 0,
                  df_eigvecs_pmrev.loc['PC1'][i]*arrow_adjust,
                  df_eigvecs_pmrev.loc['PC2'][i]*arrow_adjust,
                  color='r', head_width=0.15, head_length=0.2)
        plt.text(df_eigvecs_pmrev.loc['PC1'][i]*arrow_adjust*text_adjust,
                 df_eigvecs_pmrev.loc['PC2'][i]*arrow_adjust*text_adjust,
                 X.columns[i],
                 color='r')
plt.axhline(0, ls="--", color="y")
plt.axvline(0, ls="--", color="y")
plt.xlabel('PC1')
plt.ylabel('PC2')

ベクトルの方向については特徴は同じです。

# 主成分負荷量
pc_loading = np.sqrt(eigvals[::-1]).reshape(-1,1) * df_eigvecs / np.sqrt(np.diag(S))
pc_loading.index = df_eigvecs.index
print(f'主成分負荷量：')
display(pc_loading)

ベクトルの方向については特徴は同じです。

ーーーーーーーーーーーーーーーーーーーーーーーー
Python3ではじめるシステムトレード【第2版】環境構築と売買戦略

「画像をクリックしていただくとpanrollingのホームページから書籍を購入していただけます。

Python3ではじめるシステムトレード改訂版　アマゾンから購入いただけます。