여러 가지 도구를 이용한 행렬 계산 / 행렬 곱셈, 전치행렬, 행렬식, 역행렬 등

        A = [1, 2; 3, 4]
        B = [2, 1; -1, 2]

로 주어진 두 2x2 행렬 A, B 에 대하여

        AB, BA. transpose(A), transpose(A) B
        det(A), inverse(A),
        inverseA) A, A inverse(A)

를 각각 구하는 방법을 여러 가지 도구

        Mathematica, Maxima, Octave, Scipy

들을 각각 이용하여 알아본다.
,


* Mathematica 를 이용하여 행렬 계산하기 (행렬 곱셈, 전치행렬, 행렬식, 역행렬 등)
  (참고: Mathematica 에서 행렬 곱샘 연산자는 .(점) 이다.)

 

* Maxima 를 이용하여 행렬 계산하기 (행렬 곱셈, 전치행렬, 행렬식, 역행렬 등)
  (참고: Maxima 애서도 Mathematica 에서와 마찬가지로 행렬 곱샘 연산자는 .(점) 이다.)
** 행렬의 생성과 곱셈

 

** 전치행렬, 행렬식, 역행렬

 

* Octave 를 이용하여 행렬 계산하기 (스칼라배, 합, 차 등)
** Octave 소스

A = [1, 2; 3, 4]
B = [2, 1; -1, 2]
A*B, B*A
A', A'*B
transpose(A), transpose(A)*B
det(A), inverse(A), A^(-1), A^-1
A^-1 * A, A * A^-1
inverse(A) * A, A * inverse(A)

** Octave 명령 창에서 행렬 A, B 생성하고 행렬 곱셈 계산하기

 

** Octave 명령 창에서 전치행렬, 행렬식, 역행렬 구하기

** Octave 명령 창에서 행렬 곱셈으로 역행렬 검증하기

 

* Python 2.7 & Scipy를 이용하여 행렬 계산하기 (행렬 곱셈, 전치행렬, 행렬식, 역행렬 등)
** 파이썬 소스

import numpy as np
import scipy as sp
from scipy import mat
from scipy import linalg

A = mat( "[1, 2; 3, 4]" )
B = mat( "[2, 1; -1, 2]" )
A; B
A*B; B*A
np.transpose(A); np.transpose(A)*B
linalg.det(A); linalg.inv(A); A.I
A*(A.I); (A.I)*A

** 위의 소스를 파일로 저장하지 말고, Python 인터프리터 실행 창에서 한 줄씩 입력한다.

 

* Python 2.6 & 2.7 에서 동작하는 분수만으로 이루어진 정방행렬의
   행렬식, 역행렬 구하는 함수를 구현해 보았다.
   (참고: scipy.linalg 모듈은 분수를 부동소수점수로 바꾸어 계산하므로
            역행렬에서  분수 성질이 없어진다.)

"""
  Filename: TestFraction.py

            Test fraction calculations, wgich works on Python 2.6 or 2.7.

   Execute: python TestFraction.py
      Date: 2011. 10. 10
"""

from fractions import *
import numpy as np
import scipy as sp
from scipy import *

def evalFraction(x):
 return x.numerator/float(x.denominator)

def printMatrix(msg, m):
 if msg != None and len(msg) > 0:
  print msg
 print "%s" % m

def minorMatrix(m, x, y):
 nr, nc = m.shape
 t = []
 for i in range(0, nr - 1):
  tmp = []
  for j in range(0, nc - 1):
   if i < x:
    if j < y:
     u = m[i, j]
     tmp.append(u)
    else:
     u = m[i, j+1]
     tmp.append(u)
   else:
    if j < y:
     u = m[i+1, j]
     tmp.append(u)
    else:
     u = Fraction(m[i+1, j+1])
     tmp.append(u)
  t.append(tmp)
 t2 = matrix(t)
 return t2
 
