Pythonで1から行列の固有値全て求めてみた~べき乗法~
0. はじめに
Pythonで固有値を求めたい.
人生でそう思う時はよくありますよね.
学校で主成分分析をすることになったのですが, そのためにPythonで固有値を求める必要がありました.
例えば$${\bm{A} = \left( \begin{array}{cc} 1 & 2 \\ 2 & 4 \end{array} \right)}$$
の行列式を求めたい.
そんな時使えるコードがこちらです.
import numpy as np
import numpy.linalg as LA
A = np.array([[1,2],[2,4])
LA.eig(A)これで求まりました！
いやこんなの認めない.
こんなパッケージに頼るわけにいかない.
そんな人のための記事です.
ここから, 以下の順序で話していきます.
また今回考える行列は正定値対称行列だとします.
(動機が主成分分析だったためです. 拡張は頑張るとできると思います.)
冪乗法とは？
n個の固有値
実際のコード
より良いアルゴリズム
最初に断っておきますが, 実用性は皆無です.
固有値求めたいなら, 上で書いたようなパッケージを使いましょう.
1. べき乗法とは
べき乗法とはその名の通り, 行列のべき乗(2乗,3乗,…,)を考えていくという手法です.
今, 固有値を求めたい正定値対称行列$${\bm{A} \in \mathbb{R}^{n\times n}}$$
とします.
この時, 固有値を$${\lambda_1,…,\lambda_n}$$
($${\lambda_1 > \cdots > \lambda_n}$$ とします),
固有ベクトルを$${\bm{v}_1,…,\bm{v}_n}$$(これらは全て直交する単位行列とします)とすると,
$${A = \lambda_1 \bm{v}_1 \bm{v}_1^\top + \cdots + \lambda_n \bm{v}_n \bm{v}_n^\top}$$
と書けます.(固有値分解については以下を参照してください)https://qiita.com/KodaiS/items/b72ed968c480ee0f6aef
ここで, あるベクトル$${\bm{x} \in \mathbb{R}^n}$$が$${c_1,…,c_n\in \mathbb{R}}$$ として,
$$
\bm{x} = c_1\bm{v}_1 + \cdots c_n\bm{v}_n
$$
で表されるとします.
この時,
$$
\bm{Ax} = \lambda_1\ c_1 \bm{x}_1 + \cdots + \lambda_n c_n \bm{x}_n
$$
となります.
よって,
$$
\bm{A}^k\bm{x} = \lambda_1^k\ c_1 \bm{x}_1 + \cdots + \lambda_n^k c_n \bm{x}_n
$$
となります.
ここで,$${\bm{y}^k = \frac{\bm{A}^k\bm{x}}{\|\bm{A}^k\bm{x}\|} }$$とします.
すると, $${c_1\neq 0}$$ならば,
$$
\bm{y}^k =\lambda_1^k c_1' \left\{ \bm{v}_1 +\left( \frac{\lambda_2}{\lambda_1} \right)^k \frac{c_2'}{c_1'}\bm{v_2} \cdots + \left( \frac{\lambda_n}{\lambda_1} \right)^k \frac{c_2'}{c_1'}\bm{v_n} \right\}
$$
が成り立ちます.
$${\left\| \frac{\lambda_i}{\lambda_1} \right\| < 1\ (i>1)}$$なので,
$$
\lim_{k\rightarrow \infty} \left( \frac{\lambda_n}{\lambda_1} \right)^k = 0
$$
が成り立ちます. よって, $k$を大きくしていくと, $${y^k = \bm{v}_1}$$に成ります. $${(\bm{y}^k ,\bm{Ay}^k)}$$を求めれば, $${\lambda_1}$$も求まります
2. n個の固有値
前章では, 一番大きな固有値を求めました.
今, n個の固有値を求めたいとします.
ここでは, 方法の一つとして, 線型空間を限定します.
すなはち, $${\bm{v}_k}$$を求める時には, 常に, 考えるベクトルが, $${\bm{v}_1,…,\bm{v}_{k-1}}$$と直交するようにしておきます.
今回は$${\bm{x}}$$が$${\bm{v}_1,…,\bm{v}_{k-1}}$$と直交するためには
$$
\bm{x} = \bm{x} - \sum_{i=1}^{k-1} (\bm{x},\bm{v}_i)\bm{v}_i
$$
とすれば良いです. (本当は一回この操作をすれば, 良いのですが, 数値誤差などもあるので, 毎回, この作業をします.)
3. 実際のコード
# k個の固有値を求める関数
import numpy as np
def power_method(A,k):
if A.shape[0] != A.shape[1]:
print("Matrix is not a hermitian matrix")
return 0
n = A.shape[0]
if A.shape[0] < k:
print("Matrix is smaller than the number")
return 0
eps = 1e-8
lam_multiple = []
v_multiple = []
for i in range(k):
#ベクトルの初期化
v = np.zeros(n)
v[0] = 1
#直交化
for s in range(i):
v -= np.dot(v,v_multiple[s])*v_multiple[s]
T = 100
lam = 0
for j in range(T):
#更新
v_next = np.dot(A,v)
lam_next = np.dot(v_next,v_next)/np.dot(v_next,v)
#収束
if np.abs(lam_next - lam) < eps:
lam_multiple.append(lam)
v_multiple.append(v)
break
else:
#直交化
for s in range(i):
v_next -= np.dot(v_next,v_multiple[s])*v_multiple[s]
#正規化
v = v_next/np.sqrt(np.dot(v_next,v_next))
lam = lam_next
if j == T -1:
print("error")
return 0
return lam_multiple,v_multiple これを用いて実際に確認してみましょう.
np.random.seed(1)
A = np.random.rand(10,10)
A = np.dot(A.T,A)
lam_multiple,v_multiple = power_method(A,10)
for i in range(10):
print(lam_multiple[i],np.sort(LA.eig(A)[0])[9-i])とすると,
25.140071139463146 25.140071147995904
1.901219696015664 1.9012197232353796
1.625216786296971 1.6252167697648858
1.256183924521971 1.2561839575241638
1.1399416975318308 1.1399416539026153
0.6726468744662356 0.6726468715225401
0.30344342638193034 0.30344343181034394
0.1990857245214645 0.19908571381457515
0.024128617199257688 0.024128629377396238
0.013025928571743344 0.013025916022419945が得られました！
4 . より良いアルゴリズム
今回のアルゴリズムは$${c_1 = 0}$$の場合を無視しています.
おそらくそんなことは滅多に起きないのですが,
乱数などを何度か振って対応すれば良いのかなぁと思っています.
他にも ハウスホルダ変換やヤコビ法など色々あるようです. 暇な時に実装したいと思います.
この記事が参加している募集
この記事が気に入ったらサポートをしてみませんか？