ぱたへね

はてなダイアリーはrustの色分けができないのでこっちに来た

Gauche-CVで行列の計算を行う

Gauche-CVで簡単な行列計算をする方法です。当たり前ですがnumpyと同じ結果になります。

Python numpy で行列計算

Pythonではnumpyを使えば、簡単に行列の計算ができます。

# -*- coding: utf-8 -*-
__author__ = 'Natsutani'

import numpy as np

a = np.array([[ 1.0, -2.0, 3.0], [4.0, 5.0, 6.0], [-7.0, 8.0, -9.0]])
b = np.array([[10.0, -9.0, 8.0], [7.0, 6.0, 4.0], [-3.0, 2.0, -1.0]])

print("A + B")
print(a + b)

print("A - B")
print(a - b)

print("A * B (要素同士の積)")
print(a * b)

print("Aの逆行列")
print(np.linalg.inv(a))

print("A ・ B")
print(np.dot(a, b))

print("det A")
print(np.linalg.det(a))

print("転置 A")
print(np.transpose(a))

print("A + 5")
print(a + 5)

print("Aの5倍")
print(a * 5)

# 3次元ベクトルの内積・外積
c = np.array([ 1.0, -2.0, 3.0],)
d = np.array([ 4.0, -5.0, 6.0],)

print("C・D")
print(np.dot(c, d))

print("C×D")
print(np.cross(c, d))

実行結果

A + B
[[ 11. -11.  11.]
 [ 11.  11.  10.]
 [-10.  10. -10.]]
A - B
[[-9.  7. -5.]
 [-3. -1.  2.]
 [-4.  6. -8.]]
A * B (要素同士の積)
[[ 10.  18.  24.]
 [ 28.  30.  24.]
 [ 21.  16.   9.]]
Aの逆行列
[[-0.775       0.05       -0.225     ]
 [-0.05        0.1         0.05      ]
 [ 0.55833333  0.05        0.10833333]]
A ・ B
[[-13. -15.  -3.]
 [ 57.   6.  46.]
 [ 13.  93. -15.]]
det A
120.0
転置 A
[[ 1.  4. -7.]
 [-2.  5.  8.]
 [ 3.  6. -9.]]
A + 5
[[  6.   3.   8.]
 [  9.  10.  11.]
 [ -2.  13.  -4.]]
Aの5倍
[[  5. -10.  15.]
 [ 20.  25.  30.]
 [-35.  40. -45.]]
C・D
32.0
C×D
[ 3.  6.  3.]

Gauche-CVで行列計算

行列とスカラー値の和、積はOpenCVに見つからなかったので自作しました。パフォーマンスを気にしなければ簡単です。

(use cv)

;; 行列の表示
(define print-2d-matrix
  (lambda (m . args)
    (let-keywords args ((name ""))
                  (unless (string=? name "")
                    (format #t "~A = " name))
                  (let ((rows (slot-ref m 'rows))
                        (cols (slot-ref m 'cols)))
                    (display "[")
                    (dotimes (r rows)
                             (display " [")
                             (dotimes (c cols)
                                      (format #t " ~A"  (cv-get-real2d m r c)))
                             (display " ]"))
                    (display "]\n")))))

;; 行列とスカラー値の和  
(define cv-mat-adds-2d
  (lambda (src value dst)
    (cv-mat-operate-2d src dst (lambda (x) (+ x value)))))

;; 行列とスカラー値の和  
(define cv-mat-muls-2d
  (lambda (src value dst)
    (cv-mat-operate-2d src dst (lambda (x) (* x value)))))

;; 行列の各要素に演算を行う                                        
(define cv-mat-operate-2d
  (lambda (src dst op)
    (let ((rows (slot-ref src 'rows))
          (cols (slot-ref src 'cols)))
      (dotimes (r rows)
               (dotimes (c cols)
                        (cv-set-real2d dst r c
                                       (op (cv-get-real2d src r c))))))))

(define a (make-cv-mat-from-uvector 3 3 1 #f64(1.0 -2.0 3.0 4.0 5.0 6.0 -7.0 8.0 -9.0)))
(define b (make-cv-mat-from-uvector 3 3 1 #f64(10.0 -9.0 8.0 7.0 6.0 4.0 -3.0 2.0 -1.0)))
(define x (cv-clone-mat a))

(print-2d-matrix a :name "A")
(print-2d-matrix b :name "B")

(cv-add a b x)
(print-2d-matrix x :name "A + B")

(cv-sub a b x)
(print-2d-matrix x :name "A - B")

(cv-mul a b x)
(print-2d-matrix x :name " A mul B")

(cv-invert a x)
(print-2d-matrix x :name "inv(A)")

(cv-mat-mul a b x)
(print-2d-matrix x :name "A・B")

;(cv-cross-product a b x)
;(print-2d-matrix x :name "a×b")

(format #t "det A = ~A~%" (cv-det a))

(cv-transpose a x)
(print-2d-matrix x  :name "転置 A")

(cv-mat-adds-2d a 5 x)
(print-2d-matrix x :name "A+5")

(cv-mat-muls-2d a 5 x)
(print-2d-matrix x :name "A*5")

;; 3次元ベクトルの内積、外積
(define c (make-cv-mat-from-uvector 3 1 1 #f64(1.0 -2.0 3.0)))
(define d (make-cv-mat-from-uvector 3 1 1 #f64(4.0 -5.0 6.0)))
(define x2 (cv-clone-mat c))
 
(format #t "C・D = ~A~%" (cv-dot-product c d))
(cv-cross-product c d x2)
(print-2d-matrix x2 :name "A×B")

実行結果

A = [ [ 1.0 -2.0 3.0 ] [ 4.0 5.0 6.0 ] [ -7.0 8.0 -9.0 ]]
B = [ [ 10.0 -9.0 8.0 ] [ 7.0 6.0 4.0 ] [ -3.0 2.0 -1.0 ]]
A + B = [ [ 11.0 -11.0 11.0 ] [ 11.0 11.0 10.0 ] [ -10.0 10.0 -10.0 ]]
A - B = [ [ -9.0 7.0 -5.0 ] [ -3.0 -1.0 2.0 ] [ -4.0 6.0 -8.0 ]]
 A mul B = [ [ 10.0 18.0 24.0 ] [ 28.0 30.0 24.0 ] [ 21.0 16.0 9.0 ]]
inv(A) = [ [ -0.775 0.05 -0.225 ] [ -0.05 0.1 0.05 ] [ 0.5583333333333333 0.05 0.10833333333333334 ]]
A・B = [ [ -13.0 -15.0 -3.0 ] [ 57.0 6.0 46.0 ] [ 13.0 93.0 -15.0 ]]
det A = 120.0
転置 A = [ [ 1.0 4.0 -7.0 ] [ -2.0 5.0 8.0 ] [ 3.0 6.0 -9.0 ]]
A+5 = [ [ 6.0 3.0 8.0 ] [ 9.0 10.0 11.0 ] [ -2.0 13.0 -4.0 ]]
A*5 = [ [ 5.0 -10.0 15.0 ] [ 20.0 25.0 30.0 ] [ -35.0 40.0 -45.0 ]]
C・D = 32.0
A×B = [ [ 3.0 ] [ 6.0 ] [ 3.0 ]]