再帰のLU分解

資源の無駄使い

lu_decompの返り値が(転置されたL, U)なのであとから転置する。

# coding: utf8
# 正方行列のLU分解
# 再帰で書いたLU分解

def lu_decomp(A):
    if len(A)==1:
        return [ [1] ], [ [A[0][0]] ]
    else:
        li = [row[0] / A[0][0] for row in A]
        ui = A[0][:] # copy
        for i in xrange(len(li)):
            for j in range(len(ui)):
                A[i][j] -= li[i] * ui[j]
        nextA = [row[1:] for row in A[1:]]
        nextL, nextU = lu_decomp(nextA)
        L = [li] + [[0]+row for row in nextL]
        U = [ui] + [[0]+row for row in nextU]
        return L, U

def transpose(m):
    for i in xrange(len(m)):
        for j in xrange(i, len(m[0])):
            m[i][j], m[j][i] = m[j][i], m[i][j]
    return m

def main():
    import copy
    A = [[3.0, 1.0, -2.0],
        [8.0, 4.0, 6.0],
        [5.0, 1.0, 1.0],
        ]
    L,U = lu_decomp(copy.deepcopy(A))
    L = transpose(L)
    
    # print
    ms = [A,L,U]
    rep = ["A", "L", "U"]
    for k in range(3):
        m = ms[k]
        print rep[k],":"
        for i in xrange(len(m)):
            for j in xrange(len(m[0])):
                print "%3.1lf" % m[i][j],
            print ""
        print ""
    
    print "LU:"
    for i in xrange(3):
        for j in range(3):
            x = 0.0
            for k in range(3):
                x += L[i][k] * U[k][j]
            A[i][j] = x
    for i in range(3):
        for j in range(3):
            print A[i][j],
        print ""
    print ""
                

if __name__ == '__main__':
    main()

実行結果：

A :
3.0 1.0 -2.0 
8.0 4.0 6.0 
5.0 1.0 1.0 

L :
1.0 0.0 0.0 
2.7 1.0 0.0 
1.7 -0.5 1.0 

U :
3.0 1.0 -2.0 
0.0 1.3 11.3 
0.0 0.0 10.0 

LU:
3.0 1.0 -2.0 
8.0 4.0 6.0 
5.0 1.0 1.0