授業の課題をPythonで実装してみた

授業の課題で

が対称正定値行列での時、

minimize
subject to

を最急降下法、加速勾配法、共役勾配法のプログラムを実装して解くという課題が出ました。対称正定値行列Qとベクトルqはランダムに作ります。
対称正定値行列（固有値全てが正である行列）を作るのがなかなか面倒で、特殊な方法を使ったりすると一応作れるのですが、実際に固有値を求め、全ての固有値が正であることをチェックするし、条件を満たさなければ行列を作り直すという仕様にしました。

始めは近頃使っているRubyで実装しようと思いましたが、Rubyにはどうもデフォルトで固有値計算のライブラリがないようだったので、Pythonで実装することに。ちなみに研究室の同期の人たちはRやMATLABで実装していました。

まずは環境作り。IDEはRubyの時にインストールしたaptana studioにPyDevというEclipseプラグインが入っていたのでそれをそのまま使用。僕のMac LeopardのPythonはバージョン2.5.1でちょっと古かったので最新版Pythonをインストールしてみましたが、固有値計算をするライブラリであるNumpyがビルド出来なかったので、2.5.1を使うことに。

Python入門の本とネットの力を借りつつ作成。以下がコードです。

experiment = False #True when experiment
sumTime = 0
sumIter = 0
sumF = 0

from numpy import *
import time
#import random

def makeVector(n):
    '''make vector (size n)'''
    v = zeros(n)
    for i in range(n):
        #v[i] = random.random();
        v[i] = random.randint(0,100);
    return v

def makeSPDMatrix(n, eps):
    '''make symmetric positive definite matrix'''
    a = identity(n);
    while True:
        for i in range(n):
            for j in range(i,n):
                a[i,j] = a[j,i] = random.randint(0,100);
                
        a = dot(a,a)
    
        #l : engenvalue, v : eigenvector
        (l,v) = linalg.eig( a );
        #print "eigenvalue = ", l
        bOK = True;
        for i in range(n):
            if((0 - eps < l[i] and l[i] < 0 + eps) or l[i] < 0):
                bOK = False;
        if bOK:
            break
    return a

def doCGM(Q,q,init,eps):
    '''Conjugate gradient method'''
    if(experiment == False): print "----- Conjugate Gradient Method -----"
    x = init
    p = -gradf(Q,q,x)
    k = 0
    sttime = time.clock()
    while True:
        a = dot(-gradf(Q,q,x),p) / dot(dot(p,Q), p) #exact line search
        x1 = x + a*p
        #print "x = ", x, ". x1 = ", x1
        if(diffVector(x,x1) < eps): break
        gf1 = gradf(Q, q, x1)
        gf = gradf(Q, q, x)
        #beta = dot(gf1, gf1) / dot(gf, gf) #Fletcher-Reeves
        beta = dot(gf1, gf1 - gf) / dot(p, gf1 - gf) #Hestenes-Stiefel
        #print beta
        p = -gf1 + beta * p
        x = x1
        k += 1
    entime = time.clock()
    if(experiment):
        global sumTime,sumIter, sumF
        sumTime += entime - sttime
        sumIter += k
        sumF += f(Q, q, x)
    else:
        print "  iteration = ", k, ", time = ", entime - sttime, ", x = ", x, ", f(x) = ", f(Q, q, x)

def doAGM(Q, q, init, eps):
    '''accelerated gradient method'''
    if(experiment == False): print "----- Accelerated Gradient Method -----"
    y = init
    k = 0
    (l,v) = linalg.eig( Q );
    L = max(l)
    #print "l = ", l, ", L = ", L
    ssum = zeros(q.size)
    sttime = time.clock()
    while True:
        gf = gradf(Q,q,y)
        x = y - gf / L
        ssum += gf
        z = init - ssum / L
        y1 = (2.0/(k+3))*z + ((k+1)*1.0/(k+3))*x
        #print "y = ", y, "y1 = ", y1
        if(diffVector(y,y1) < eps): break
        y = y1
        k += 1
    entime = time.clock()
    if(experiment):
        global sumTime,sumIter, sumF
        sumTime += entime - sttime
        sumIter += k
        sumF += f(Q, q, x)
    else:
        print "  iteration = ", k, ", time = ", entime - sttime, ", x = ", x, ", f(x) = ", f(Q, q, x)

