Pythonによる「べき乗法（Power Method）」の実装 - 最大特異値・特異ベクトルを求める

はじめに

　理工学の分野において、行列の固有値と固有ベクトルを求める問題に直面することが多々あります。一般に、固有値や固有ベクトルは固有値分解という行列分解で求められますが、これは正方行列に対してのみ定義されています。特異値分解は、固有値分解を長方形行列に拡張したものであり、より適用範囲が広く、画像圧縮や自然言語処理における潜在意味解析などにも応用されています。特異値分解は比較的計算コストが高いため、例えば絶対値最大の特異値やその特異ベクトルを求めたいだけの場合は、他の手法を使った方が計算効率が高いです。
　今回はべき乗法（Power method）と呼ばれる数値計算アルゴリズムをPythonで実装し、行列の最大特異値および特異ベクトルを求める方法について書いていきます。

Pythonによる実装

　べき乗法の理論については、下記リンクの資料（P5〜）を参考にしました。

http://www.cs.yale.edu/homes/el327/datamining2013aFiles/07_singular_value_decomposition.pdf

import numpy as np

def power_method(A, iter_num=1):
    """
    Calculate the first singular vector/value of a target matrix based on the power method.
    Parameters
    ----------
    A : numpy array
        Target matrix
    iter_num : int
               Number of iterations

    Returns
    -------
    u : numpy array
        first left singular vector of A
    s : float
        first singular value of A
    v : numpy array
        first right singular vector of A
    """
    # set initial vector q
    q = np.random.normal(size=A.shape[1])
    q = q / np.linalg.norm(q)

    for i in range(iter_num):
        q = np.dot(np.dot(A.T, A), q)

    v = q / np.linalg.norm(q)
    Av = np.dot(A, v)
    s = np.linalg.norm(Av)
    u = Av / s

    return u, s, v

特異値分解との比較

　特異値分解の結果と一致しているか確かめます。Pythonで特異値分解を行うにはNumPyのlinalg.svgを使うとよいでしょう。
以下のように4×3行列を作成します。

A = np.array([[1,2,3],[4,5,6],[7,8,9],[10,11,12]])

power_method関数を使うと、

u, s, v = power_method(A, iter_num=1)

u
#-> array([-0.14091245, -0.34396479, -0.54701712, -0.75006945])
s
#-> 25.462398134435947
v
#-> array([-0.50388184, -0.57446659, -0.64505134])

NumPyで特異値分解を行うと

U, S, V = np.linalg.svd(A, full_matrices=False)

U
#-> array([[-0.14087668,  0.82471435,  0.54771037],
#         [-0.34394629,  0.42626394, -0.72755711],
#         [-0.54701591,  0.02781353, -0.1880169 ],
#         [-0.75008553, -0.37063688,  0.36786364]])
S
#-> array([  2.54624074e+01,   1.29066168e+00,   1.97614008e-15])
V
#-> array([[-0.50453315, -0.5745157 , -0.64449826],
#         [-0.76077568, -0.05714052,  0.64649464],
#         [-0.40824829,  0.81649658, -0.40824829]])

両者を比較すると、左特異ベクトル、特異値、右特異ベクトルがほぼ一致していることがわかります。