自作数値計算プログラム ほげ

私はmatlab環境がないため、scheme数値計算環境を作り始めた。

一方、世間にはmatlab以外にもscilabpythonのnumpy,scipy等
すぐれた数値計算ライブラリがある。それらは、すぐれている。
それらを使わずに、数値計算環境を自作するのは、とてつもなく非生産的だ
という自覚はある。

初歩的で恥ずかしいのだが、とりあえず、できたところまで
schemeのコードをブログに書くことにした。
今日のブログも、主な目的は私の備忘録である。

;;行列を表示。デバッグ用に使用。
(define (disp-mat a) 
  (begin (newline) (for-each (lambda (x) (display x) (newline)) a) (newline)))


(define (my-map-item2 op a b)
  (define (my-map-item2-p op a b tmp)
    (if (and (null? a) (null? b))
      tmp
      (my-map-item2-p op
                      (cdr a)
                      (cdr b)
                      (append tmp (list (op (car a) (car b)))))))
  (my-map-item2-p op a b '() ))

(define (l-accum2 op1 op2 ini a b)
  (define (l-accum2-p op1 op2 ini a b)
    (if (and (null? a) (null? b))
      ini
      (l-accum2 op1 op2 (op1 ini
                             (op2 (car a) (car b)))
                (cdr a) (cdr b))))
  (l-accum2-p op1 op2 ini a b))

(define (do-n-times op ini n)
  (define (do-n-times-p op ini n)
    (if (= n 0)
      ini
      (do-n-times-p op (op ini) (- n 1))))
  (do-n-times-p op ini n))

(define (find-1stitem-from-list filter? list-a)
  (define (find-1stitem-from-list-p filter? counter list-a)
    (cond ((null? list-a) #f) ;; 見つからなかったときはfalseを出力
          ((filter? (car list-a))
           (cons counter (car list-a))) ;;(位置. 非ゼロの数)を出力
          (else (find-1stitem-from-list-p filter? (+ counter 1) (cdr list-a)))))
  (find-1stitem-from-list-p filter? 1 list-a)) ;;先頭の位置を1としている

(define (zero-line n)
  (do-n-times (lambda (x) (append x '(0))) '() n))

;;mXnの零行列
(define (zero-mat m n)
  (do-n-times (lambda (x) (append x (list (zero-line n)))) '() m))

(define (do-n-times-with-filter-counter op ini n filter? op-c counter)
  (define (do-n-times-with-filter-counter-p op ini n filter? op-c counter)
    (if (= n 0)
      ini
      (do-n-times-with-filter-counter-p op
                                        (op ini filter? counter)
                                        (- n 1)
                                        filter?
                                        op-c
                                        (op-c counter))))
  (do-n-times-with-filter-counter-p op ini n filter? op-c counter))

(define (eye-line n loc1)
  (do-n-times-with-filter-counter (lambda (ini filter? counter)
                                    (if (filter? counter)
                                      (append ini '(1))
                                      (append ini '(0))))
                                  '()
                                  n
                                  (lambda (counter) (= counter 1))
                                  (lambda (counter) (- counter 1))
                                  loc1))

;;単位行列
(define (eye-mat n)
  (do-n-times-with-filter-counter (lambda (ini filter counter)
                                    (append ini (list (filter counter))))
                                  '()
                                  n
                                  (lambda (counter) (eye-line n counter))
                                  (lambda (counter) (+ counter 1))
                                  1))

;;和
(define (add-mat a1 a2)
  (my-map-item2 (lambda (s t) (my-map-item2 + s t)) a1 a2))
;;スカラー倍
(define (scalar-mat a AA)
  (define (scalar-each a a1) (map (lambda (x) (* a x)) a1))
  (map (lambda (x) (scalar-each a x)) AA))
;;引き算
(define (sub-mat A1 A2) (add-mat A1 (scalar-mat -1 A2)))

(define (sum-of-mul-each a1 a2) (l-accum2 + * 0 a1 a2))
;;転置行列
(define (transposed-mat a) (apply map list a))
;;行列の積
(define (mul-mat A1 A2)
  (define (mul-oneline-mat oneline A2)
    (map (lambda (x) (sum-of-mul-each oneline x)) (transposed-mat A2)))
  (define (mul-mat-pre A1 A2 tmp)
    (if (null? A1)
      tmp
      (mul-mat-pre (cdr A1) A2
                   (append tmp
                           (list (mul-oneline-mat (car A1) A2))))))
  (mul-mat-pre A1 A2 '() ))

;;行列のサイズ
;;(行数. 列数)
(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 (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-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 (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 '()))

;;リストのx番目の要素
(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 (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 '()))

;;cdrのdをn個繰り返す
(define (cdr*n x n)
  (define (cdr*n-p x n)
    (if (= n 0)
      x
      (cdr*n-p (cdr x) (- n 1))))
  (cdr*n-p x n))

;;nondiag-biggest-abs
;;入力:対称行列A
;;非対角成分の絶対値が最大となる成分をa_{ij}としたとき
;;出力:(a_{ij} i . j)
(define (nondiag-biggest-abs A)
  ;;内部定義
  (define (lbai l r)
    (define (lbai-p l tmp r c)
      (if (null? l)
        tmp
        (lbai-p (cdr l)
                (if (> (abs (car l)) (abs (car tmp)))
                  (cons (car l) (cons r c))
                  tmp)
                r
                (+ c 1))))
    (lbai-p l (cons 0 0) r (+ r 1)))
  
  (define (nondiag-biggest-abs-p A tmp row)
    (if (= (car (size-mat A)) 1)
      tmp
      (nondiag-biggest-abs-p (cdr A)
                             (if (> (abs (car (lbai (cdr*n (car A) row) row)))
                                    (abs (car tmp)))
                               (lbai (cdr*n (car A) row) row)
                               tmp)
                             (+ row 1))))
  
  (nondiag-biggest-abs-p A (cons 0 "all-nondiag") 1))

;;ヤコビ法
;;入力:実対称行列A, 正の実数tolerance
;;出力:(対角化された行列 . 固有ベクトルを並べた行列)
(define (jacobi-method A tolerance)
  ;;内部定義
  (define (norm-plane a b)
    (sqrt (+ (* a a) (* b b))))
  
  ;;ギブンスの回転行列(準備)
  ;;入力:行列A
  ;;出力: '(k m cos(psi) sin(psi))
  ;;ただしk,m は、いずれも回転する行と列
  (define (Pdat-from-mat A)
    (let ((nba (nondiag-biggest-abs A)))
      (let ((k (cadr nba))
            (m (cddr nba))
            (akm (car nba)))
        (let ((akk (caar (extract-mat A k k k k)))
              (amm (caar (extract-mat A m m m m))))
          (let ((s (/ (- akk amm) 2.0)))
            (if (>= s 0)
              (let ((c-psi (sqrt (/ (+ 1 (/ s (norm-plane s akm))) 2))))
                (let ((s-psi (/ akm (* 2 (norm-plane s akm) c-psi))))
                  (list k m c-psi s-psi)))
              (let ((c-psi (sqrt (/ (- 1 (/ s (norm-plane s akm))) 2))))
                (let ((s-psi (/ (* -1.0 akm) (* 2 (norm-plane s akm) c-psi))))
                  (list k m c-psi s-psi)))))))))
    
  ;;入力: '(k m cos(psi) sin(psi)), 行列A
  ;;出力: ギブンスの回転行列
  (define (make-P-from-Pdat-A Pdat A)
    (let ((k (car Pdat))
          (m (car (cdr*n Pdat 1)))
          (c-psi (car (cdr*n Pdat 2)))
          (s-psi (car (cdr*n Pdat 3)))
          (n (car (size-mat A))))
      (let ((m1 (- k 1))
            (m2 (- m k 1))
            (m3 (- n m)))
        (let ((els (lambda (lines counter-ini)
                     (do-n-times-with-filter-counter (lambda (ini filter counter)
                                                       (append ini (list (filter counter))))
                                                     '()
                                                     lines
                                                     (lambda (counter) (eye-line n counter))
                                                     (lambda (counter) (+ counter 1))
                                                     counter-ini))))
          (let ((em1 (els m1 1))
                (em2 (els m2 (+ k 1)))
                (em3 (els m3 (+ m 1))))
            (let ((linek (map +
                              (map (lambda (x) (* c-psi x))
                                   (eye-line n k))
                              (map (lambda (x) (* -1 s-psi x))
                                   (eye-line n m))))
                  (linem (map +
                              (map (lambda (x) (* s-psi x))
                                   (eye-line n k))
                              (map (lambda (x) (* c-psi x))
                                   (eye-line n m)))))
              (append em1 (list linek) em2 (list linem) em3)))))))
  
  (define (jacobi-method-p A tolerance tmpA tmpP)
    (if (< (abs (car (nondiag-biggest-abs tmpA))) tolerance)
      (cons tmpA tmpP)
;      (disp-mat tmpA) ;;デバッグ用
      (jacobi-method-p A
                       tolerance
                       (mul-mat (mul-mat (transposed-mat (make-P-from-Pdat-A (Pdat-from-mat tmpA) tmpA))
                                         tmpA)
                                (make-P-from-Pdat-A (Pdat-from-mat tmpA) tmpA))
                       (mul-mat tmpP
                                (make-P-from-Pdat-A (Pdat-from-mat tmpA) tmpA)))));) ;;displayデバック時

  ;;ここから本体 ヤコビ法
  (if (equal? (transposed-mat A) A)
    (jacobi-method-p A tolerance A (eye-mat (car (size-mat A))))
    (error "jacobi-method: (jacobi-method A tolerance) A should be a symetric matrix!!")))


;べき乗法
(define (power-iter A ini-v counter-limit)
  (define (power-iter-p A v1 v2 counter counter-limit)
    (if (>= counter counter-limit)
      (let ((l1 (car (transposed-mat v1)))
            (l2 (car (transposed-mat v2))))
        (let ((max-eiglist (map (lambda (x1 x2)
                                  (/ x2 x1))
                                l1
                                l2)))
          (begin
            (display "\n絶対値最大の固有値の近似値\n")
            (display max-eiglist)
            (display "\n\nそれに対応する固有ベクトル\n")
            (disp-mat (transposed-mat (list
                                       (map (lambda (x y)
                                              (/ x (expt y counter)))
                                            l2
                                            max-eiglist)))))))
      (power-iter-p A v2 (mul-mat A v2) (+ counter 1) counter-limit)))
  (power-iter-p A '() ini-v 0 counter-limit))


(define (first-nonzero list-a)
  (let ((u (find-1stitem-from-list
            (lambda (x)
              ;(not (= x 0)))
              (> (abs x) 0.00000000000001))
            list-a)))
    (if (not (equal? #f u))
      u
      (cons (+ (length list-a) 1) 0)))) ;;非ゼロが見つからなかったとき、位置をlength+1とする

;;ガウス消去法における前進消去
;n=列数
(define (gauss-elim-forward tmp-mat-up tmp-mat-low b-line zeros n option)
  (if (null? tmp-mat-low)
    (begin 
      (disp-mat (append tmp-mat-up zeros))
      (cond
       ;;出力 (rank . 行操作された行列)
       ((equal? option "solv") (cons (length tmp-mat-up) (append tmp-mat-up)))
       ((equal? option "rank") (cons (length tmp-mat-up) (append tmp-mat-up zeros)))
       (else (error "unknown option"))))
    ;ここからがelse
    (if (not (equal? (cdr (first-nonzero (car tmp-mat-low)))
                     0))
      (begin
        (set! tmp-mat-low (sort tmp-mat-low (lambda (x y)
                                              (< (car (first-nonzero x))
                                                 (car (first-nonzero y))))))
        (set! b-line (map (lambda (a) (/ a (cdr (first-nonzero (car tmp-mat-low)))))
                          (car tmp-mat-low))) ;; normalize
        (cond
         ((equal? option "solv")
          (if (= (car (first-nonzero b-line)) (length b-line))
            (error "GAUSS-ELIM-FORWARD: No Solution!")))
         ((equal? option "rank") (newline))
         (else (error "unknown option")))

        (disp-mat (append tmp-mat-up tmp-mat-low zeros))    ;デバッグ用
        (set! tmp-mat-up (append tmp-mat-up (list b-line)))
        (if (not (null? (cdr tmp-mat-low)))
          (set! tmp-mat-low (sub-mat (cdr tmp-mat-low)
                                     (x-th-scaled-list-to-each (car (first-nonzero (car tmp-mat-low)))
                                                               b-line (cdr tmp-mat-low))))
          (set! tmp-mat-low (cdr tmp-mat-low)))
        (gauss-elim-forward tmp-mat-up tmp-mat-low b-line zeros n option))	
      (gauss-elim-forward tmp-mat-up (cdr tmp-mat-low) b-line (append zeros (zero-mat 1 n)) n option))))

;;行列のランク
(define (rank A)
  (car (gauss-elim-forward '() A '() '() (cdr (size-mat A)) "rank")))

;;連立1次方程式Ax=bをガウスの消去法で解く
(define (div-gauss-mat A b)
  (define (gauss-elim-back tmp-mat-up tmp-mat-low)
    (if (null? tmp-mat-low)
      (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)))
          (set! tmp-mat-low (sub-mat (cdr tmp-mat-low)
                                     (x-th-scaled-list-to-each (car (first-nonzero (car tmp-mat-low)))
                                                               (car tmp-mat-low)
                                                               (cdr tmp-mat-low))))
          (set! tmp-mat-low (cdr tmp-mat-low)))
        (gauss-elim-back tmp-mat-up tmp-mat-low))))
  (let ((Ab (connect-mat-side A b)))
    (gauss-elim-back '() (reverse (cdr (gauss-elim-forward '() Ab '() '() (cdr (size-mat Ab)) "solv"))))))

;;逆行列
;;ガウスの消去法で求める
(define (inv A)
  (define (inv-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)))))
          (hakidashi-ahead (+ i 1) tmp-mat-up tmp-mat-low b-line size-mat-AE))))
    (define (hakidashi-back i tmp-mat-up tmp-mat-low size-mat-AE)                                                                                           
      (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)))

  (let ((N (car (size-mat A))))
    (extract-mat (inv-hakidashi (connect-mat-side A (eye-mat N)))
                 1 N (+ N 1) (* 2 N))))

和add-mat , スカラー倍scalar-mat , 引き算sub-mat , 積mul-mat , 転置行列transposed-mat
単位行列eye-mat

この部分の使い方に関しては、ほとんど変わっていない。これが一番最初の記事。作成日記 Schemeで数値計算(1) - nibosiiwasi’s blog

ただ、転置行列に関しては下記の変更を行った。
転置行列のコード - nibosiiwasi’s blog

その後、下記の数値計算法のコードを追加した。

ヤコビ法jacobi-method(実対称行列の固有値計算)

http://nibosiiwasi.hatenablog.com/entry/2017/03/22/230510

ガウスの消去法

(rank, 連立1次方程式の求解div-gauss-mat , 逆行列 inv)

ガウスの消去法における、前進消去の際、
first-nonzeroでxがゼロでないことを
(> (abs x) 0.00000000000001)
とした。
浮動小数点の演算の際、厳密にxが0にならない場合があるためである。

逆行列のコードについては、以前のブログ
http://nibosiiwasi.hatenablog.com/entry/2016/10/02/230615
と基本的に同じである。

今回、行列Aの逆行列
(inv A)
と書くことにした。



以下、rankと連立1次方程式の求解の例を挙げる。

行列のランク

(define A
  '((1 2 3 4)
    (2 1 3 0)
    (-2 3 1 1)))

(rank A)

実行結果
前進消去の計算過程も表示。

(1 2 3 4)
(2 1 3 0)
(-2 3 1 1)



(1 2 3 4)
(0 -3 -3 -8)
(0 7 7 9)



(1 2 3 4)
(0 1 1 8/3)
(0 0 0 -29/3)


(1 2 3 4)
(0 1 1 8/3)
(0 0 0 1)

3
連立1次方程式の求解

線形方程式Ax=bを解く。ガウスの消去法を使い、行列の形で答えとした。

(define A
  '((1 -1 0)
    (3 -1 1)
    (2 -1 2)
    (0 1 -1)))

(define b (transposed-mat '((-2 -2 -1 1))))

(disp-mat (div-gauss-mat A b))

実行結果
前進消去の計算過程も表示。(後退代入は表示させていない。)

(1 -1 0 -2)
(3 -1 1 -2)
(2 -1 2 -1)
(0 1 -1 1)


(1 -1 0 -2)
(0 2 1 4)
(0 1 2 3)
(0 1 -1 1)


(1 -1 0 -2)
(0 1 1/2 2)
(0 0 3/2 1)
(0 0 -3/2 -1)


(1 -1 0 -2)
(0 1 1/2 2)
(0 0 1 2/3)
(0 0 0 0)


(1 0 0 -1/3)
(0 1 0 5/3)
(0 0 1 2/3)

#<undef>


本ページのコードのアルゴリズムは全て、下記の本を参考にしました。
・洲之内著,石渡改訂:数値計算,サイエンス社

コードはご自由にお使い下さい。
コードの使用により受けた損害について、著者は一切の責任を負いません。