def doSDM(Q, q, init, eps):
    '''Steepest descent method
    init: initial value, eps: tolerance'''
    if(experiment == False): print "----- Steepest Descent Method -----"
    x = init
    k = 0
    sttime = time.clock()
    while True:
        d = -gradf(Q,q,x)
        a = dot(d,d) / dot(dot(d,Q), d) #exact line search
        x1 = x + a*d
        if(diffVector(x,x1) < eps): break
        #print "doSDM : k = ", k, "x = ",x, "x1 = ", x1, "a = ", a, "d = ", d
        x = x1
        k += 1
    entime = time.clock()
    if(experiment):
        global sumTime,sumIter, sumF
        sumTime += entime - sttime
        sumIter += k
        sumF += f(Q, q, x)
    else:
        print "  iteration = ", k, ", time = ", entime - sttime, ", x = ", x, ", f(x) = ", f(Q, q, x)
        
def diffVector(v1, v2):
    '''calculate Euclid norm'''
    return sqrt(sum([x*x for x in v1 - v2]))

def f(Q,q,x):
    '''objective function'''
    return dot(dot(x, Q), x) / 2 + dot(q, x)

def gradf(Q,q,x):
    '''gradient of f'''
    return dot(Q, x) + q

def exp(n, Q, q, init, eps):
    '''funtion for experiment
        n : loop iteration'''
    global experiment
    experiment = True
    ary = [doSDM, doAGM, doCGM]
    for i in range(len(ary)):
        global sumTime, sumIter, sumF
        sumTime = 0
        sumIter = 0
        sumF = 0
        print ary[i], "start."
        for j in range(n):
            ary[i](Q, q, init, eps)
        print "  average time = ", sumTime / n, ", average iteration = ", sumIter / n
    
    experiment = False

def doIt(n):
    eps = 1.0e-9;
    #n = 10 #size of matrix, vector
    it = 1 #num of iteration when experiment
    bExperiment = True #True if experient
    print "make start."
    Q = makeSPDMatrix(n, eps)
    print "made."
    q = makeVector(n)
    init = makeVector(n)
    print "n = ", n
    print "Q = %s" % Q
    print "q = %s" % q
    print "init = %s" % init
    invQ = linalg.inv(Q)
    trueX = -dot(invQ, q)
    ans = -(dot(dot(q, invQ), q)) / 2
    print "trueX = ", trueX, "true ans = %s" % ans
    if(bExperiment):
        exp(it, Q, q, init, eps)
    else:
        doSDM(Q, q, init, eps)
        doAGM(Q, q, init, eps)
        doCGM(Q, q, init, eps)

def doIt2():
    '''exec doIt func 5 times.'''
    print "doIt2 start."
    num = [10]
    for i in range(len(num)):
        for j in range(1):
            doIt(num[i])
    print "doIt2 end."

if __name__ == '__main__':
    doIt2()

Pythonでの実装の際にいくつか気になった点を書いておきます。

・コメントについて
1行コメントは#、複数行コメントは「'''」で囲みます。

・グローバル変数について
普通に変数名を書けばOK。関数などからの参照は普通に変数名で出来るが、グローバル変数への書き込みは出来ない。グローバル変数に値を書き込むときには

global g_Valuable
g_Valuable = 2

のようにグローバル変数の使用を宣言してから書き込む。

・ライブラリのインポートについて
ライブラリをインポートして関数を使う際に、

import (ライブラリ名) #import
（ライブラリ名）.（関数名）#use function

で使う方法と、

from (ライブラリ名) import (関数名) #import
（関数名）#use function

で使う方法があります。下の方法では関数名に*を用いるとライブラリの全ての関数をインポートできます。ただ、*を使ったインポートはあまり推奨されていないようで、小さいプログラム以外では避けた方がいいようです。

・関数について
定義は以下のように書きます。

def （関数名）(引数1, 引数2,...):
（関数本体）
return （戻り値）

関数の始めに簡単な関数のコメントを書くのが推奨されているみたいです。複数の値を返せるとか・・・
また、関数型言語のように関数を配列に入れたり、関数の引数や戻り値に出来たりします。これを高階関数というらしいです。便利です。

・ブロックについて
上の関数定義やif文などのブロックの始まりでは:を付ける。if文やfor文でよく付け忘れて怒られましたorz
また、pythonはブロックの範囲をインデントで判断します。インデントをずらすとまずいことに。

・構文糖衣について
複数の代入や、数値を簡単に取り扱えます。

eps = 1.0e-9; #eps = 10の-9乗
(a, b) = [1,2] #a = 1, b = 2

僕の印象ですが、pythonはシンタックス、能力とRubyのような使い心地でした。ある人に言わせればPythonとRubyは哲学が全然違うらしいですが・・・（RubyはTIMTOWTDI（there's more than one way to do it - あることをするのにいくつものやり方がある）、Pythonは対照的に言語自身の機能をなるべく抑えている）。
前にRuby使いの人が「Pythonはスクリプト言語の癖に速い」と言っていたので、スクリプト言語の中ではコンテスト向けかもしれません。