2009年4月13日月曜日

九元連立方程式求解機

情報処理学会歴史特別委員会では, 本年3月2日に情報処理技術の貴重な物件に対して「情報処理技術遺産」を認定した. その詳細は学会のページを見て欲しい.

その遺産の中に, 九元連立方程式求解機というものがある. 要するに連立方程式を解くアナログ計算機である. アナログ計算機と聞くと, なとなく眉に唾をつけたくなるのが, 現代の計算機屋の人情だが, なかなか面白い面もある.

その1つの課題が, 精度の高め方だ. われわれはアナコンと聞くや, 精度は悪いと思い勝ちである. まずこの機械の目盛リは, マイクロメーターで設定し, 滑車にはボールベアリングを使うという気の使いようではあるが, それでも程度は知れていると思う.

この計算機の使い方を読むと, 逐次的に機械を使うことで, 精度を高めていくと書いてあったので, その理由を考えてみた.

2元連立方程式
a00x + a01y = c0
a10x + a11y = c1 (0)
を解きたいとする. そして上のアナコンで, x',y'が得られたとする. その時
a00x' + a01y' = c0'
a10x' + a11y' = c1' (1)
と書くと, x'=x, y'=y (つまり正確な値が得られた) なら c0'=c0, c1'=c1 である.
違うなら, (0)-(1)
a00(x-x') + a01(y-y') = c0 - c0'
a10(x-x') + a11(y-y') = c1 - c1' (2)
を作る.

右辺の定数を c0 - c0', c1 - c1'として, もう一度この連立方程式を解くと, 前のx'、y'の代りに, x-x', y-y'が得られる. これをそれぞれ前に得られた値, x', y'に加えると, x, y が得られるということで, これを繰り返すらしい.

手元に計算機があるから, さっそくやってみることにする.

添字が面倒なので, ここでは
ax+by=c
dx+ey=f
と書くと
x=(ce-bf)/(ae-bd), y=(af-cd)/(ae-bd)
と解ける. これにわざと誤差を入れてみる.

Schemeのプログラムはこう書いた.

(define (solv a b c d e f)
(let* ((det (- (* a e) (* b d))) ;分母の行列式
(x (* (/ (- (* c e) (* b f)) det) 1.1)) ;誤差を入れる
(y (* (/ (- (* a f) (* c d)) det) 0.9)) ;誤差を入れる
(c0 (+ (* a x) (* b y)))
(f0 (+ (* d x) (* e y))))
(list (list x y) a b (- c c0) d e (- f f0))))

c0とf0は上のc0'とc1'のこと.
結果は最初に新しいxとyのリスト. 続いて次に解くための引数をリストにして出す.

x=2, y=3と決め. a00=2, a01=1, a10=1, a11=2とすると,

c0=a00x+a01y=2*2+1*3=7
c1=a10x+a11y=1*2+2*3=8

だから

2x+y=7
x+2y=8

を解くことになる
(solv 2 1 7 1 2 8)=>
((2.2 2.7) 2 1 -.1 1 2 .4)

最初の近似解は x=2.2, y=2.7である. また
c0'-c0=-.1,
c1'-c1=.4.

以後次のように進行する.
(solv 2 1 -.1 1 2 .4)=>
((-.22 .27) 2 1 .07 1 2 .08)

(+ 2.2 -.22)=>1.98
(+ 2.7 .27)=>2.97

(solv 2 1 .07 1 2 .08)=>
((2.2e-2 2.7e-2) 2 1 -1.e-3 1 2 4.e-3)

(+ 1.98 2.2e-2)=>2.002
(+ 2.97 2.7e-2)=>2.997

(solv 2 1 -1.e-3 1 2 4.e-3)=>
((-2.2e-3 2.7e-3) 2 1 7.e-4 1 2 8.e-4)

(+ 2.002 -2.2e-3)=>1.9998
(+ 2.997 2.7e-3)=>2.9997

(solv 2 1 7.e-4 1 2 8.e-4)=>
((2.2e-4 2.7e-4) 2 1 -1.e-5 1 2 4.e-5)

(+ 1.9998 2.2e-4)=>2.00002
(+ 2.9997 2.7e-4)=>2.99997

従って
xは2.2 1.98 2.002 1.9998 2.00002
yは2.7 2.97 2.997 2.9997 2.99997
のように2と3に近づく.

我々は計算機パワーを大活用して, 一気に計算しようとするが, アナコンのような低精度の計算機を使っても, 工夫次第では有効な計算が出来ると改めて思う.

1 件のコメント:

Shitohichi さんのコメント...

こんにちは.いつも楽しく拝見させて頂いています.

ところでこの解法は gradient 法の一種だと思います.現在のデジタル計算機でも大規模な疎行列の solver ではこの方法に似た方法が利用されているようです.というのは,逆行列を求めるような直接法では O(n^3) の計算量(n はたとえば正方行列の次元)が通常必要なのに対し,このような iterative な方法は行列の乗算を繰り返すことで解を求めることができるからです.疎な行列に限れば,問題によっては疎行列の性質,(ほとんどの要素が 0 であること,non-zero の要素の位置があらじめわかっている場合など)を生かして,O(n) 程度の行列乗算関数を書くことができるためです.

ここで示された方法は多次元空間での勾配はあまり気にしていないようですが,もし共役な勾配を計算できる場合にはもっと収束を早くすることができると思います.たとえば,Numerical Recipe という本での Conjugate gradient methodが私にはわかりやすい解説でした.http://en.wikipedia.org/wiki/Conjugate_gradient の図(http://en.wikipedia.org/wiki/File:Conjugate_gradient_illustration.svg)には残差を単に小さくしようとした場合(緑)と,共役勾配方向に進んだ場合(赤)の収束の違いが示されています.

数学的基礎は同じなので,アナログ計算機で利用されている面白いアルゴリズムがあればデジタル計算機でも利用可能だと思うので,他にも面白いアルゴリスムがあれば楽しみにしています.