python で べき乗法

2016年1月8日

パソコンで固有値求めるのってどうやってるんかな?と思ったら、色々方法があるらしい。
古典的なものが「べき乗法」らしくて、基礎として理解しておくことは大事らしい。
なので、とりあえず べき乗法 – Wikipediaを見て、pythonでちょっとばかり実験。

# -*- coding:utf-8 -*-
import numpy as np
import sys


def power_method_tekitou(A, x, loop_count, a_err=sys.float_info.min, r_err=0.0):
    '''
    A,x,loop_count : Ax -> Ax^2 -> Ax^3 ... -> Ax^{loop_count}
    a_err: tollerable absolute error
    r_err: tollerable relative error
    '''
    q = [x, 0]  # q[1] = A * q[0]
    koyuuchi = [0, 0]  # lambda_{k+1}=kyuuchi[1] / lambda_{k}=koyuuchi[0]
    for i in xrange(loop_count):
        #power_method
        q[0] = q[0] / np.linalg.norm(q[0])
        q[1] = np.dot(A, q[0])

        koyuuchi[1] = np.dot(q[1].T, q[1]) / np.dot(q[1].T, q[0])

        # exit condition
        if np.abs(koyuuchi[0] - koyuuchi[1]) < a_err + r_err*(np.abs(koyuuchi[0]) + np.abs(koyuuchi[1])):
            return koyuuchi[1][0][0], i

        # next loop preparation
        q[0] = q[1]
        koyuuchi[0] = koyuuchi[1]

    return koyuuchi[1][0][0], i

if __name__ == "__main__":
    A = np.array([[2., 1., 1.], [0., 1., 2.], [3., 3., 1.]])
    x = np.array([[1.], [1.], [1.]])
    loop_count = 100

    ans, i = power_method_tekitou(A, x, loop_count)
    print "max eig(jibun):", ans, "loop_count:", i
    print "max eig(numpy):", np.max(np.linalg.eig(A)[0])

おっ、適当に書いただけだけどできた!この例だと30回程度で終了なんですね。本当はwhileで無限ループ作って回す方がよいのでしょうけれど…

「行列がどういう時にベクトルが収束しないのか?」の終了条件を数式からよく考えなきゃいけないので、for分で繰り返し回数を指定しちゃった。無限ループって怖いね

 

でも誤差判定いれるだけで一気に紛らわしくなりますね…自分で軽く実験するだけなら

# -*- coding:utf-8 -*-
import numpy as np


def power_method_tekitou(A, x, loop_count):
    '''
    q[1]=Ax -> q[1]=Ax^2 -> q[1]=Ax^3 ... -> q[1]=Ax^{loop_count}
    q[0]= x -> q[0]=Ax   -> q[0]=Ax^2 ... -> q[0]=Ax^{loop_count-1}
    '''
    q = [x, 0]
    for i in xrange(loop_count):
        # power_method
        q[0] = q[0] / np.linalg.norm(q[0])  # measure against overflow/underflow
        q[1] = np.dot(A, q[0])
        koyuuchi = np.dot(q[1].T, q[1]) / np.dot(q[1].T, q[0])

        # next loop preparation
        q[0] = q[1]

    return koyuuchi[0][0]

if __name__ == "__main__":
    A = np.array([[2., 1., 1.], [0., 1., 2.], [3., 3., 1.]])
    x = np.array([[1.], [1.], [1.]])
    loop_count = 100

    print "max eig(jibun):", power_method_tekitou(A, x, loop_count)
    print "max eig(numpy):", np.max(np.linalg.eig(A)[0])

これで十分かな?収束しても、loop_countで指定した回数だけ回っちゃうので大分無駄だけど、アルゴリズムという面だとこう書いた方がわかりやすかったですね・・。

while分で回したいけど、収束しない場合はどういう時だろ~。とりあえず最大固有値が0だと明らかにマズいですよね。あと、最大固有値が複数ある時とかだろうけれど、ジョルダン系とかだと固有ベクトル変わってくる? マズいんかな。

行列論はだいぶ専門的にやってたはずなのにそこらへん忘れちゃった\(^o^)/

numpy, scipy, matlplotlibが有名だと思ってたんですが、sympyやsageっていうのも有名みたいですね。こういうのって海外の大学教授が主導になって開発してる印象あるけど、日本の大学教授ってやらないのかな。ゼミでわけわかんない本読むより、ゼミの学生たちつかって勉強も兼ねてソース開発に貢献すればいいのにw

自分がモロに理系院生だった時代はかなり数学理論に頭回ったんだけど、ちょっと離れるともうこれだもん~


PAGE TOP