作成日記　Schemeで数値計算(3) - 逆行列

今日は、逆行列を計算するプログラムを書いてみた。このページの後ろの方で、動作確認してみる。
まず逆行列を単純にガウスの消去法（＝掃き出し法）で計算するプログラム。

ページ最後の追記をご参照下さい。

(define (x-th-item x list-a) (if (null? list-a) (error "x-th-item") (if (= x 1) (car list-a) (x-th-item (- x 1) (cdr list-a))))) (define (extract-from-list x1 x2 list-a) (define (extract-from-list-p x1 x2 list-a tmp count) (cond ((> count x2) tmp) ((and (>= count x1) (<= count x2)) (extract-from-list-p x1 x2 (cdr list-a) (append tmp (list (car list-a))) (+ count 1))) ((< count x1) (extract-from-list-p x1 x2 (cdr list-a) tmp (+ count 1))))) (if (and (> x1 0) (>= x2 x1) (>= (length list-a) x2)) (extract-from-list-p x1 x2 list-a '() 1) (error "EXTRACT-FROM-LIST: (extract-from-list x1 x2 list-a) WRONG X1 X2!"))) (define (size-mat A) (define (N-mat A j) (if (null? (car A)) j (N-mat (list (cdr (car A))) (+ j 1)))) (define (M-mat A i) (if (null? A) i (M-mat (cdr A) (+ i 1)))) (cons (M-mat A 0) (N-mat A 0))) (define (x-th-scaled-list-to-each x list-a A) (define (x-th-scaled-list-to-each-p x list-a A tmp) (if (null? A) tmp (x-th-scaled-list-to-each-p x list-a (cdr A) (append tmp (scalar-mat (/ (x-th-item x (car A)) (x-th-item x list-a)) (list list-a)))))) (x-th-scaled-list-to-each-p x list-a A '())) (define (gyakugyouretu-hakidashi AE) (define (hakidashi-ahead i tmp-mat-up tmp-mat-low b-line size-mat-AE) (if (> i size-mat-AE) tmp-mat-up (begin (set! b-line (map (lambda (a) (/ a (x-th-item i (car tmp-mat-low)))) (car tmp-mat-low))) ;; normalize (set! tmp-mat-up (append tmp-mat-up (list b-line))) (if (not (null? (cdr tmp-mat-low))) ;;or (< i size-mat-AE) (set! tmp-mat-low (sub-mat (cdr tmp-mat-low) (x-th-scaled-list-to-each i b-line (cdr tmp-mat-low))))) ;; When i=N cdr \ tml is null. (hakidashi-ahead (+ i 1) tmp-mat-up tmp-mat-low b-line size-mat-AE)))) ;; i = from 1 to N (define (hakidashi-back i tmp-mat-up tmp-mat-low size-mat-AE) ;; When using hakidashi-back, input-matrix must be reve\ rsed . and output-matrix too. (if (> i size-mat-AE) (reverse tmp-mat-up) (begin (set! tmp-mat-up (append tmp-mat-up (list (car tmp-mat-low)))) (if (not (null? (cdr tmp-mat-low))) ;;or (< i size-mat-AE) (set! tmp-mat-low (sub-mat (cdr tmp-mat-low) (x-th-scaled-list-to-each (- (+ size-mat-AE 1) i) (car tmp-mat-low) (cdr tmp-mat-low))))) (hakidashi-back (+ i 1) tmp-mat-up tmp-mat-low size-mat-AE)))) (let ((N (car (size-mat AE)))) (hakidashi-back 1 '() (reverse (hakidashi-ahead 1 '() AE '() N)) N))) (define (connect-mat-side A B) (define (connect-mat-side-p A B tmp) (if (null? A) tmp (connect-mat-side-p (cdr A) (cdr B) (append tmp (list (append (car A) (car B))))))) (if (= (car (size-mat A)) (car (size-mat B))) (connect-mat-side-p A B '()) (error "CONNECT-MAT-SIDE: (connect-mat-side A B) WRONG SIZE MATRIXES !"))) (define (extract-mat A i j k l) (define (extract-mat-p D k l tmp) (if (null? D) tmp (extract-mat-p (cdr D) k l (append tmp (list (extract-from-list k l (car D))))))) (extract-mat-p (extract-from-list i j A) k l '())) (define (gyakugyouretu A) (let ((N (car (size-mat A)))) (extract-mat (gyakugyouretu-hakidashi (connect-mat-side A (eye-mat N))) 1 N (+ N 1) (* 2 N))))

Schemeでは、分数の計算が簡単にできる。Schemeは途中の演算で小数点を使わなければ自動的に分数計算する。以前のブログ
http://nibosiiwasi.hatenablog.com/entry/2016/08/04/231338
のソースを一部変更して、分数計算するようにする。

(define (sub-mat A1 A2) (add-mat A1 (scalar-mat -1 A2))) ;; 1.0 => 1に変更

(define (eye-line len loc1)
  (define (eye-line-p len loc1 tmp)
    (if (<= len 0)
        tmp
        (eye-line-p (- len 1)
                    (- loc1 1)
                    (append tmp (list (if (= loc1 1) 1 0))))))  ;; 1.0 => 1に変更
  (eye-line-p len loc1 '() ))

それでは、ためしに以前の日記
http://nibosiiwasi.hatenablog.com/entry/2016/08/20/221037
で書いたヒルベルト行列の逆行列を計算してみる。
hilbert-matrix-r はヒルベルト行列を生成する。

(define (hilbert-matrix-r n)
  (define (hilbert-matrix-r-p op ini count n)
    (if (= count n)
        '()
        (cons (op ini n)
              (hilbert-matrix-r-p op (+ ini 1) (+ count 1) n))))
  (define (hilbert-matrix-r-each-line ini n)
    (if (= n 0)
        '()
        (cons (/ 1 ini) (hilbert-matrix-r-each-line (+ ini 1) (- n 1)))))
  (hilbert-matrix-r-p hilbert-matrix-r-each-line 1 0 n))

実験として、20x20のヒルベルト行列の逆行列を分数を用いて計算してみる。結構大きいが、結果はすぐに返ってきた。（スペースを取るので実行結果は下に書く。）はやかった。自分用で使う分には、このレベルで十分強力な気がする。
逆行列の計算にかかった時間を印字してみた。

;(time (gyakugyouretu hil20))
; real   0.029
; user   0.030
; sys    0.000

逆行列の右下の方を見ると、かなり桁数が大きいので、メモリーはだいぶ食ってると思う。メモリーがどういう風に使われたのかどうすれば確認できるんだろう？

(define hil20 (hilbert-matrix-r 20))
hil20
=>
((1 1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20) (1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21) (1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22) (1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23) (1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24) (1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25) (1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26) (1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27) (1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28) (1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29) (1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30) (1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31) (1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32) (1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33) (1/15 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33 1/34) (1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33 1/34 1/35) (1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33 1/34 1/35 1/36) (1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33 1/34 1/35 1/36 1/37) (1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33 1/34 1/35 1/36 1/37 1/38) (1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30 1/31 1/32 1/33 1/34 1/35 1/36 1/37 1/38 1/39))

(define gh20 (gyakugyouretu hil20))

gh20
=>
((400 -79800 5266800 -171609900 3294910080 -41186376000 356948592000 -2237302782000 10440746316000 -37006645275600 100927214388000 -213323430411000 350069219136000 -444318624288000 431623806451200 -314725692204000 166619484108000 -60440401098000 13431200244000 -1378465288200) (-79800 21226800 -1576089900 54777880080 -1095557601600 14085740592000 -124619677182000 793496720016000 -3749272002075600 13423319513604000 -36914128662411000 78568660369836000 -129700645689888000 165464255684851200 -161454280100652000 118188754060608000 -62787775594698000 22846471615044000 -5091096452488200 523816809516000) (5266800 -1576089900 124826320080 -4519175106600 92965887907200 -1220177278782000 10966531592016000 -70700557753425600 337434480186804000 -1218166245859563000 3373383450072636000 -7222704706855638000 11984339661745651200 -15357151230750252000 15043739981143104000 -11050648504666848000 5888832426829044000 -2148710655394888200 480017665520316000 -49500688499262000) (-171609900 54777880080 -4519175106600 168285473017200 -3533994933361200 47119932444816000 -428791385247825600 2792314957736304000 -13438015734105963000 48851589962162988000 -136086572037454038000 292867300483909351200 -488112167473182252000 627936850324010304000 -617257652189248224000 454821427928919744000 -243045200549516488200 88904324471894316000 -19906187129228862000 2057028610969332000) (3294910080 -1095557601600 92965887907200 -3533994933361200 75391891911705600 -1017790540808025600 9355448405407104000 -61430929070198688000 297703733186347488000 -1088692576299632304000 3048339213638970451200 -6589514260887960402000 11025592488805999104000 -14233235274010900224000 14034489776092380672000 -10369928556779370163200 5555318869703234016000 -2036717251537942512000 456976817575340832000 -47311658052294636000) (-41186376000 14085740592000 -1220177278782000 47119932444816000 -1017790540808025600 13878961920109440000 -128637415574347680000 850582094818135680000 -4146587712238411440000 15241696068194852256000 -42867270191798021970000 93028436624300617440000 -156195893591418320640000 202261764420154897920000 -199991479309316424576000 148141836525419573760000 -79542065634387214320000 29222464913370479520000 -6569041752645524460000 681287875953042758400) (356948592000 -124619677182000 10966531592016000 -428791385247825600 9355448405407104000 -128637415574347680000 1200615878693911680000 -7986020779125829440000 39131501817716564256000 -144478577313097037010000 407939512413450457440000 -888364144801191698640000 1496192243875691281920000 -1942836614902487880576000 1925843874830454458880000 -1429793179798367704320000 769290895652382623520000 -283160273442983396460000 63763498612345890758400 -6623632127321249040000) (-2237302782000 793496720016000 -70700557753425600 2792314957736304000 -61430929070198688000 850582094818135680000 -7986020779125829440000 53392253209012688256000 -262789996263109325010000 974059652089259255520000 -2759835680919567890640000 6028662370412383621920000 -10181740892252025672576000 13254337256186068922880000 -13168270131145899644160000 9796684058920137899520000 -5281025000511636836460000 1947218076313218066758400 -439187362891157921040000 45689544061930248480000) (10440746316000 -3749272002075600 337434480186804000 -13438015734105963000 297703733186347488000 -4146587712238411440000 39131501817716564256000 -262789996263109325010000 1298491746241246076520000 -4829712441609243808620000 13726551149836798192920000 -30067953572431763314326000 50908704461260128362880000 -66422303749750640852160000 66127617397710930821760000 -49289566671441943806960000 26616366002578649655758400 -9829706635234995048540000 2220336112394187267480000 -231303316813521882930000) (-37006645275600 13423319513604000 -1218166245859563000 48851589962162988000 -1088692576299632304000 15241696068194852256000 -144478577313097037010000 974059652089259255520000 -4829712441609243808620000 18019628875711681578360000 -51355942295778292498326000 112776921688485978803880000 -191379018622885297364160000 250215634844739665861760000 -249577329653809207530480000 186351072808177541622758400 -100790844667884487656540000 37278274729144512543480000 -8431990712544592122930000 879523723192157283240000) (100927214388000 -36914128662411000 3373383450072636000 -136086572037454038000 3048339213638970451200 -42867270191798021970000 407939512413450457440000 -2759835680919567890640000 13726551149836798192920000 -51355942295778292498326000 146731263702223692852360000 -322952093926118939302020000 549174575178714331566720000 -719369950178626539352560000 718782709402970517687782400 -537551171562050600834880000 291173551262777408785560000 -107840723323596625572210000 24423697236336059942280000 -2550618797257256121396000) (-213323430411000 78568660369836000 -7222704706855638000 292867300483909351200 -6589514260887960402000 93028436624300617440000 -888364144801191698640000 6028662370412383621920000 -30067953572431763314326000 112776921688485978803880000 -322952093926118939302020000 712281693323218846365720000 -1213516958995113590104560000 1592364028560688696695782400 -1593613921832714355013440000 1193570920878823097170560000 -647405655744540206724210000 240083137818163818458280000 -54438604501155664185396000 5691463431896356634520000) (350069219136000 -129700645689888000 11984339661745651200 -488112167473182252000 11025592488805999104000 -156195893591418320640000 1496192243875691281920000 -10181740892252025672576000 50908704461260128362880000 -191379018622885297364160000 549174575178714331566720000 -1213516958995113590104560000 2071068943351660527111782400 -2721989792411433669565440000 2728162104865881364961280000 -2046121578649411023720960000 1111255684956145642193280000 -412587318324548191720896000 93657814195536626555520000 -9801964799377058648340000) (-444318624288000 165464255684851200 -15357151230750252000 627936850324010304000 -14233235274010900224000 202261764420154897920000 -1942836614902487880576000 13254337256186068922880000 -66422303749750640852160000 250215634844739665861760000 -719369950178626539352560000 1592364028560688696695782400 -2721989792411433669565440000 3582789983174023804385280000 -3595846797981509313876480000 2700329396185347990497280000 -1468304109175782969832896000 545758754722107075515520000 -124016939248194856280340000 12991953343553024480640000) (431623806451200 -161454280100652000 15043739981143104000 -617257652189248224000 14034489776092380672000 -199991479309316424576000 1925843874830454458880000 -13168270131145899644160000 66127617397710930821760000 -249577329653809207530480000 718782709402970517687782400 -1593613921832714355013440000 2728162104865881364961280000 -3595846797981509313876480000 3613560329006048768624640000 -2716862025141584814928896000 1478936989492394959739520000 -550283540316104136728340000 125167374677213361440640000 -13124524296038259424320000) (-314725692204000 118188754060608000 -11050648504666848000 454821427928919744000 -10369928556779370163200 148141836525419573760000 -1429793179798367704320000 9796684058920137899520000 -49289566671441943806960000 186351072808177541622758400 -537551171562050600834880000 1193570920878823097170560000 -2046121578649411023720960000 2700329396185347990497280000 -2716862025141584814928896000 2044949911396891796183040000 -1114337939999478146748180000 415028663403391672145280000 -94489096570053223832640000 9916307245895573787264000) (166619484108000 -62787775594698000 5888832426829044000 -243045200549516488200 5555318869703234016000 -79542065634387214320000 769290895652382623520000 -5281025000511636836460000 26616366002578649655758400 -100790844667884487656540000 291173551262777408785560000 -647405655744540206724210000 1111255684956145642193280000 -1468304109175782969832896000 1478936989492394959739520000 -1114337939999478146748180000 607820694545169898226280000 -226587340130160526888140000 51631542054350511594264000 -5422980525099141914910000) (-60440401098000 22846471615044000 -2148710655394888200 88904324471894316000 -2036717251537942512000 29222464913370479520000 -283160273442983396460000 1947218076313218066758400 -9829706635234995048540000 37278274729144512543480000 -107840723323596625572210000 240083137818163818458280000 -412587318324548191720896000 545758754722107075515520000 -550283540316104136728340000 415028663403391672145280000 -226587340130160526888140000 84541831107387625158264000 -19279944336904242362910000 2026580957476495940520000) (13431200244000 -5091096452488200 480017665520316000 -19906187129228862000 456976817575340832000 -6569041752645524460000 63763498612345890758400 -439187362891157921040000 2220336112394187267480000 -8431990712544592122930000 24423697236336059942280000 -54438604501155664185396000 93657814195536626555520000 -124016939248194856280340000 125167374677213361440640000 -94489096570053223832640000 51631542054350511594264000 -19279944336904242362910000 4400227536350517776520000 -462861082880434258020000) (-1378465288200 523816809516000 -49500688499262000 2057028610969332000 -47311658052294636000 681287875953042758400 -6623632127321249040000 45689544061930248480000 -231303316813521882930000 879523723192157283240000 -2550618797257256121396000 5691463431896356634520000 -9801964799377058648340000 12991953343553024480640000 -13124524296038259424320000 9916307245895573787264000 -5422980525099141914910000 2026580957476495940520000 -462861082880434258020000 48722219250572027160000))

(mul-mat  hil20 gh20)
=>
((1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0) (0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1))

追記 2016年9月6日

私のラップトップPCでは、40x40のヒルベルト行列の逆行列も計算できた。誤差無しであった。
この計算結果は、あまりにも長いので一番右下の(40,40)成分だけ書かさせて頂く。 58520505972825894943544903398826670092825440000。

ここまで計算結果の桁数が大きくても、『普通』に動作することに対して、心配になってきた。作業を一時中断することにした。
【お詫び】「Schemeで数値計算」中断します - nibosiiwasi’s blog

追記 2018年1月14日

ガウスの消去法（掃き出し法）の前進消去で、行の交換をしていませんでした。
ヒルベルト行列の逆行列計算は、それでもできましたが、
例えば、

'((0 1) (1 0))

のような行列の逆行列を計算できません。

前進消去で行の交換をするように訂正したコードはこちらのページに書きました。
nibosiiwasi.hatenablog.com
申し訳ありません。