여러 가지 도구를 이용한 행렬 계산 / LU 분해(decomposition)

여기서는 2x2 정방행렬

        A = [4, 3; 6, 3]

의 LU 분해(LU decomposition)를 위해 여러 가지 도구

        Mathematica, Maxima, Octave, Scipy, 
        Jama 패키지를 이용한 Java 애플리케이션,
        Jama 패키지를 이용한 Groovy 애플리케이션

로는 각각 어떻게 해결하는지 알아보고자 한다.
,


* Mathematica 를 이용하여 행렬 계산하기
  (참고: Mathematica 에서는 LUDecompostion 이라는 함수가 준비되어 있다.)

* Maxima 를 이용하여 행렬 계산하기 (행렬 곱셈, 전치행렬, 행렬식, 역행렬 등)
  (참고: Maxima 애서는 변수에 수식 계산에 의한 값을 할당할 때 쓰는 할당 연샂자가 등호(=)가
           아니라 콜론(:)임에 유의한다.   Maxima 에서는 LU 분해를 위해 lu_factor 함수와 
           get_lufactors 함수를 연이어 사용한다,)

* Octave 를 이용하여 정방행렬을  LU 분해하기
** Octave 소스
   (참고: 함수 lu 를 호출하면 하삼각행렬 l, 상삼각행렬 u, 열교환 행혈 p 를 한꺼번에
            구할 수 있다.)

a = [4, 3; 6, 3]
[l, u, p] = lu(a)
p*l*u
p*l*u == a

** Octave 명령 창에서 행렬 A 의 LU 분해를 구하고 검증하기

* Python 2.7 & Numpy & Scipy를 이용하여 LU 분해 구하기
  (참고: 함수 lu_factor 를 불러 써야 하는데 이 함수는 입력 데이터와 출력 데이터를 모두
            array 타입으로 하므로, matrix 타입과 array 터입을 상호 변환하는 과정이 팔요하다.
                     matrixType = matrix(arrayType)
                     arrayType = array(matrixType)
  )

** 파이썬 소스

#!/usr/bin/env python

"""
Find a LU decomposition of a given square matrix.

Use the function 'lu_factor' of the module 'scipy.linalg.decomp_lu'.

Copyright  (c) 2011/09/18  (Sun)    PKim     pkim01 (AT) paran (DOT) com
"""

import numpy as np
import scipy as sp
from scipy import linalg
from numpy import *
from scipy.linalg.decomp_lu import *

def printWithType (x):
    print x, "   as", type(x)
   
def getLU(lu, piv):
    """
        Input
                lu      LU array
                piv    pivot array
        Output
                l       array for the lower diagonal matrix
                u      array for the upper diagonal matrix
                p      array for the permutation
    """
    l = lu.copy()
    u = lu.copy()
    for i in range(len(lu)):
        for j in range(len(lu[0])):
          if i == j:
              l[i][j] = 1
          elif i > j:
             u[i][j] = 0
          else:
             l[i][j] = 0
    p = np.eye(len(lu))
    for i in range(len(piv)):
        if i != piv[i]:
            tmp = p[piv[i]].copy()
            p[piv[i]] =  p[i].copy()
            p[i] = tmp
    return l, u, p

a = array([[4, 3], [6, 3]])
# a = array([[1, 2, 5], [2, 1, 2], [7, 1, 1]])   # This shoud work, too!
A = matrix(a)
print "A ="
printWithType(A)
print "a ="
printWithType(a)

lu, piv =lu_factor(a, False)
print "lu ="
printWithType(lu)
print "piv ="
printWithType(piv)

l, u, p = getLU(lu, piv)
print "l ="
printWithType(l)
print "u ="
printWithType(u)
print "p ="
printWithType(p)

L = matrix(l)
U = matrix(u)
P = matrix(p)

print "L ="
printWithType(L)
print "U ="
printWithType(U)
print "P ="
printWithType(P)

print "P*A =", P*A
print "L*U =", L*U
print P*A == L*U

** 위의 소스를 실행한 결과

A =
[[4 3]
 [6 3]]     as <class 'numpy.matrixlib.defmatrix.matrix'>
a =
[[4 3]
 [6 3]]     as <type 'numpy.ndarray'>
lu =
[[ 6.          3.        ]
 [ 0.66666667  1.        ]]     as <type 'numpy.ndarray'>
piv =
[1 1]     as <type 'numpy.ndarray'>
l =
[[ 1.          0.        ]
 [ 0.66666667  1.        ]]     as <type 'numpy.ndarray'>
