逆行列を求める

逆行列を計算するプログラムです。

まず共通処理から。

(define nil '())

(define *stack* nil)
(define level-ratio 100)

(define (last-pushed)
  (if (empty?)
      nil
      (car *stack*)))

(define (after-popped)
  (if (empty?)
      nil
      (cdr *stack*)))

(define (init)
  (set! *stack* nil))

(define (empty?)
  (null? *stack*))

(define (push x)
  (set! *stack* (cons x *stack*)))

(define (pop)
  (if (empty?)
      nil
      (let ((result (last-pushed)))
        (set! *stack* (after-popped))
        result)))

(define (pop-all level)
  (define (iter operators)
    (if (or (empty?) (< (priority (last-pushed)) (* level level-ratio)))
        operators
        (iter (append operators (cons (operator (pop)) nil)))))
  (iter nil))

;; r.p.n. procedures

(define (make-priority op level)
  (let ((priority
         (cond ((eq? op '**) 30)
               ((eq? op '*)  20)
               ((eq? op '/)  20)
               ((eq? op '+)  10)
               ((eq? op '-)  10)
               (else 0))))
    (cons op (+ priority (* level level-ratio)))))

(define (operator p)
  (car p))

(define (priority p)
  (if (null? p)
      0
      (cdr p)))

(define (operator? x)
  (or (eq? x '+) (eq? x '-) (eq? x '*) (eq? x '/) (eq? x '**)))

(define (exp->rpn exp)
  (define (iter e r level)
    (if (null? e)
        (append r (pop-all level))
        (if (list? (car e))
            (iter (cdr e) (iter (car e) r (+ level 1)) level)
            (cond ((operator? (car e))
                   (cond ((> (priority (last-pushed)) (priority (make-priority (car e) level)))
                          (let ((p (list (operator (pop)))))
                            (push (make-priority (car e) level))
                            (iter (cdr e) (append r p) level)))
                         (else
                          (push (make-priority (car e) level))
                          (iter (cdr e) r level))))
                  (else (iter (cdr e) (append r (list (car e))) level))))))
  (init)
  (iter exp nil 0))

(define (rpn->s-exp rpn)
  (define (iter r)
      (if (null? r)
          (pop)
          (let ((e (car r)))
            (cond ((eq? e '+) (let ((augend (pop)))
                                (push (make-sum (pop) augend))))
		  ((eq? e '-) (let ((subtrahend (pop)))
				(push (make-sub (pop) subtrahend))))
                  ((eq? e '*) (let ((multiplicand (pop)))
                                (push (make-product (pop) multiplicand))))
		  ((eq? e '/) (let ((divisor (pop)))
				(if (and (number? divisor) (= divisor 0))
				    (error "attempt to calculate a division by zero")
				    (push (make-division (pop) divisor)))))
                  ((eq? e '**) (let ((exponent (pop)))
                                 (push (make-exponentiation (pop) exponent))))
                  (else (push e)))
            (iter (cdr r)))))
  (init)
  (iter rpn))

次は、記号微分の処理。今回は微分はしないですし、いらない処理がたくさんあります。

(define (deriv exp var)
  (define (deriv-iter exp var)
    (cond ((number? exp) 0)
          ((variable? exp)
           (if (same-variable? exp var) 1 0))
          ((sum? exp)
           (make-sum (deriv-iter (addend exp) var)
                     (deriv-iter (augend exp) var)))
          ((product? exp)
           (make-sum
            (make-product (multiplier exp)
                          (deriv-iter (multiplicand exp) var))
            (make-product (deriv-iter (multiplier exp) var)
                          (multiplicand exp))))
          ((exponentiation? exp)
           (make-product
            (make-product (exponent exp)
                          (make-exponentiation (base exp) (make-sum (exponent exp) -1)))
            (deriv-iter (base exp) var)))
          (else
           (error "unknown expression type -- DERIV" exp))))
  (deriv-iter (rpn->s-exp (exp->rpn exp)) var))

(define (variable? x) (symbol? x))

(define (same-variable? v1 v2)
  (and (variable? v1) (variable? v2) (eq? v1 v2)))

(define (=number? exp num)
  (and (number? exp) (= exp num)))

(define (make-sum a1 a2)
  (cond ((=number? a1 0) a2)
        ((=number? a2 0) a1)
        ((and (number? a1) (number? a2)) (+ a1 a2))
        (else (list '+ a1 a2))))

(define (make-sub a1 a2)
  (cond ((=number? a1 0) (* -1 a2))
        ((=number? a2 0) a1)
        ((and (number? a1) (number? a2)) (- a1 a2))
        (else (list '- a1 a2))))

(define (make-product m1 m2)
  (cond ((or (=number? m1 0) (=number? m2 0)) 0)
        ((=number? m1 1) m2)
        ((=number? m2 1) m1)
        ((and (number? m1) (number? m2)) (* m1 m2))
        (else (list '* m1 m2))))

(define (make-division m1 m2)
  (cond ((=number? m1 0) 0)
	((=number? m2 0) (error "attempt to calculate a division by zero"))
        ((=number? m2 1) m1)
        ((and (number? m1) (number? m2)) (/ m1 m2))
        (else (list '/ m1 m2))))

(define (make-exponentiation e1 e2)
  (cond ((=number? e2 0) 1)
        ((=number? e2 1) e1)
        (else (list '** e1 e2))))

(define (sum? x)
  (and (pair? x) (eq? (car x) '+)))

(define (addend s) (cadr s))

(define (augend s) (caddr s))

(define (product? x)
  (and (pair? x) (eq? (car x) '*)))

(define (multiplier p) (cadr p))

(define (multiplicand p) (caddr p))

(define (exponentiation? x)
  (and (pair? x) (eq? (car x) '**)))

(define (base e) (cadr e))

(define (exponent e) (caddr e))

最後に逆行列を計算するプログラム。

(define (make-key r c)
  `(,r ,c))

(define (square-matrix? matrix)
  (let ((row (length matrix)))
    (equal? (make-list row row) (map length matrix))))
  
(define (make-unit-matrix matrix)
  (define (iter i lst mat)
    (if (= i 0)
	mat
	(iter (- i 1) (append (cdr lst) (cons 0 nil)) (cons lst mat))))
  (let ((len (length matrix)))
    (iter len (append (make-list (- len 1) 0) '(1)) nil)))

(define (matrix->hash-table matrix)
  (define ht (make-hash-table 'equal?))
  (define (col-iter r c cols)
    (if (null? cols)
	'done
	(begin
	  (hash-table-set! ht (make-key r c) (car cols))
	  (col-iter r (+ c 1) (cdr cols)))))
  (define (row-iter r rows)
    (if (null? rows)
	ht
	(begin
	  (col-iter r 1 (car rows))
	  (row-iter (+ r 1) (cdr rows)))))
  (row-iter 1 matrix))

(define (hash-table->row ht size row)
  (define (iter c row-list)
    (if (> c size)
        row-list
        (iter (+ c 1) (append row-list (cons (hash-table-get ht (make-key row c)) nil)))))
  (iter 1 nil))

(define (hash-table->matrix ht size)
  (define (r-iter r matrix)
    (if (> r size)
        matrix
        (r-iter (+ r 1) (append matrix (cons (hash-table->row ht size r) nil)))))
  (r-iter 1 nil))

(define (solve-main ht1 ht2 size)
  (define (change-row row1 row2)
    (define (iter c r1 r2 r3 r4)
      (if (> c size)
          'done
          (begin
            (hash-table-set! ht1 (make-key row1 c) (car r2))
            (hash-table-set! ht1 (make-key row2 c) (car r1))
            (hash-table-set! ht2 (make-key row1 c) (car r4))
            (hash-table-set! ht2 (make-key row2 c) (car r3))
            (iter (+ c 1) (cdr r1) (cdr r2) (cdr r3) (cdr r4)))))
    (let ((r1 (hash-table->row ht1 size row1))
          (r2 (hash-table->row ht1 size row2))
          (r3 (hash-table->row ht2 size row1))
          (r4 (hash-table->row ht2 size row2)))
      (iter 1 r1 r2 r3 r4)))
    
  (define (irekae pivot-col)
    (define (r-iter r max-r max-val)
      (if (> r size)
          (change-row pivot-col max-r)
          (let ((val (hash-table-get ht1 (make-key r pivot-col))))
            (if (and (number? val) (> (abs val) max-val))
                (r-iter (+ r 1) r val)
                (r-iter (+ r 1) max-r max-val)))))
    (r-iter pivot-col pivot-col (hash-table-get ht1 (make-key pivot-col pivot-col))))
  
  (define (taikaku pivot-col)
    (define (iter val c)
      (if (> c size)
	  'done
          (begin
	    (hash-table-set! ht1 (make-key pivot-col c) (make-division (hash-table-get ht1 (make-key pivot-col c)) val))
	    (hash-table-set! ht2 (make-key pivot-col c) (make-division (hash-table-get ht2 (make-key pivot-col c)) val))
	    (iter val (+ c 1)))))
    (let ((val (hash-table-get ht1 (make-key pivot-col pivot-col))))
      (and (not (= val 1)) (iter val 1))))

  (define (hitaikaku pivot-col)
    (define (c-iter val row col)
      (if (> col size)
          'done
          (begin
            (hash-table-set! ht1
                             (make-key row col)
                             (make-sum (hash-table-get ht1 (make-key row col))
                                       (make-product (hash-table-get ht1 (make-key pivot-col col))
                                                     val)))
            (hash-table-set! ht2
                             (make-key row col)
                             (make-sum (hash-table-get ht2 (make-key row col))
                                       (make-product (hash-table-get ht2 (make-key pivot-col col))
                                                     val)))
            (c-iter val row (+ col 1)))))
    
    (define (r-iter row)
      (if (> row size)
          'done
          (begin
            (let ((val (hash-table-get ht1 (make-key row pivot-col))))
              (and (not (= row pivot-col)) (c-iter (make-product -1 val) row 1)))
            (r-iter (+ row 1)))))
    (r-iter 1))
  
  (define (iter pivot-col)
    (if (> pivot-col size)
        ht2
        (begin
          (irekae pivot-col)
          (taikaku pivot-col)
          (hitaikaku pivot-col)
          (iter (+ pivot-col 1)))))
  (iter 1)
  )

(define (solve matrix)
  (let ((size (length matrix)))
    (hash-table->matrix (solve-main (matrix->hash-table matrix)
                                    (matrix->hash-table (make-unit-matrix matrix))
                                    size)
                        size)))

(define (inverse-matrix matrix)
  (if (not (square-matrix? matrix))
      (error "Require square matrix.")
      (solve matrix)))

掃き出し法です。部分ピボット機能付き。
早速、逆行列を求めてみます。

gosh> (inverse-matrix '((1 3 2) (-1 0 1) (2 3 0)))
((1 -2 -1) (-2/3 4/3 1) (1 -1 -1))

検算してみます。
行列の逆行列の逆行列は、元の行列になります。

gosh> (inverse-matrix (inverse-matrix '((1 3 2) (-1 0 1) (2 3 0))))
((1 3 2) (-1 0 1) (2 3 0))

はい。元に戻りました。問題なさそうです。

gosh> (inverse-matrix '((1 2 7 6) (2 4 4 2) (1 8 5 2) (2 4 3 3)))
((-1/6 7/12 -1/3 1/6) (-1/15 -13/60 1/6 1/6) (1/10 9/20 0 -1/2) (1/10 -11/20 0 1/2))

n x n の行列も計算できます。
Gauche は有理数をサポートするので結果が見やすく誤差が生じず素晴らしいです。

なお、デバッグで以下のサイトを利用させていただきました。
ありがとうございます。
keisan.casio.jp

応用編
数字だけじゃつまんないですよね?
というかこっちがメインです。

例えば、技術士試験 平成 30 年度の基礎科目 I-3-3 に出題されたような行列の逆行列を求めてみます。
技術士試験の過去問題は、日本技術士会のホームページで公開されています。
www.engineer.or.jp

gosh> (inverse-matrix '((1 0 0) (a 1 0) (b c 1)))
((1 0 0) ((* -1 a) 1 0) ((+ (* -1 b) (* (* -1 a) (* -1 c))) (* -1 c) 1))

式の最適化に少し難があり課題が残りますが解けました。
結果は常に前置記法です。プログラムが LISP ですから。

なお、中間記法は以下のように前置記法に変換できます。
(注 LISP ではべき乗は「**」ではないので、この式は変数 x を定義しても LISP 処理系では実行できません。)

gosh> (rpn->s-exp (exp->rpn '(x ** 2 + 4 * x + 3)))
(+ (** x 2) (+ (* 4 x) 3))

しかし、こんなプログラムを書いていたのでは、試験の勉強が捗りませんな!
まったく困ったものです。

以下蛇足
式の微分も出来ます。

gosh> (deriv '(x ** 2 + 4 * x + 3) 'x)
(+ (* 2 x) 4)

詳しくは SICP の 2.3.2 節「例：記号微分」を参照ください。
https://sicp.iijlab.net/fulltext/xcont#s232