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

For a given square matrix A, if there exist two matrices Q and R such that

                                         A = QR

where Q is an orthogonal matrix (i.e transpose(Q)=inverse(Q)) and R is an upper traingular matrix. then the product QR is called a QR decomposition or QR fatorization of A


* Mathematica 에서 QR 분해하기 (분수로 계산)

* Mathematica 에서 QR 분해하기 (부동소수점수 계산)

* 64bit Ubuntu 의 Octave 에서 QR 분해하기

octave:1> a = [3, 4; 1, 2];
octave:2> [q, r] = qr(a)
q =

   0.94868  -0.31623
   0.31623   0.94868

r =

   3.16228   4.42719
   0.00000   0.63246

octave:3> q*r
ans =

   3.0000   4.0000
   1.0000   2.0000

octave:4> [q, r, p] = qr(a)
q =

   0.89443   0.44721
   0.44721  -0.89443

r =

   4.47214   3.13050
   0.00000   0.44721

p =

   0   1
   1   0

octave:5> q*r
ans =

   4   3
   2   1

octave:6> a
a =

   3   4
   1   2

octave:7> q*r*p - a
ans =

   0   0
   0   0

* 다음은 Octave 실행 창에서 help 명령으로 qr 사용법을 알아본 것이다.

octave:1> help qr
adable Function: [Q, R, P] = qr (A)
     Compute the QR factorization of A, using standard LAPACK
     subroutines.  For example, given the matrix `a = [1, 2; 3, 4]',

          [q, r] = qr (a)

     returns

          q =

            -0.31623  -0.94868
            -0.94868   0.31623

          r =

            -3.16228  -4.42719
             0.00000  -0.63246

     The `qr' factorization has applications in the solution of least
     squares problems

          `min norm(A x - b)'

     for overdetermined systems of equations (i.e., `a'  is a tall,
     thin matrix).  The QR factorization is `q * r = a' where `q' is an
     orthogonal matrix and `r' is upper triangular.

     The permuted QR factorization `[Q, R, P] = qr (A)' forms the QR
     factorization such that the diagonal entries of `r' are decreasing
     in magnitude order.  For example, given the matrix `a = [1, 2; 3,
     4]',

          [q, r, p] = qr(a)

     returns

          q =

            -0.44721  -0.89443
            -0.89443   0.44721

          r =

            -4.47214  -3.13050
             0.00000   0.44721

          p =

             0  1
             1  0

     The permuted `qr' factorization `[q, r, p] = qr (a)' factorization
     allows the construction of an orthogonal basis of `span (a)'.

Overloaded function:

   spqr (sparse matrix, ...)

   spqr (sparse complex matrix, ...)

   spqr (sparse bool matrix, ...)

qr is a built-in function

Additional help for built-in functions and operators is
available in the on-line version of the manual.  Use the command
`doc <topic>' to search the manual index.

Help and information about Octave is also available on the WWW
at http://www.octave.org and via the help@octave.org
mailing list.

* 64bit Ubuntu 의 Python 2.6.5 에서 QR 분해하기 (아래의 진한 글자만 입력)

