作成日記 Schemeで数値計算(1)

線形代数周辺の問題を解けるようにすることは、私の大きな目標の１つである。数値計算をする上で、線形代数、行列計算は避けて通れない。

まず、行列計算の本当に基本的な部分をlist操作を使って書く。Gaucheにはarrayもあるのは知っているが、あえて使わない方針をとる。

Schemeらしく、高階関数やlambda,mapを使い、短く書くことに努める。
mapとほぼ同じであるが、my-map-item2を自作した。これは引数の２つの行列のサイズが異なるときに、エラーを吐かせるためである。

行列の転置 transposed-mat、和 add-mat、スカラー倍 scalar-mat、積 mul-mat、減 sub-mat、単位行列eye-mat を定義する。
col1sep-matは行列の1列目とそれ以外を分離する関数である。col2lineは1列を1行に変換する。col1sep-matとcol2lineを使って行列の転置を定義した。

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

(define (add-mat A1 A2)
  (define (add-each a1 a2) (my-map-item2 + a1 a2))
  (my-map-item2 add-each A1 A2))

(define (scalar-mat a A)
  (define (scalar-each a a1) (map (lambda (x) (* a x)) a1))
  (map (lambda (x) (scalar-each a x)) A))

(define (sub-mat A1 A2) (add-mat A1 (scalar-mat -1.0 A2)))

(define (l-accum2 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))))

(define (sum-of-mul-each a1 a2) (l-accum2 + * 0 a1 a2))

(define (col1sep-mat A)
  (define (col1sep-mat-pre  A tmp1 tmp2)
    (if (null? A)
	(cons tmp1 tmp2)
	(col1sep-mat-pre (cdr A)
				(append tmp1 (list (list (caar A))))
				(append tmp2 (list (cdr (car A)))))))
  (col1sep-mat-pre A '() '() ))

(define (col2line col)
  (define (col2line-p col tmp)
    (if (null? col)
	tmp
	(col2line-p (cdr col)
		    (append tmp (list (caar col))))))
  (col2line-p col '() ))

(define (transposed-mat A)
  (define (transposed-mat-p A tmp)
    (if (null? (car A))
	tmp
	(transposed-mat-p (cdr (col1sep-mat A))
			  (append tmp
				  (list (col2line (car (col1sep-mat A))))))))
  (transposed-mat-p A '() ))

(define (mul-oneline-mat oneline A2)
  (map (lambda (x) (sum-of-mul-each oneline x)) (transposed-mat A2)))

(define (mul-mat A1 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 (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 0))))))
  (eye-line-p len loc1 '() ))

(define (eye-mat N)
  (define (eye-mat-pre N i tmp)
    (if (> i N)
	tmp
	(eye-mat-pre N (+ i 1) (append tmp (list (eye-line N i))))))
  (eye-mat-pre N 1 '()))

動作確認

(define A1 '( (1 2 3)
	      (4 5 6)))

;;A1
;;=> ((1 2 3) (4 5 6))


(define tA1 (transposed-mat A1))

;;tA1
;;=> ((1 4) (2 5) (3 6))


(define eye3 (eye-mat 3))

;;eye3
;;=> ((1.0 0 0) (0 1.0 0) (0 0 1.0))


;;(add-mat A1 A1)
;;=> ((2 4 6) (8 10 12))

;;(scalar-mat 3 A1)
;;=> ((3 6 9) (12 15 18))


;;(mul-mat A1 (eye-mat 3))
;;=>((1.0 2.0 3.0) (4.0 5.0 6.0))

;;(mul-mat A1 tA1)
;;=>((14 32) (32 77))

つづく