def det(m):
 nr, nc = m.shape
 if nr != nc:
      print "Caught TypeError"
      return None
 if nr == 1:
  return m[0, 0]
 elif nr == 2:
          return m[0, 0]*m[1, 1] - m[0, 1]*m[1, 0]
 else:
  sum = 0
  sign = 1
  for j in range(0, nc):
   minor = minorMatrix(m, 0, j)
   sum = sum + sign*m[0, j]*det(minor)
   sign = -sign
  return sum
 return None

def inverse(m):
 nr, nc = m.shape
 if nr != nc:
      print "Caught TypeError"
      return None
 if nr == 1:
  return matrix([[1/m[0, 0]]])
 elif nr == 2:
          d = [1, 1] - m[0, 1]*m[1, 0]
          return matrix( [[m[1, 1]/d, m[0,1]/d], [m[1, 0]/d, m[1, 1]/d]] )
 else:
  d = det(m)
  sign = 1
  t = []
  for j in range(0, nc):
   tmp = []
   for i in range(0, nr):
       minor = minorMatrix(m, i, j)
       sign = (-1)**(i + j)
       v = sign*det(minor)
       tmp.append(v/d)
   t.append(tmp)
  t2 = matrix(t)
  return t2
 return None

def hilbertMatrix(n):
 if n <= 0:
  return None
 t = []
 for i in range(0, n):
  tmp = []
  for j in range(0, n):
      tmp.append(Fraction(1, i + j + 1))
  t.append(tmp)
 t2 = matrix(t)
 return t2

z = Fraction(1, 2)
w = Fraction("5/4")

A = mat([[Fraction(1,5), Fraction(3,5)], [Fraction(3,5), Fraction(4,5)]])
B = mat([[Fraction(4,3), Fraction(-1,3)], [Fraction(2,3), Fraction(7,3)]])
printMatrix("A = ", A)
printMatrix("B = ", B)
printMatrix("A + B = ", A + B)
printMatrix("A - B = ", A - B)
printMatrix("A * B = ", A * B)
printMatrix("A / B = ", A / B)
printMatrix("A**2 = ", A**2)
printMatrix("A**3 = ", A**3)
printMatrix("B**-1 = ", B**-1)
printMatrix("B*B**-1 = ", B*B**-1)
printMatrix("B.transpose() = ", B.transpose())

print "B[0, 0] = %s" % (B[0, 0])
print "B[0, 1] = %s" % (B[0, 1])
print "B[1, 0] = %s" % (B[1, 0])
print "B[1, 1] = %s" % (B[1, 1])

B[0][0] = Fraction(5, 1) 
B[0][0] = 50
printMatrix("Changed: B = ", B)
print "  ---> det(B) = %s" % det(B)
B[0,0] = 80
printMatrix("Changed again: B = ", B)
print "  ---> det(B) = %s" % det(B)

C = mat([ [Fraction(1,1), Fraction(1,2), Fraction(1,3)],
          [Fraction(1,2), Fraction(1,3), Fraction(1,4)], 
          [Fraction(1,3), Fraction(1,4), Fraction(1,5)] ])
printMatrix("C = ", C)
v = det(C)
print "  ---> det(C) = %s" % v
print "       evalFraction(det(C)) = %s" % evalFraction(v)

from scipy import linalg
print "Cehck again by using scipy.linalg(): "
print "       linalg.det(C) = %s" % linalg.det(C)

D = inverse(C)
print "inverse(C) = \n%s" % D
printMatrix("inverse(C)*C = ", D*C)
printMatrix("C*inverse(C) = ", C*D)
printMatrix("linalg.inv(C) = ", linalg.inv(C))
printMatrix("linalg.inv(C)*C = ", linalg.inv(C)*C)