u =
[[ 6.  3.]
 [ 0.  1.]]     as <type 'numpy.ndarray'>
p =
[[ 0.  1.]
 [ 1.  0.]]     as <type 'numpy.ndarray'>
L =
[[ 1.          0.        ]
 [ 0.66666667  1.        ]]     as <class 'numpy.matrixlib.defmatrix.matrix'>
U =
[[ 6.  3.]
 [ 0.  1.]]     as <class 'numpy.matrixlib.defmatrix.matrix'>
P =
[[ 0.  1.]
 [ 1.  0.]]     as <class 'numpy.matrixlib.defmatrix.matrix'>
P*A = [[ 6.  3.]
 [ 4.  3.]]
L*U = [[ 6.  3.]
 [ 4.  3.]]
[[ True  True]
 [ True  True]]

* Jama 패키지를 이용한 Java 애플리케이션
  (참고: JAMA is a basic linear algebra package for Java. )

** Java 소스

/**
 * Filename: TestLuDecompose.java
 *
 *  Putpose
 *            Test the LU decomposition of the Jama 1.0.2 package,
 *                 which is a famous Java Matrix package downloadable
 *                 at t http://math.nist.gov/javanumerics/jama/
 *
 *  Compile: javac -cp .;Jama-1.0.2.jar -d . TestLuDecompose.java
 *  Execute: java -classpath .;Jama-1.0.2.jar scripting.tistory.examples.Jama.TestLuDecompose
 *
 *  Author:  2011/09/19 (Mon)  Copyright (C)   pkim01 (AT) paran (DOT) com
 */

package scripting.tistory.examples.Jama;

import Jama.*;
import java.util.Date;

/*---   Example testing the LU decomposition of the JAMA package.  ---*/

public class TestLuDecompose {

   private static void print (String s) {
      System.out.print(s);
   }

   private static void printMatrix(String msg, Matrix A) {
       if (msg != null && msg.length() > 0)
           print(msg + "\n");
       int m = A.getRowDimension();
       int n = A.getColumnDimension();
       for (int i =0 ; i < m; i++) {
           for (int j =0 ; j < n; j++) {
               print(fixedWidthDoubletoString(A.get(i, j), 12, 3));
           }
           print("\n");
       }
   }

   private static void printPivot(String msg, int[] pvt) {
       if (msg != null && msg.length() > 0)
           print(msg + "[ ");
       int k = pvt.length;
       if (k > 0)
            print("" + pvt[0]);
       for (int i =1 ; i < k; i++) {
            print(", " + pvt[i]);
       }
       print(" ]\n");
   }

   private static Matrix applyPivot(Matrix A, int[] pvt) {
       int k = pvt.length;
       int m = A.getRowDimension();
       int n = A.getColumnDimension();
       double[][] t = new double[m][n];
       for (int i =0 ; i < m; i++) {
           for (int j =0 ; j < n; j++) {
              t[pvt[i]][j] = A.get(i, j);
           }
       }
        return new Matrix(t);
   }

   /** Format double with Fw.d. **/
   public static String fixedWidthDoubletoString (double x, int w, int d) {
      java.text.DecimalFormat fmt = new java.text.DecimalFormat();
      fmt.setMaximumFractionDigits(d);
      fmt.setMinimumFractionDigits(d);
      fmt.setGroupingUsed(false);
      String s = fmt.format(x);
      while (s.length() < w) {
         s = " " + s;
      }
      return s;
   }

   /** Format integer with Iw. **/
   public static String fixedWidthIntegertoString (int n, int w) {
      String s = Integer.toString(n);
      while (s.length() < w) {
         s = " " + s;
      }
      return s;
   }

   public static void main (String argv[]) {
      double eps = Math.pow(2.0,-52.0);
      int n = 2;

      double[][] a = new double[n][n];
      a[0][0] = 4;
      a[0][1] = 3;
      a[1][0] = 6;
      a[1][1] = 3;

      Matrix A = new Matrix (a);
      printMatrix("A = ", A);

      LUDecomposition LU = new LUDecomposition(A);
      Matrix L = LU.getL();
      Matrix U = LU.getU();
      int[] p = LU.getPivot();

      printMatrix("L = ", L);
      printMatrix("U = ", U);
      printPivot("pivot = ", p);
      printMatrix("L*U = ", L.times(U));

      // This should eqal to A.
      printMatrix("applyPuvot(L*U) = ", applyPivot(L.times(U), p));

      // This should eqal to A., too.
      printMatrix("t(L*U).getMatrix(p,0,n-1) = ", (L.times(U)).getMatrix(p,0,n-1));

      Matrix R = L.times(U).minus(A.getMatrix(p,0,n-1));
      printMatrix("R = ", R);

      double res = R.norm1()/(n*eps);
      print("The error is about ");
      print(fixedWidthDoubletoString(res,12, 3));
      print(".\n");
   }
}

