概要
少々ニッチな内容になりますが、
Lanczos アルゴリズム
について調査してみました。
Lanczos アルゴリズムとは
ここにあるように、
Lanczosアルゴリズムは、行列Aに対して低ランクへの分解を行うことが可能です。
この様にA=VTV* と展開できたあとは、Tの性質を利用して、
固有値と固有ベクトルを計算することができます。
Tは Tridiagonal matrix と呼ばれ、
以下の様な形をしています。
この形の行列の固有値、固有ベクトルの求め方は上のWikiにも書かれています。
こちらにも記載がありました。
Lanczos アルゴリズムのPythonプログラム
以下の様なlanczos.py プログラムを実装してみます。
また、比較としてnumpy を持ちいて固有値などを求める方法と見比べて、
たしかにlanczosアルゴリズムが動いていることを確かめてみます。
使用するライブラリ
- Lanczos によるイテレーションアルゴリズム（numpy)
- pylanczos ライブラリ（https://github.com/mrcdr/pylanczos）
- SVD(Singular Value Decomposition)
np.linalg.svd - Numpy の線形演算エンジン（
np.linalg.eigh)
main.py
import numpy as np
from scipy.sparse.linalg import eigs
from scipy.sparse import csc_matrix
from pylanczos import PyLanczos
from sklearn.decomposition import PCA
from lanczos import lanczos_algorithm
# Create a symmetric matrix (replace this with your own matrix)
matrix = np.array([[1.0, 2.0, 3.0],
[2.0, 4.0, 5.0],
[3.0, 5.0, 6.0]])
# Ensure the matrix is symmetric
matrix = (matrix + matrix.T) / 2.0
# Convert the matrix to a sparse matrix
sparse_matrix = csc_matrix(matrix)
# Number of eigenvalues and eigenvectors to compute
num_eigenvalues = 3
# Convert the sparse matrix to a dense matrix and use scipy.linalg.eigh
dense_matrix = sparse_matrix.toarray()
eigenvalues, eigenvectors = np.linalg.eigh(dense_matrix)
# 'eigenvalues' will contain the real eigenvalues
# 'eigenvectors' will contain the corresponding eigenvectors
## cov then eigen
#X = matrix - matrix.mean(axis=0)
#cov =np.cov(X,rowvar=False)
#l, v = np.linalg.eig(cov)
#print("\n\nNumpy.linalg.eig of cov================================")
#print("Eigenvalue: {}".format(l))
#print("Eigenvector:")
#print(v)
pca = PCA(n_components=2)
pca.fit(matrix)
print(f"singular_values_: {pca.singular_values_}")
print("\n\nNumpy.linalg.eigh================================")
print("Eigenvalue: {}".format(eigenvalues))
print("Eigenvector:")
print(eigenvectors)
U, s, V = np.linalg.svd(matrix, full_matrices=True)
print("\n\nNumpy.linalg.svd================================")
print(U)
print(s)
print(V)
engine = PyLanczos(matrix, True, 2) # Find 2 maximum eigenpairs
eigenvalues, eigenvectors = engine.run()
print("\n\nLanczos================================")
print("Eigenvalue: {}".format(eigenvalues))
print("Eigenvector:")
print(eigenvectors)
eigenvalues, eigenvectors, Q, T = lanczos_algorithm(matrix, 2)
print("\n\nLanczos scratch================================")
print("Eigenvalues:", eigenvalues)
print("Eigenvectors (in columns):")
print(eigenvectors)
Lanczos によるイテレーションアルゴリズム（numpy)
の中身(lanczos.py) は以下の通りです。
lanczos.py
import numpy as np
def lanczos_algorithm(A, k):
n = A.shape[0]
# Initialize variables
Q = np.zeros((n, k))
alpha = np.zeros(k)
beta = np.zeros(k-1)
# Random initial vector
Q[:, 0] = np.random.rand(n)
Q[:, 0] /= np.linalg.norm(Q[:, 0])
for i in range(k):
# Arnoldi iteration
if i == 0:
V = A @ Q[:, 0]
else:
V = A @ Q[:, i] - beta[i-1] * Q[:, i-1]
alpha[i] = np.dot(Q[:, i], V)
V -= alpha[i] * Q[:, i]
if i < k - 1:
beta[i] = np.linalg.norm(V)
if beta[i] < 1e-10:
break
Q[:, i+1] = V / beta[i]
# Tridiagonalize the result
T = np.diag(alpha) + np.diag(beta, 1) + np.diag(beta, -1)
# Compute eigenvalues and eigenvectors of T
eigenvalues, eigenvectors = np.linalg.eig(T)
return eigenvalues, eigenvectors, Q, T
# Example usage
if __name__ == "__main__":
# Create a random symmetric matrix for testing
n = 5
A = np.random.rand(n, n)
A = A + A.T # Make it symmetric
k = 3 # Number of Lanczos iterations
eigenvalues, eigenvectors, Q, T = lanczos_algorithm(A, k)
print("Eigenvalues:", eigenvalues)
print("Eigenvectors (in columns):")
print(eigenvectors)
print("\n\n")
print(A)
print("\n")
print("Q")
print(Q)
print("\n")
print("T")
print(T)
最終的にTについて np.linalg.eig(T) を使ってeigenvalues, eigenvectors
を計算していますが、ここで大事なのは、
もともとAの固有値、固有ベクトルを計算していたのですが、それが
Lanczosアルゴリズムを使って A = (L^t)TL
と変形したことで、
Tの固有値、固有ベクトルを計算する問題に置き換わったと言うことです。
実行結果
main.pyの実行結果は以下の様になります。
Numpy.linalg.eigh================================
Eigenvalue: [-0.51572947 0.17091519 11.34481428]
Eigenvector:
[[-0.73697623 0.59100905 -0.32798528]
[-0.32798528 -0.73697623 -0.59100905]
[ 0.59100905 0.32798528 -0.73697623]]
Numpy.linalg.svd================================
[[-0.32798528 -0.73697623 -0.59100905]
[-0.59100905 -0.32798528 0.73697623]
[-0.73697623 0.59100905 -0.32798528]]
[11.34481428 0.51572947 0.17091519]
[[-0.32798528 -0.59100905 -0.73697623]
[ 0.73697623 0.32798528 -0.59100905]
[-0.59100905 0.73697623 -0.32798528]]
Lanczos================================
Eigenvalue: [11.34481428 0.17091519]
Eigenvector:
[[ 0.32798528 0.59100905]
[ 0.59100905 -0.73697623]
[ 0.73697623 0.32798528]]
Lanczos scratch================================
Eigenvalues: [11.34274132 -0.03293642]
Eigenvectors (in columns):
[[ 0.89613077 -0.44379009]
[ 0.44379009 0.89613077]]
lanczos scratchの結果が少し微妙ですが、
その他のライブラリでは固有値の値が揃っていることがわかります。