H = hilbertMatrix(4)
printMatrix("H = ", H)
v = det(H)
print "  ---> det(H) = %s" % v
print "       evalFraction(det(H)) = %s" % evalFraction(v)
print "Cehck again by using scipy.linalg(): "
print "       linalg.det(H) = %s" % linalg.det(H)
G = inverse(H)
print "inverse(H) = \n%s" % G
printMatrix("inverse(H)*H = ", G*H)
printMatrix("H*inverse(H) = ", H*G)
printMatrix("linalg.inv(H) = ", linalg.inv(H))
printMatrix("linalg.inv(H)*H = ", linalg.inv(H)*H)

"""
Execute Result:
A =
[[1/5 3/5]
 [3/5 4/5]]
B =
[[4/3 -1/3]
 [2/3 7/3]]
A + B =
[[23/15 4/15]
 [19/15 47/15]]
A - B =
[[-17/15 14/15]
 [-1/15 -23/15]]
A * B =
[[2/3 4/3]
 [4/3 5/3]]
A / B =
[[3/20 -9/5]
 [9/10 12/35]]
A**2 =
[[2/5 3/5]
 [3/5 1]]
A**3 =
[[11/25 18/25]
 [18/25 29/25]]
B**-1 =
[[ 0.7  0.1]
 [-0.2  0.4]]
B*B**-1 =
[[1.0 0.0]
 [-5.55111512313e-17 1.0]]
B.transpose() =
[[4/3 2/3]
 [-1/3 7/3]]
B[0, 0] = 4/3
B[0, 1] = -1/3
B[1, 0] = 2/3
B[1, 1] = 7/3
Changed: B =
[[50 50]
 [2/3 7/3]]
  ---> det(B) = 250/3
Changed again: B =
[[80 50]
 [2/3 7/3]]
  ---> det(B) = 460/3
C =
[[1 1/2 1/3]
 [1/2 1/3 1/4]
 [1/3 1/4 1/5]]
  ---> det(C) = 1/2160
       evalFraction(det(C)) = 0.000462962962963
Cehck again by using scipy.linalg():
       linalg.det(C) = 0.000462962962963
inverse(C) =
[[9 -36 30]
 [-36 192 -180]
 [30 -180 180]]
inverse(C)*C =
[[1 0 0]
 [0 1 0]
 [0 0 1]]
C*inverse(C) =
[[1 0 0]
 [0 1 0]
 [0 0 1]]
linalg.inv(C) =
[[   9.  -36.   30.]
 [ -36.  192. -180.]
 [  30. -180.  180.]]
linalg.inv(C)*C =
[[1.0 0.0 0.0]
 [0.0 1.0 0.0]
 [0.0 0.0 1.0]]
H =
[[1 1/2 1/3 1/4]
 [1/2 1/3 1/4 1/5]
 [1/3 1/4 1/5 1/6]
 [1/4 1/5 1/6 1/7]]
  ---> det(H) = 1/6048000
       evalFraction(det(H)) = 1.65343915344e-07
Cehck again by using scipy.linalg():
       linalg.det(H) = 1.65343915344e-07
inverse(H) =
[[16 -120 240 -140]
 [-120 1200 -2700 1680]
 [240 -2700 6480 -4200]
 [-140 1680 -4200 2800]]
inverse(H)*H =
[[1 0 0 0]
 [0 1 0 0]
 [0 0 1 0]
 [0 0 0 1]]
H*inverse(H) =
[[1 0 0 0]
 [0 1 0 0]
 [0 0 1 0]
 [0 0 0 1]]
linalg.inv(H) =
[[   16.  -120.   240.  -140.]
 [ -120.  1200. -2700.  1680.]
 [  240. -2700.  6480. -4200.]
 [ -140.  1680. -4200.  2800.]]
linalg.inv(H)*H =
[[1.0 0.0 0.0 -3.5527136788e-15]
 [0.0 1.0 -5.68434188608e-14 2.84217094304e-14]
 [0.0 0.0 1.0 1.13686837722e-13]
 [0.0 0.0 5.68434188608e-14 1.0]]
"""