/***********************************
  Result output:
 A =
        4.000       3.000
        6.000       3.000
 L =
        1.000       0.000
        0.667       1.000
 U =
        6.000       3.000
        0.000       1.000
 pivot = [ 1, 0 ]
 L*U =
        6.000       3.000
        4.000       3.000
 applyPuvot(L*U) =
        4.000       3.000
        6.000       3.000
 t(L*U).getMatrix(p,0,n-1) =
        4.000       3.000
        6.000       3.000
 R =
        0.000       0.000
        0.000       0.000
 The error is about        0.000.
***************************************/

* Jama 패키지를 이용한 Groovy 애플리케이션
  (참고: Groovy 언어는 Java오언어와 달리 퍼레이터 오바로딩을 지원하므로 ,
           두 행렬의 합과 곱, 또 정방행렬의 멱승을 구현하고 테스트하는 루틴을
            추가하였다. 행렬을 출력하는 printMatrix 함수도 개선하였다.)
**  Groovy 소스

/**
 * Filename: testLuDecompose.groovy
 *
 *  Putpose
 *            Test the LU decomposition of the Jama 1.0.2 package,
 *                 which is a famous Java Matrix package downloadable
 *                 at t http://math.nist.gov/javanumerics/jama/
 *
 *  Execute: groovy testLuDecompose.groovy
 *
 *  Author:  2011/09/19 (Mon)  Copyright (C)   pkim01 (AT) paran (DOT) com
 */

package scripting.tistory.examples.Jama;

import Jama.*;
import java.util.Date;

/*---   Example testing the LU decomposition of the JAMA package.  ---*/

public class MyMatrix extends Matrix {

      public MyMatrix(double[][] a) {
          super(a);
      }

      public MyMatrix(Matrix M) {
          this(M.getArray());
      }

      public Matrix multiply(Matrix B) {
          return (this as Matrix).times(B);
      }

      public Matrix multiply(double c) {
          return (this as Matrix).times(c);
      }

      public MyMatrix power(int k) {
          Matrix M = super.copy();
          if (k == 0) {
         int m = super.getRowDimension();
          int n = super.getColumnDimension();
              return new MyMatrix(Matrix.identity(m, n));
          }
          else if (k == -1) {
              if (M.det() != 0.0) {
                  return new MyMatrix(M.inverse());
              }
              return null;
          }
          else if (k < -1) {
              if (M.det() != 0.0) {
                  Matrix N =  M.inverse();
                  Matrix K =  N.copy();
              for (int i = 1 ; i < -k; i++) {
                   K = K.times(N);
               }
                  return new MyMatrix(K);
              }
              return null;
          }
          else {
       for (int i = 1 ; i < k; i++) {
               M = M.times(this);
           }
              return new MyMatrix(M);
          }
      }
}

def printMatrix(String msg, Matrix A, width=12, dwidth=3) {
   if (msg != null && msg.length() > 0)
       print(msg + "\n");
   int m = A.getRowDimension();
   int n = A.getColumnDimension();
   for (int i =0 ; i < m; i++) {
       for (int j =0 ; j < n; j++) {
           print(fixedWidthDoubletoString(A.get(i, j), width, dwidth));
       }
       print("\n");
   }
}

def printPivot(String msg, int[] pvt) {
   if (msg != null && msg.length() > 0)
       print(msg + "[ ");
   int k = pvt.length;
   if (k > 0)
    print("" + pvt[0]);
   for (int i =1 ; i < k; i++) {
    print(", " + pvt[i]);
   }
     print(" ]\n");
}

def applyPivot(Matrix A, int[] pvt) {
   int k = pvt.length;
   int m = A.getRowDimension();
   int n = A.getColumnDimension();
   double[][] t = new double[m][n];
   for (int i =0 ; i < m; i++) {
       for (int j =0 ; j < n; j++) {
           t[pvt[i]][j] = A.get(i, j);
        }
   }
   return new Matrix(t);
}

