ベクトル電卓プログラムを改良しました。Array クラスの演算子の意味を変えてしまうような乱暴なことはやめて、きちんと Matrix クラスと Vector クラスを定義しました。使い方は次の説明で大体分かると思います。機能はあまりありませんがスクリプトをできるだけ簡単な構造にしましたので、連立方程式の解法プログラムや固有値のプログラムを勉強するときの土台に利用できると思います。
使い方は次のようにします。まず irb を --noinspect オプションをつけて起動します。続いて、require 'vector.rb' として、irb に 'vector.rb' を読み込みます。これで、ベクトル電卓を使うことができるようになります。
irb(main):001:0> require 'vector.rb' true
行列を作るのには次のようにします。
irb(main):002:0> a = Matrix[[1,2],[3,4]] [[1, 2] [3, 4]]
ベクトルを作るには次のようにします。
irb(main):003:0> b = Vector[1, 2] [1, 2]
行列とベクトルの積は * で計算します。
irb(main):004:0> a * b [5, 11]
行列と行列の積や和、ベクトルの積や和、行列やベクトルとスカラーとの積も計算できます。
irb(main):005:0> a * a [[7, 10] [15, 22]] irb(main):006:0> b + b [2, 4] irb(main):007:0> a * 2 [[2, 4] [6, 8]]
行列の列と列の交換や演算は次のようにします。行列 a = [[1,2],[3,4]] の逆行列を計算してみましょう。まずこの行列と単位行列を並べた行列を作ります。
irb(main):008:0> c = Matrix[[1,2,1,0],[3,4,0,1]] [[1, 2, 1, 0] [3, 4, 0, 1]]
次に第2行の行ベクトルから、第1行の行ベクトルを3倍したものを引きます。
irb(main):009:0> c.rs!(2, 1, 3) [[1, 2, 1, 0] [0, -2, -3, 1]]
次に第1行の行ベクトルに、第2行の行ベクトルを足します。
irb(main):010:0> c.ra!(1,2) [[1, 0, -2, 1] [0, -2, -3, 1]]
第2行を2で割ります。
irb(main):011:0> c.rd!(2, -2) [[1, 0, -2, 1] [-0.0, 1.0, 1.5, -0.5]]
元の行列が単位行列になったので、残りの行列が逆行列です。これを d として a * d を計算して確認します。
irb(main):012:0> d = Matrix[[-2,1],[1.5, -0.5]] [[-2, 1] [1.5, -0.5]] irb(main):013:0> a * d [[1.0, 0.0] [0.0, 1.0]]
間違いないようです。
vector.rb は次のようになります。ここを左クリックするとファイルをダウンロードできます。
ソースを見てもらうと分かりますが、Matrix クラスも、Vector クラスもインスタンス変数としては二次元配列 @a をひとつ持っているだけです。ただし、Vector クラスの場合は @a は多行一列の配列、たとえば、[[1], [2], [3]]の形式でなければなりません。したがって、2次元配列を引数にとり2次元配列を戻り値にとるようなメソッドを Matrix クラスのクラスメソッドとして作れば、行列処理の機能の拡張は簡単です。新しい二次元配列のオブジェクトが必要なときは Matrix.dim( row, column )で簡単に作ることができます。またどちらのクラスのコンストラクターの initialize メソッドも引数に二次元配列を一個とるだけですから、Matrix.new( arry ) や Vector.new( arry )とすれば、簡単に Matrix や Vector のオブジェクトを作り出すことができます。
実際には Ruby に添付されている matrix.rb モジュールの使い方を勉強したほうが速いかも知れません。文書が乏しい状況では、他の人の作ったプログラムを解析するか、自分で作るか、選択が難しい所です。
# vector.rb version 0.1, 2002/4/21
# simple vector caliculator
# by Tomokiyo Nomura
#
# Usage: this program is supposed to be used on the irb with --noinspect option.
# after irb prompt appeared, type 'require "vector.rb"'. To create a Matrix,
# type as follows.
# irb> a = Matrix[[1,2],[3,4]]
# or to create a vector,
# irb> b = Vector[[1,2]]
# in this case, try 'a + a', 'a * a', 'a * 2', 'b + b', 'b * 2', 'a * b'
# on the irb prompt. To calculate inner_product, use Vector.inner_product
# class method. ex,
# irb> Vector.inner_product( b, b )
# To transpose matrix or vector, use transpose method.
# irb> a.transpose
# To convert a vector to the matrix, use to_matrix method.
# irb> b.to_matrix
#
# Each row or column of the matrix can be processed, to exchange rows,
# irb> a.row.exchange(1,2)
# The following methods are usable for row or column calculation.
# row_exchange!(i, j) has 'ce!' shortcut
# row_add!( i, j, scalar ) has 'ca!' shortcut
# row_subtract!( i, j, scalar ) has 'cs!' shortcut
# row_multiply!( i, scalar ) has 'cm!' shortcut
# row_divide!( i, scalar ) has 'cd!' shortcut
# column_exchange! has 'ce!' shortcut
# column_add!( i, j, scalar ) has 'ca!' shortcut
# column_subtract!( i, j, scalar ) has 'cs!' shortcut
# column_multiply!( i, scalar ) has 'cm!' shortcut
# column_divide!( i, scalar ) has 'cd!' shortcut
class Matrix
def initialize( a )
@a = a
end
attr_reader :a
def to_s
b = "["
@a.each {|c| b += c.inspect + "\n"}
b.chop + "]"
end
def Matrix.[]( *a )
Matrix.new( a )
end
def Matrix.dim( row, column )
(0...row).collect { Array.new(column).fill(0) }
end
def Matrix.add( x, y )
row, column = x.size, x[0].size
c = Matrix.dim(row, column)
for i in 0...row
for j in 0...column
c[i][j] = x[i][j] + y[i][j]
end
end
c
end
def Matrix.sub( x, y )
row, column = x.size, x[0].size
c = Matrix.dim(row, column)
for i in 0...row
for j in 0...column
c[i][j] = x[i][j] - y[i][j]
end
end
c
end
def Matrix.multiply( x, y )
x_row, x_column = x.size, x[0].size
y_column = y[0].size
c = Matrix.dim(x_row, y_column)
for i in 0...x_row
for j in 0...y_column
for k in 0...x_column
c[i][j] += x[i][k] * y[k][j]
end
end
end
c
end
def Matrix.multiply_by_scalar( x, a )
x_row, x_column = x.size, x[0].size
c = Matrix.dim(x_row, x_column)
for i in 0...x_row
for j in 0...x_column
c[i][j] = x[i][j] * a
end
end
c
end
def Matrix.transpose( arry )
row, column = arry.size, arry[0].size
c = Matrix.dim( column, row )
for i in 0...row
for j in 0...column
c[j][i] = arry[i][j]
end
end
c
end
def +( other )
Matrix.new( Matrix.add( self.a, other.a ) )
end
def -( other )
Matrix.new( Matrix.sub( self.a, other.a ) )
end
def *( other )
if other.kind_of?( Numeric )
Matrix.new( Matrix.multiply_by_scalar( self.a, other ) )
elsif other.type == Vector
Vector.new( Matrix.multiply( self.a, other.a ) )
else
Matrix.new( Matrix.multiply( self.a, other.a ) )
end
end
def row_exchange!( i, j )
i -= 1; j -= 1
@a[i], @a[j] = @a[j], @a[i]
self
end
def row_add!( i, j, c = 1 )
i -=1; j -= 1
@a[i].each_index {|k| @a[i][k] += @a[j][k] * c}
self
end
def row_subtract!( i, j, c = 1 )
i -=1; j -= 1
@a[i].each_index {|k| @a[i][k] -= @a[j][k] * c}
self
end
def row_multiply!( i, c )
i -=1; column = @a.size
(0...column).each {|k| @a[i][k] *= c}
self
end
def row_divide!( i, c )
i -=1; column = @a[0].size
(0...column).each {|k| @a[i][k] /= c.to_f}
self
end
def transpose
Matrix.new( Matrix.transpose( @a ) )
end
alias re! row_exchange!
alias ra! row_add!
alias rs! row_subtract!
alias rm! row_multiply!
alias rd! row_divide!
def column_exchange!( i, j )
row = @a.size
i -= 1; j -= 1
for k in 0...row
@a[k][i], @a[k][j] = @a[k][j], @a[k][i]
end
self
end
def column_add!( i, j, c = 1 )
row = @a.size
i -= 1; j -= 1
for k in 0...row
@a[k][i] += @a[k][j] * c
end
self
end
def column_subtract!( i, j, c = 1)
row = @a.size
i -= 1; j -= 1
for k in 0...row
@a[k][i] -= @a[k][j] * c
end
self
end
def column_multiply!( i, c )
row = @a.size
i -= 1
for k in 0...row
@a[k][i] *= c
end
self
end
def column_divide!( i, c )
row = @a.size
i -= 1
for k in 0...row
@a[k][i] /= c.to_f
end
self
end
alias ce! column_exchange!
alias ca! column_add!
alias cs! column_subtract!
alias cm! column_multiply!
alias cd! column_divide!
end
class Vector
def initialize( a )
@a = a
end
attr_reader :a
def to_s
@a.collect {|c| c[0]}.inspect
end
def Vector.[]( *a )
Vector.new( a.collect{|c| c.to_a} )
end
def +( other )
Vector.new( Matrix.add( self.a, other.a ) )
end
def -( other )
Vector.new( Matrix.sub( self.a, other.a ) )
end
def *( other )
if other.kind_of?( Numeric )
Vector.new( Matrix.multiply_by_scalar( self.a, other ) )
else
Vector.new( Matrix.multiply( self.a, other.a ) )
end
end
def to_matrix
Matrix.new( @a )
end
def Vector.inner_product( x, y )
p = 0
x.a.each_index {|i| p += x.a[i][0] * y.a[i][0]}
p
end
end
vector.rb の拡張法の例としてガウス・ジョルダン法で連立一次方程式を解くメソッドを作りました。下に示したスクリプトを gaussjordan.rb という名前で作成して irb から起動するだけです。実行例は次のようになります。
$ irb --noinspect irb(main):001:0> require 'gaussjordan.rb' true irb(main):002:0> a = Matrix[[1,2],[3,4]] [[1, 2] [3, 4]] irb(main):003:0> a.inv [[-2.0, 1.0] [1.5, -0.5]] irb(main):004:0> b = Vector[3, 7] [3, 7] irb(main):005:0> Matrix.solver( a, b ) [1.0, 1.0]
次のスクリプトを見れば分かるように、スクリプトの始めに require 'vector.rb' と記述するだけで vector.rb を部品のように組み込むことができます。Matrix クラスの拡張は gaussjordan.rb の class Matrix で定義できます。機能を拡張するのに vector.rb に一切手をいれないですみます。Ruby でプログラムを作るのは非常に楽です。
gaussjordan.rb のソースリスト
require 'vector.rb'
class Matrix
def Matrix.gaussjordan( arry )
matrix = []
arry.each{|c| matrix.push( c.dup )}
n = matrix.size
m = matrix[0].size
for k in 0...n
for i in k.succ...n
if matrix[i][k].to_f.abs > matrix[k][k].to_f.abs
matrix[i], matrix[k] = matrix[k], matrix[i]
end
end
coeff = matrix[k][k].to_f
for j in k.succ...m
matrix[k][j] /= coeff
end
for i in 0...n
if i == k
next
else
coeff = matrix[i][k].to_f
for j in k.succ...m
matrix[i][j] -= matrix[k][j] * coeff
end
end
end
end
matrix.collect{|c| c[n...m]}
end
def Matrix.unit( n )
matrix = Matrix.dim( n, n )
for i in 0...n
matrix[i][i] = 1
end
matrix
end
def Matrix.solver( matrix, vector )
arry = []
matrix.a.each{|c| arry.push( c.dup )}
v = vector.a
v.each_index {|i| arry[i].push( v[i][0] )}
Vector.new( Matrix.gaussjordan( arry ) )
end
def inv
n = @a.size
unit = Matrix.unit( n )
arry = @a.collect{|c| c + unit.shift}
Matrix.new( Matrix.gaussjordan( arry ) )
end
end
行列式の計算をするメソッドです。今のところ determinant.rb というファイルに作成しています。固有値と固有ベクトルを計算するメソッドまで完成したら vector.rb に組み込むつもりですが、気力がまだありません。使い方は次のようにします。
$ irb --noinspect irb(main):001:0> require 'determinant.rb' true irb(main):002:0> a = Matrix[[1,2],[3,4]] [[1, 2] [3, 4]] irb(main):003:0> a.det -2.0
determinat.rb のソースリスト
require 'vector.rb'
class Matrix
def determinant
matrix = []
@a.each{|c| matrix.push( c.dup )}
n = matrix.size
det = 1
for k in 0...n
if matrix[k][k] == 0
for i in k.succ...n
if matrix[i][k] != 0
matrix[i], matrix[k] = matrix[k], matrix[i]
det *= -1
break
end
end
end
akk = matrix[k][k].to_f
if akk == 0
return 0
else
det *= akk
end
for i in k.succ...n
coeff = matrix[i][k] / akk
for j in k.succ...n
matrix[i][j] -= matrix[k][j] * coeff
end
end
end
det
end
alias det determinant
end
ヤコビ法を使って対象行列の固有値と固有ベクトルを計算するメソッドです。行列 a に対し、eigenvalue というメッセージを送ると固有値のリストが帰って来ます。また、eigenvector というメッセージを送ると固有ベクトルのリストが帰って来ます。使い方は次のようにします。固有値計算プログラムのファイル名は jacobi.rb です。
$ irb --noinspect irb(main):001:0> require 'jacobi.rb' true irb(main):002:0> a = Matrix[[1,2],[2,2]] [[1, 2] [2, 2]] irb(main):003:0> b = a.eigenvalue -0.56155281283.561552813 irb(main):004:0> c = a.eigenvector [0.788205438, -0.6154122094][0.6154122094, 0.788205438] irb(main):005:0> d = b.shift -0.5615528128 irb(main):006:0> e = c.shift [0.788205438, -0.6154122094] irb(main):007:0> a * e [-0.4426189808, 0.3455864572] irb(main):008:0> e * d [-0.4426189808, 0.3455864572]
相関係数の行列のように対角線上の要素の値が同じ行列はこのプログラムでは正しい答えがでません。例えば a = Matrix[[1,2],[2,1]]では正しい答えがでませんが、a = Matrix[[1.000000001,2],[2,1]] のように対角線上の要素の値を僅かに変えてやると計算が可能になるようです。
jacobi.rb のソースは次のようになります。
require 'vector.rb'
class Matrix
Iterate_Max = 50
Epsilon = 1e-15
def Matrix.jacobi( arry )
a = arry.collect{|c| c.dup }
n = a.size
w = (0...n).collect{ Array.new(n).fill(0) }
w.each_index{|i| w[i][i] = 1}
for iterate in 0..Iterate_Max
amax = 0
for i in 0...n
for j in i.succ...n
if amax < a[i][j].to_f.abs
amax = a[i][j].to_f.abs
imax = i
jmax = j
end
end
end
if amax < Epsilon
break
end
alpha = (a[imax][imax] - a[jmax][jmax])/2.0
beta = Math.sqrt( alpha * alpha + amax * amax )
sgn = (a[imax][imax] - a[jmax][jmax]) <=> 0
cos = Math.sqrt( ( 1 + sgn * alpha / beta) / 2 )
sin = sgn * amax / 2.0 / beta / cos
for j in 0...n
x, y = a[imax][j], a[jmax][j]
a[imax][j] = x * cos + y * sin
a[jmax][j] = -1 * x * sin + y * cos
end
for i in 0...n
x, y = a[i][imax], a[i][jmax]
a[i][imax] = x * cos + y * sin
a[i][jmax] = -1 * x * sin + y * cos
end
for j in 0...n
x, y = w[imax][j], w[jmax][j]
w[imax][j] = x * cos + y * sin
w[jmax][j] = -1 * x * sin + y * cos
end
end
[ a, w ]
end
def eigenvalue
a, w = Matrix.jacobi( @a )
n = a.size
(0...n).collect{|i| a[i][i] }
end
def eigenvector
a, w = Matrix.jacobi( @a )
n = w.size
ret_val = []
for i in 0...n
v = (0...n).collect{ w[i].shift.to_a }
ret_val.push( Vector.new( v ) )
end
ret_val
end
end
ここで説明したスクリプトを一本にまとめたスクリプト vector2.rb もダウンロードできます。ここから左クリックでダウンロードしてください。(2002/04/29)