行列方程式を解く（２）

Python

昨日の問題設定で行列方程式Ax=bを解いてみる。
numpyモジュールのarrayで行列を作り、solveメソッドがあるのでとりあえずそれを使う。計算速度を計測する。読み込むデータはあらかじめ作成しておく。
ソースは以下の通り。

#                               (2009/01/15)
# solving matrix equation Ax=b
#
from numpy import *
from numpy.linalg import inv, solve
import time

bb = list()

# --- reading b data
input = open('BB', 'r')
while True:
    try:
       bb.append(float(input.readline().split(',')[1]))
    except:
        break
print "info> Data B has been read completely.\
 number of lines is", len(bb)
input.close()

# --- defining matrix list
matlen = len(bb)
aa = [[0 for j in xrange(matlen)] for i in xrange(matlen)]

# --- reading matrix data (diagonal)
input = open('DD', 'r')
while True:
    try:
        tmp = input.readline().split(',')
        i = int(tmp[0])
        aa[i-1][i-1] = float(tmp[1])
    except:
        break
print "info> Diagonal components have been read completely.\
 number of lines is", i
input.close()

# --- reading matrix data (upper triangular)
input = open('AU', 'r')
while True:
    try:
        tmp = input.readline().split(',')
        i = int(tmp[0])
        j = int(tmp[1])
        aa[i-1][j-1] = float(tmp[2])
    except:
        break
print "info> Upper-triangular components have been read completely."
input.close()

# --- reading matrix data (lower triangular)
input = open('AL', 'r')
while True:
    try:
        tmp = input.readline().split(',')
        i = int(tmp[0])
        j = int(tmp[1])
        aa[i-1][j-1] = float(tmp[2])
    except:
        break
print "info> Lower-triangular components have been read completely."
input.close()

begin_time = time.clock()

xx = array(aa)
yy = array(bb)
zz = solve(xx, yy)

print zz
print dot(xx, zz)
print yy

end_time = time.clock()

print " *** Calculation time [sec]:", end_time-begin_time

Windows XP（Intel Core2 CPU 1.86GHz, 2.00GB RAM)で複数回実行した場合、計算時間のばらつきは、0.017秒から0.02秒の間となった。

実行結果の例：

info> Data B has been read completely. number of lines is 48
info> Diagonal components have been read completely. number of lines is 48
info> Upper-triangular components have been read completely.
info> Lower-triangular components have been read completely.
[ 47.79181255  49.14570428  50.9394021   52.29329383  48.74925936
  50.10315109  51.89684891  53.25074064  49.70670617  51.0605979
  52.85429572  54.20818745  42.48047401  43.74845139  45.42176138
  46.68973876  43.39536763  44.663345    46.336655    47.60463237
  44.31026124  45.57823862  47.25154861  48.51952599  30.98626448
  32.03097226  33.39455966  34.43926744  31.77349852  32.8182063
  34.1817937   35.22650148  32.56073256  33.60544034  34.96902774
  36.01373552  12.66011312  13.20737948  13.8990035   14.44626986
  13.10692163  13.65418799  14.34581201  14.89307837  13.55373014
  14.1009965   14.79262052  15.33988688]
[ -3.  -4.  -5.  -6.  -4.  -5.  -6.  -7.  -5.  -6.  -7.  -8.  -4.  -5.  -6.
  -7.  -5.  -6.  -7.  -8.  -6.  -7.  -8.  -9.  -5.  -6.  -7.  -8.  -6.  -7.
  -8.  -9.  -7.  -8.  -9. -10.  -6.  -7.  -8.  -9.  -7.  -8.  -9. -10.  -8.
  -9. -10. -11.]
[ -3.  -4.  -5.  -6.  -4.  -5.  -6.  -7.  -5.  -6.  -7.  -8.  -4.  -5.  -6.
  -7.  -5.  -6.  -7.  -8.  -6.  -7.  -8.  -9.  -5.  -6.  -7.  -8.  -6.  -7.
  -8.  -9.  -7.  -8.  -9. -10.  -6.  -7.  -8.  -9.  -7.  -8.  -9. -10.  -8.
  -9. -10. -11.]
 *** Calculation time [sec]: 0.0184352468261

思ったより速いか。これをもう少し速くできるか検討する。

  • やること
    • 他の解法があるか検討
    • solveではどういうアルゴリズムを用いているか調べる