/** Format double with Fw.d. **/
def fixedWidthDoubletoString (double x, int w, int d) {
      java.text.DecimalFormat fmt = new java.text.DecimalFormat();
      fmt.setMaximumFractionDigits(d);
      fmt.setMinimumFractionDigits(d);
      fmt.setGroupingUsed(false);
      String s = fmt.format(x);
      while (s.length() < w) {
         s = " " + s;
      }
      return s;
}

/** Format integer with Iw. **/
def fixedWidthIntegertoString (int n, int w) {
      String s = Integer.toString(n);
      while (s.length() < w) {
         s = " " + s;
      }
      return s;
}

double eps = Math.pow(2.0,-52.0);
int n = 2;

double[][] a = new double[n][n];
a[0][0] = 4;
a[0][1] = 3;
a[1][0] = 6;
a[1][1] = 3;

MyMatrix A = a as MyMatrix;
printMatrix("A = ", A);

LUDecomposition LU = new LUDecomposition(A);
Matrix L = LU.getL();
Matrix U = LU.getU();
int[] p = LU.getPivot();

printMatrix("L = ", L);
printMatrix("U = ", U);
printPivot("pivot = ", p);
printMatrix("L*U = ", L.times(U));

// This should eqal to A.
printMatrix("applyPuvot(L*U) = ", applyPivot(L.times(U), p));

// This should eqal to A., too.

Matrix R = L.times(U).minus(A.getMatrix(p,0,n-1));
printMatrix("R = ", R);

double res = R.norm1()/(n*eps);
print("The error is about ");
print(fixedWidthDoubletoString(res,12, 3));
print(".\n");

printMatrix("A + A = ", A + A);
printMatrix("A*A = ", A*A);
printMatrix("A.inverse() = ", A.inverse());
printMatrix("A*A.inverse() = ", A*A.inverse());
print("A.det() = " + A.det() + "\n");
printMatrix("A*3 = ", A*3);
printMatrix("A = ", A);
printMatrix("A**2 = ", A**2);
printMatrix("A**2 = ", A**2);
printMatrix("A**0 = ", A**0);
printMatrix("A**3 = ", A**3);
printMatrix("A**4 = ", A**4);
printMatrix("A**5 = ", A**5);
printMatrix("A**6 = ", A**6);
printMatrix("A**7 = ", A**7);
printMatrix("A**8 = ", A**8, 15, 3);
printMatrix("A**9 = ", A**9, 15, 3);
printMatrix("A**10 = ", A**10, 15, 3);
printMatrix("A**(-1) = ", A**(-1));
printMatrix("A**(-2) = ", A**(-2));
printMatrix("A**(-2)*A**2 = ", A**(-2)*A**2);

/*------------------------------------
  Execution Result:
A =
       4.000       3.000
       6.000       3.000
L =
       1.000       0.000
       0.667       1.000
U =
       6.000       3.000
       0.000       1.000
pivot = [ 1, 0 ]
L*U =
       6.000       3.000
       4.000       3.000
applyPuvot(L*U) =
       4.000       3.000
       6.000       3.000
R =
       0.000       0.000
       0.000       0.000
The error is about        0.000.
A + A =
       8.000       6.000
      12.000       6.000
A*A =
      34.000      21.000
      42.000      27.000
A.inverse() =
      -0.500       0.500
       1.000      -0.667
A*A.inverse() =
       1.000       0.000
       0.000       1.000
A.det() = -6.0
A*3 =
      12.000       9.000
      18.000       9.000
A =
       4.000       3.000
       6.000       3.000
A**2 =
      34.000      21.000
      42.000      27.000
A**2 =
      34.000      21.000
      42.000      27.000
A**0 =
       1.000       0.000
       0.000       1.000
A**3 =
     262.000     165.000
     330.000     207.000
A**4 =
    2038.000    1281.000
    2562.000    1611.000
A**5 =
   15838.000    9957.000
   19914.000   12519.000
A**6 =
  123094.000   77385.000
  154770.000   97299.000
A**7 =
  956686.000  601437.000
 1202874.000  756207.000
A**8 =
    7435366.000    4674369.000
    9348738.000    5877243.000
A**9 =
   57787678.000   36329205.000
   72658410.000   45677943.000
A**10 =
  449125942.000  282350649.000
  564701298.000  355009059.000
A**(-1) =
      -0.500       0.500
       1.000      -0.667
A**(-2) =
       0.750      -0.583
      -1.167       0.944
A**(-2)*A**2 =
       1.000       0.000
       0.000       1.000
-------------------------------------------*/