注)Rubyにはちゃんと逆行列はき出す関数あるので，間違っても以下のcodeを実装なんかしないように．
Rubyで逆行列演算のプログラム．ソースコードは下の方に．
原理は単純に，元行列の各成分の余因子を行列式で割り，その転置行列を求める，というだけのもの． 念頭に“次元の制約をつけない，構造は出来るだけ単純に”というのがあったので， 計算速度はあまり速くないと思われる．もう少し綺麗なコードがかけないか，後でまた色々弄る予定．
で，どのくらいの速度で計算できるのか
day1 = Time.now
[main block]
day2 = Time.now
print day2 - day1以上のコードをプログラムの頭と末尾に追加，Check部分をコメントアウトして計算時間を計ってみた．以下それぞれの次数の計算速度．
| Dim | Time(s) | Dim | Time(s) | Dim | Time(s) |
|---|---|---|---|---|---|
| CPU Spec：Intel Celeron(Northwood-128K) 2.00GHz | |||||
| 4 | 0.016 | 6 | 0.375 | 8 | 26.453 |
| 5 | 0.047 | 7 | 2.969 | 9 | 256.672 |
で，これをプロット．
ごくごく単純に，アバウトに推定すれば10次で40分位かかるので，10次以上は計算しなかった． 果たして手計算で40分で10次の逆行列計算できるのか．誰かやってみて結果教えて＼（＞∀＜）／
以下，逆行列演算プログラム．
#inverse.rb（Shift_JIS：左クリックよりも右クリック→保存の方が良いかも）
def det(mt)
row = mt.size - 1
col = mt[0].size - 1
if(row == col)
val = 0
case row
when 0
val = mt[0][0]
when 1
val = mt[0][0] * mt[1][1] - mt[0][1] * mt[1][0]
else
sub_mtx = []
sub_mtx_col = []
sft = 0
for l in 0..row
for m in 0..(row - 1)
if l == m
sft = 1
end
for n in 0..(col - 1)
sub_mtx_col << mt[m + sft][n]
end
sub_mtx << sub_mtx_col
sub_mtx_col = []
end
sft = 0
val += (((-1)**(l + row)) * mt[l][row] * det(sub_mtx))
sub_mtx = []
end
end
else
val = 0
end
return(val)
end
def cof(mtb,m,n)
sub_mtx = []
sub_mtx_col = []
sub_r = mtb.size - 2
sub_c = mtb[0].size - 2
sft_r = 0
sft_c = 0
for x in 0..sub_r
if x == m
sft_r = 1
end
for y in 0..sub_c
if y == n
sft_c = 1
end
sub_mtx_col << mtb[x + sft_r][y + sft_c]
end
sft_c = 0
sub_mtx << sub_mtx_col
sub_mtx_col = []
end
return(((-1)**(m + n + 2)) * det(sub_mtx))
end
def transpose(mtt)
trs_mtx = []
trs_mtx_col = []
for a in 0..(mtt.size - 1)
for b in 0..(mtt[0].size - 1)
trs_mtx_col << mtt[b][a]
end
trs_mtx << trs_mtx_col
trs_mtx_col = []
end
return(trs_mtx)
end
#※行列サンプル（Encoding : Shift_JIS）
mtx = <適当な行列>
inv_mtx = []
inv_mtx_buf = []
inv_mtx_col = []
inv_r = mtx.size - 1
inv_c = mtx[0].size - 1
for m in 0..inv_r
for n in 0..inv_c
inv_mtx_col << cof(mtx,m,n)
end
inv_mtx_buf << inv_mtx_col
inv_mtx_col = []
end
inv_mtx = transpose(inv_mtx_buf)
######################## display part ##########################################
print("matrix:\\n")
for m in 0..(mtx.size - 1)
for n in 0..(mtx[0].size - 1)
print(mtx[m][n].to_s(10).rjust(8))
end
print "\\n"
end
print("\\ndet(matrix) : ")
print det(mtx)
print("\\n\\ninverse matrix(Elements isn't divided by determinant.) : \\n")
for m in 0..(inv_mtx.size - 1)
for n in 0..(inv_mtx[0].size - 1)
print(inv_mtx[m][n].to_s(10).rjust(8))
end
print "\\n"
end
# check
print("\\nCheck\\n")
for m in 0..(mtx.size - 1)
for n in 0..(inv_mtx[0].size - 1)
calc_buf = 0
for k in 0..(mtx[0].size - 1)
calc_buf += mtx[m][k] * inv_mtx[k][n]
end
print (calc_buf / det(mtx)).to_s(10).rjust(8)
end
print "\\n"
end