コレスキー分解をする関数

pythonで正値対称行列をコレスキー分解する関数。
例にあるように、リストのリストで行列を作成して使用する。

def cholesky_decomposition(matrix):
    size = len(matrix)
    result = [[0.0 for x in xrange(size)] for x in xrange(size)]
    
    for i in xrange(size):
        #calc non-diagonal
        for j in xrange(i):
            tmp = matrix[i][j]
            for k in xrange(j):
                tmp -= result[i][k] * result[j][k] 
            result[i][j] = tmp/result[j][j]
        
        #calc diagonal
        tmp = matrix[i][i]
        for k in xrange(i):
            tmp -= result[i][k]**2
        result[i][i] = tmp**(0.5)
    return result

if __name__ == '__main__':
    print cholesky_decomposition([[10,2,3],[2,100,5],[3,5,1000]])

実行すると、

[[3.1622776601683795, 0.0, 0.0], [0.63245553203367588, 9.979979959899719, 0.0], [0.94868329805051377, 0.44088264883091133, 31.605468237157314]]

として、リストのリスト形式で下３角行列が出力される。値はRのchol関数とチェック済。