Ruby における数値厳密対角化プログラムの複素行列への対応

Ruby 標準の Matrix クラスは "require 'matrix'" と書くだけで使える様になる便利な代物である。また Matrix クラスのメソッド -eigensystem を呼び出すことで固有対角化を行うこともできる。ところが簡単だからといって LAPACK を呼び出す手間を惜しんで固有対角化を行った際、 Matrix が複素数を係数体にとっていた場合は、以下の様なエラーメッセージを目にすることがある。

/usr/lib/ruby/2.3.0/matrix/eigenvalue_decomposition.rb:379:in `to_f': can't convert <Complex String> into Float (RangeError)
	from /usr/lib/ruby/2.3.0/matrix/eigenvalue_decomposition.rb:379:in `sqrt'
	from /usr/lib/ruby/2.3.0/matrix/eigenvalue_decomposition.rb:379:in `block in reduce_to_hessenberg'
	from /usr/lib/ruby/2.3.0/matrix/eigenvalue_decomposition.rb:362:in `upto'
	from /usr/lib/ruby/2.3.0/matrix/eigenvalue_decomposition.rb:362:in `reduce_to_hessenberg'
	from /usr/lib/ruby/2.3.0/matrix/eigenvalue_decomposition.rb:36:in `initialize'
	from /usr/lib/ruby/2.3.0/matrix.rb:1372:in `new'
	from /usr/lib/ruby/2.3.0/matrix.rb:1372:in `eigensystem'
	from ... `<main>'

但し Mac の場合は "/System/Library/Frameworks/Ruby.framework/Versions/2.3" などがパスの接頭辞としてつく。また次の様なメッセージもしばしば出力される。

/usr/lib/ruby/2.3.0/matrix/eigenvalue_decomposition.rb:380:in `block in reduce_to_hessenberg': undefined method `>' for <Complex String>:Complex (NoMethodError)
	from /usr/lib/ruby/2.3.0/matrix/eigenvalue_decomposition.rb:362:in `upto'
	from /usr/lib/ruby/2.3.0/matrix/eigenvalue_decomposition.rb:362:in `reduce_to_hessenberg'
	from /usr/lib/ruby/2.3.0/matrix/eigenvalue_decomposition.rb:36:in `initialize'
	from /usr/lib/ruby/2.3.0/matrix.rb:1372:in `new'
	from /usr/lib/ruby/2.3.0/matrix.rb:1372:in `eigensystem'
	from ... `<main>'

/usr/lib/ruby/2.3.0/matrix/eigenvalue_decomposition.rb:519:in `hessenberg_to_real_schur': undefined method `>=' for <Complex String>:Complex (NoMethodError)
	from /usr/lib/ruby/2.3.0/matrix/eigenvalue_decomposition.rb:37:in `initialize'
	from /usr/lib/ruby/2.3.0/matrix.rb:1372:in `new'
	from /usr/lib/ruby/2.3.0/matrix.rb:1372:in `eigensystem'
	from ... `<main>'

いずれも Hessenberg 行列や Shur 分解におけるエラーメッセージの様であり、実態は Complex モジュールの実装上の不備である。以上のエラーメッセージを回避するには如何にすべきであろうか。

結論から言えば、以上のエラーメッセージは "require 'matrix'" の直後に以下の通り追記することで簡単に制御することができる。

module Math
	alias :_sqrt :sqrt
	def self.sqrt(number)
		return number.kind_of?(Complex) ? number ** 0.5 : _sqrt(number)
	end
end
class Complex
	def >(z)
		r1 = real
		r2 = z.real
		return r1 > r2 ? true: ((r1 < r2) ? false: imag > z.imag)
	end
	def >=(z) return !(self < z) end
	def <(z)
		r1 = real
		r2 = z.real
		return r1 < r2 ? true: ((r1 > r2) ? false: imag < z.imag)
	end
	def <=(z) return !(self < z) end
end

Math.sqrt に関する不具合は Complex の冪乗計算で代用して計算でき、Complex でないものは旧来の Math.sqrt を呼び出すことで回避する。冪乗の速度に不満がある場合は Complex.polar(_sqrt(number.abs), number.angle*0.5) とすれば良いかもしれないが、呼び出し回数が多い為、どちらが高速であるかは未確認である。また複素数の根は一般には多価函数であるが、位相の取り方は実装に依存している（"Complex(1,0) ** (1.0/3)" は "1" のままで "x^2+x+1=0" の根とはならない）。

Complex の順序に関しては、一般には複素数に順序を定めることができないと信じられているが、実際には複素数は順序構造そのものは導入できる代数の為、解決できる場合がある。例えば Quick Sort などのアルゴリズムを採用する際に数の大小を定めねばならないことは多いが、実数の様に二乗すれば 0 以上であるといった定理は必要ではないため、自然な順序は不必要であり、順序構造を保証する定義を導入することで急場を凌ぐことが出来る。一般には上記の辞書的な定義をとれば任意の複素数の順序を決める事が出来る。以上の点は複素数と順序構造、全順序などを参考にされたい。Ruby で採用しているアルゴリズムの実装において大小比較がどの様に用いられているかは確認していないので、計算値が正しいかどうかは LAPACK などのパッケージで各自確認されたい。

Ruby における数値厳密対角化プログラムの複素行列への対応

解決策