Python 2.6.5 (r265:79063, Apr 16 2010, 13:57:41)
[GCC 4.4.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from numpy import *글
>>> import scipy as Sci
>>> import scipy.linalg
>>> a = mat("3, 4; 1, 2")
>>> q, r = scipy.linalg.qr(a)
>>> q2, r2 = matrix(q), matrix(r)
>>> q2*r2
matrix([[ 3.,  4.],
        [ 1.,  2.]])
>>> q2*r2 - a
matrix([[  0.00000000e+00,   0.00000000e+00],
        [ -1.11022302e-16,  -8.88178420e-16]])
>>> q2*r2 == a
matrix([[ True,  True],
        [False, False]], dtype=bool)

* Ububtu 에서 Python 인터프리터로 자업 중인 창을 쓰한 것

* Ubuntu 에서 Python 인터프리터로 자업 중인 화면 스샷

* Maxima 는 (아쉽지만) 아직 QR 분해를 지원하지 않는다.


* Ubuntu 에서 SciLab 을 이용하여 QR 분해하기
  (참고: SciLab 은 Ubuntu 를 설치하면 "프로그램"->"과학" 메뉴에 있는 기본 앱이다.)

a = [3, 1, 4, 2]
A = matrix(a, 2, 2)
[Q, R] = qr(A)
[Q Q'*A R]

* gsl 을 이용하여 QR 분해하는 C 소스 코드

/***********************************************
  Filename: example-05.c
  Purpose:
              Test GSL library.
              http://www.gnu.org/software/gsl/
  On Window XP & cygwin:
  On Ubuntu
        Compile: gcc -Wall -I/usr/local/include -c example-05.c
               Link: gcc -o example-05 -L/usr/local/lib example-05.o -lgsl -lgslcblas -lm

 On Mac OS X Lion:
        Compile: gcc -Wall -I/usr/local/include -c example-05.c
               Link: gcc -o example-05 -L/usr/local/lib example-05.o -lgsl -lgslcblas -lm

  Execute: ./example-05

  Result:
A is a 2 by 2 matrix
A =
[ [   3.00000000  4.00000000 ],
  [   1.00000000  2.00000000 ] ]
A =
[ [  -3.16227766 -4.42718872 ],
  [   0.16227766  0.63245553 ] ]
tau = [   1.94868330  0.00000000 ]
R =
[ [  -3.16227766 -4.42718872 ],
  [   0.00000000  0.63245553 ] ]
Q =
[ [   0.00000000  0.00000000 ],
  [   0.16227766  0.00000000 ] ]
***************************************************/

#include <stdio.h>
#include <string.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

void printVector(char *msg, gsl_vector *x, const char *fmt) {
   if (msg != NULL && strlen(msg) > 0)
       printf("%s", msg);
    size_t n = x->size;
    size_t i;
    printf("[ ");
    printf (fmt, gsl_vector_get(x, 0));
    for (i = 1; i < n; i++) {
        printf (fmt, gsl_vector_get(x, i));
    }
    printf(" ]\n");
}

void printMatrix(char *msg, gsl_matrix *a, char *fmt) {
    if (msg != NULL && strlen(msg) > 0)
        printf("%s\n", msg);
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    printf("[ [ ");
    printf (fmt, gsl_matrix_get(a, 0, 0));
    for (j = 1; j < n; j++) {
        printf (fmt, gsl_matrix_get(a, 0, j));
    }
    printf(" ]");
    if (m > 1)
        printf(",\n");
    for (i = 1; i < m; i++) {
        printf("  [ ");
        printf (fmt, gsl_matrix_get(a, i, 0));
        for (j = 1; j < n; j++) {
            printf (fmt, gsl_matrix_get(a, i, j));
    }
    printf(" ]");
    if (i < m -1)
        printf(",\n");
    }
    printf(" ]\n");
}

void product(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b) {
    size_t m = a->size1, n = a->size2, p = b->size2;
    size_t i, j, k;
    double sum;
    for (i = 0; i < m; i++) {
        for (k = 0; k < p; k++) {
            sum = 0.0;
            for (j = 0; j < n; j++) {
                sum = sum + gsl_matrix_get(a, i, j)*gsl_matrix_get(b, j, k);
            }
            gsl_matrix_set (c, i, k, sum);
        }
    }
}

void productWIthRightVector(gsl_vector *c, const gsl_matrix *a, const gsl_vector *b) {
    size_t m = a->size1, n = b->size;
    size_t i, j;
    double sum;
    for (i = 0; i < m; i++) {
        sum = 0.0;
        for (j = 0; j < n; j++) {
            sum = sum + gsl_matrix_get(a, i, j)*gsl_vector_get(b, j);
        }
        gsl_vector_set(c, i, sum);
    }
}

void productWIthLeftVector(gsl_vector *c, const gsl_vector *a, const gsl_matrix *b) {
    size_t m = a->size, n = b->size2;
    size_t i, j;
    double sum;
    for (j = 0; j < n; j++) {
        sum = 0.0;
        for (i = 0; i < m; i++) {
            sum = sum + gsl_vector_get(a, j)*gsl_matrix_get(b, i, j);
        }
        gsl_vector_set(c, j, sum);
    }
}

void add(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b) {
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    double v;
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            v = gsl_matrix_get (a, i, j) + gsl_matrix_get (b, i, j);
         gsl_matrix_set (c, i, j, v);
  }
 }
}

void subtract(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b) {
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    double v;
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            v = gsl_matrix_get (a, i, j) - gsl_matrix_get (b, i, j);
            gsl_matrix_set (c, i, j, v);
        }
    }
}

void multiply(gsl_matrix *c, double x, const gsl_matrix *a) {
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    double v;
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            v = x*gsl_matrix_get (a, i, j);
         gsl_matrix_set (c, i, j, v);
  }
 }
}

void set(gsl_matrix *c, const gsl_matrix *a) {
    gsl_matrix_memcpy(c, a);
}

double trace(const gsl_matrix *a) {
 size_t m = a->size1;
 size_t i;
 double sum = 0;
    for (i = 0; i < m; i++) {
        sum = sum + gsl_matrix_get (a, i, i);
 }
 return sum;
}

void transpose(gsl_matrix *c, const gsl_matrix *a) {
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    for (i = 0; i < n; i++) {
        for (j = 0; j < m; j++) {
         gsl_matrix_set (c, i, j, gsl_matrix_get (a, j, i));
  }
 }
}

void getQR(gsl_matrix *Q, gsl_matrix *R, const gsl_matrix *A, const gsl_vector *tau) {
   size_t m = A->size1;
   size_t n = A->size2;
   size_t i, j;
   for (i =0 ; i < m; i++) {
       for (j =0 ; j < n; j++) {
           if (i > j) {
               gsl_matrix_set(R, i, j, 0);
               gsl_matrix_set(Q, i, j, gsl_matrix_get(A, i, j));
           }
           else {
               gsl_matrix_set(R, i, j, gsl_matrix_get(A, i, j));
               gsl_matrix_set(Q, i, j, 0);
           }
       }
   }
}

int main(void)  {
    size_t n = 2;
    gsl_vector *tau = gsl_vector_alloc(n);
    // size_t i, j;
    gsl_matrix *A = gsl_matrix_alloc(n, n);
    gsl_matrix *Q = gsl_matrix_alloc(n, n);
    gsl_matrix *R = gsl_matrix_alloc(n, n);

    gsl_matrix_set (A, 0, 0, 3);
    gsl_matrix_set (A, 0, 1, 4);
    gsl_matrix_set (A, 1, 0, 1);
    gsl_matrix_set (A, 1, 1, 2);

    printf ("A is a 2 by 2 matrix\n");
    printMatrix("A = ", A, "%12.8f");

    gsl_linalg_QR_decomp(A, tau);

    printMatrix("A = ", A, "%12.8f");
    printVector("tau = ", tau, "%12.8f");

    getQR(Q, R, A, tau);
    printMatrix("R = ", R, "%12.8f");
    printMatrix("Q = ", Q, "%12.8f");

     gsl_matrix_free (Q);
     gsl_matrix_free (R);
     gsl_matrix_free (A);
     gsl_vector_free (tau);
     return 0;
}