固有値から固有ベクトルを求める

8017 ワード
Python3 線形代数 Python3 テキストリンク

今話題の固有値から固有ベクトルを求める論文を実装してみた。
https://arxiv.org/abs/1908.03795

（Twitter 上でRで実装している人をすでに観測しているので、多分 $N$ 番煎じではあることを承知してこの記事を書いてる）

正確には $n$ 次 Hermitian Matrix $A$ に対し、固有ベクトルの各要素の大きさの2乗を求める。
そのため、複素数の位相成分はこの手法では求まらない。
また、Hermitian ではない正方行列にこの手法を適用しても、固有ベクトルの各要素の大きさは求まるとは限らない。
（Hermitian の固有ベクトルは直交するという性質を証明に使っているため）

$A$ の$i$ 番目の固有ベクトル$\mathbf{v}_i$と固有値 $\lambda_i$ と定義する。
また、$A$ の$j$ 行と$j$ 列を削除した行を $M_j$ とする。

M_j = 
\left[ \begin{array}{ccc|ccc}
A_{1,1} &\cdots& A_{1,j-1} & A_{1,j+1} &\cdots& A_{1,n} \\
\vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
A_{j-1,1} &\cdots& A_{j-1,j-1} & A_{j-1,j+1} &\cdots& A_{j-1,n} \\ \hline
A_{j+1,1} &\cdots& A_{j+1,j-1} & A_{j+1,j+1} &\cdots& A_{j+1,n} \\
\vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
A_{n,1} &\cdots& A_{n,j-1} & A_{n,j+1} &\cdots& A_{n,n}
\end{array} \right]

$A_{i,j}$ は$A$ の$i$ 行$j$ 番目の要素。
また、$M_j$ の$k$番目の固有値を$\mu^{(j)}_k$ と定義する。

すると$\mathbf{v}_i$ の$j$ 番目の要素$v_{i,j}$ に対して、次のような式が成立する（論文の式(2)にあたる)。

|v_{i,j}|^2 \prod_{k=1,k\neq i}^n (\lambda_i - \lambda_k)
= \prod_{k=1}^{n-1} (\lambda_i - \mu^{(j)}_k)

この式を使えば、固有値 $\lambda_i, \mu^{(j)}_k$ から固有ベクトルの各要素の大きさの2乗$|v_{i,j}|^2$ が求まる。

これを Python3 で実装してみた。

import numpy as np;
import matplotlib.pyplot as plt;


def calc_eigenvector (A, i):
    n,_ = A.shape;
    v = np.zeros ((n,1));   # 求める固有ベクトルの平方
    lams,_ = np.linalg.eig (A);

    for j in range (n):
        x = [xx for xx in range (n)];
        x.remove (j);
        Mj = A[:,x][x,:];
        muj, _ = np.linalg.eig (Mj);

        numerator = (lams[i] - muj).prod ();
        denominator = np.delete (lams[i]-lams, i).prod ();

        v[j,0] = numerator / denominator;

    return v;

N = 10;
A = np.random.randn (N,N);
A = A + A.T;
_, vs = np.linalg.eig (A);

plt.figure ();
plt.plot([0,1],[0,1]);
v_t = vs[:,0];
v_c = calc_eigenvector (A,0);

plt.scatter (v_t**2, v_c);
plt.show ();

calc_eigenvector が先の数式に基づいて固有ベクトルを求める関数である。
これを実行する。

横軸がnumpy.linalg.eig で求めた0番目の固有ベクトルの各要素の大きさの二乗、
縦軸が関数で求めた結果。
$y=x$ の直線上にちゃんと乗ってる。すごい。

全ての固有値でやったものが以下の図になる。

これもちゃんと$y=x$に乗ってることが確認できた。