Sparse Matrix

Sparse Matrix는 어디서나 중요한 문제가 된다.

MATLAB/Octave

우선 MATLAB에서다음과 같은 sparse matrix format을 사용할 수 있다. 세 컬럼으로 이루어진 텍스트 파일이다. row index, col index, value로 구성한다. 예를 들어, 다음과 같다.

1 1 1
2 2 2
3 3 3
4 4 4
2 3 5
4 3 6
4 5 7
5 2 8
5 5 0

이 텍스트 파일을 mat.mtl라는 이름으로 저장했다고하자. load를 이용해 파일을 읽고 spconvert()를 이용해 sparse matrix로 바꾼다.

octave-3.4.0:38> load mat.mtl;
octave-3.4.0:39> x = spconvert(mat);

이렇게 만들어진 sparse matrix는 다음과 같다.

octave-3.4.0:40> x
x =

Compressed Column Sparse (rows = 5, cols = 5, nnz = 8 [32%])

  (1, 1) ->  1
  (2, 2) ->  2
  (5, 2) ->  8
  (2, 3) ->  5
  (3, 3) ->  3
  (4, 3) ->  6
  (4, 4) ->  4
  (4, 5) ->  7

일반적인 형렬 형태로 보면 다음과 같다.

octave-3.4.0:41> full(x)
ans =

   1   0   0   0   0
   0   2   5   0   0
   0   0   3   0   0
   0   0   6   4   7
   0   8   0   0   0

[Octave] Simple Linear Regression

Octave/MATLAB에서 단순선형회귀

모의실험: 난수생성

단순선형회귀의 간단한 예를 모의생성해 보자.
\[
Y_i = 0.5 + 2 X_{1i} + \epsilon_i, \qquad \epsilon_i \stackrel{iid}{\sim} N(0,1)
\]
와 같은 모형을 생각해 보자.

octave-3.4.0:171> x1 = [1:10]';
octave-3.4.0:171> n = length(x1);
octave-3.4.0:173> x = [ones(n,1), x1]
x =

    1    1
    1    2
    1    3
    1    4
    1    5
    1    6
    1    7
    1    8
    1    9
    1   10

octave-3.4.0:174> beta = [0.5; 2]
beta =

   0.50000
   2.00000

다음과 같이 randn() 함수를 이용하여 표준정규분포에서 난수를 생성할 수 있다. normrnd() 함수를 이용할 수도 있다. 여기에서는 동일한 실험 결과를 얻기 위하여 seed를 지정하여 randn으로 난수를 생성하였다.

octave-3.4.0:175> randn("seed", 12345);
octave-3.4.0:176> y = x * beta + randn(n,1)
y =

    3.9025
    3.6416
    6.3307
    8.7463
   11.6325
   12.6323
   14.5467
   16.9846
   18.5778
   21.1140

그래프는 다음과 같이 그릴 수 있다.

octave-3.4.0:177> plot(x1, y, "o")

회귀계수의 추정: 최소제곱법

최소제곱법(Least Square Method)를 이용한 회귀계수 추청치는 $b_1 = S_{XY}/S_{XX}$, $b_0 = \overline Y - b_1 \overline X$라는 공식을 이용하여 구할 수 있다.

octave-3.4.0:50> Sxy = sum((x1 - mean(x1)) .* (y - mean(y)))
Sxy =  165.56
octave-3.4.0:51> Sxx = sum(((x1 - mean(x1)).^2))
Sxx =  82.500
octave-3.4.0:52> b1 = Sxy / Sxx
ans =  2.0068
octave-3.4.0:53> b0 = mean(y) - b1 * mean(x1)
b0 =  0.77332

적합된 직선 $y= 0.77332 + 2.0068x$를 원데이터와 함께 그리면 다음과 같다.

octave-3.4.0:59> plot(x1, y, "o");
octave-3.4.0:60> hold on;
octave-3.4.0:61> plot(x1, b0 + b1 * x1);

행렬 표현을 이용하면 $b = (X'X)^{-1}X'y$이고 다음과 같이 계산할 수 있다.

octave-3.4.0:206> b = inv(x' * x) * x' * y
ans =

   0.77332
   2.00683

Octave에서는 Ordinary Least Square 문제를 풀 수 있는 ols(Y, X) 함수를 제공한다.

octave-3.4.0:207> ols(y, x)
ans =

   0.77332
   2.00683

적합치와 잔차

반응변수의 적합치(fitted values)는 회귀계수를 구하였으니 당연히

octave-3.4.0:54> yhat = b0 + b1 * x1;

으로 구할 수 있다. 잔차(residuals)는 관측치와 기대치의 차이이므로

octave-3.4.0:55> res = y - yhat;

로 구할 수 있다.

행렬표현을 이용하면 $\hat y = Xb=X(X'X)^{-1}X'y$로 구할 수 있다. 이때 $H=X(X'X)^{-1}X'$을 모자행렬(hat matrix)라고 한다.

octave-3.4.0:52> H = x * inv(x' * x) * x';
octave-3.4.0:52> yhat = H * y
ans =

    2.7802
    4.7870
    6.7938
    8.8007
   10.8075
   12.8143
   14.8212
   16.8280
   18.8348
   20.8417

잔차(residuals)는 $e = (I-H)y$로 구할 수 있다.

octave-3.4.0:53> n
n =  10
octave-3.4.0:54> res = (eye(n) - H) * y
ans =

   1.122329
  -1.145389
  -0.463069
  -0.054352
   0.824982
  -0.182018
  -0.274477
   0.156663
  -0.256991
   0.272321

결정계수

적합된 회귀식이 얼마나 타당한가를 알아보기 위한 것으로 결정계수(coefficient of determination)가 사용된다. 결졍계수는
\[
R^2 = \frac{SSR}{SST} = \frac{SST -SSE}{SST}
\]
이고 여기서 총제곱합(SST), 회귀제곱합(SSR), 오차제곱합(SSE)은 $SST=SSR+SSE$ 관계를 가지며 다음과 같다.
\[
SST = \sum_{i=1}^n (Y_i - \overline Y)^2, \qquad SSR = \sum_{i=1}^n (\hat Y_i - \overline Y)^2, \qquad SSE=\sum_{i=1}^n(Y_i-\hat Y_i)^2
\]
각각을 계산하면 다음과 같다.

octave-3.4.0:68> SST = sum((y - mean(y)).^2)
SST =  336.00
octave-3.4.0:70> yhat = b0 + b1*x1;
octave-3.4.0:71> SSR = sum((yhat - mean(y)).^2)
SSR =  332.26
octave-3.4.0:72> SSE = sum((y - yhat).^2)
SSE =  3.7427

결정계수 $R^2$는 다음과 같다.

octave-3.4.0:73> R2 = SSR / SST
R2 =  0.98886

분산분석

단순선형회귀의 분산분석은 귀무가설 $H_0: \beta_1 =0$ 대 $H_1: \beta_1 \neq 0$을 검정하기 위하여 귀무가설 하에
\[
F_0 = \frac{(SSR/\sigma^2)/1}{(SSE/\sigma^2)/(n-2)} = \frac{MSR}{MSE} \sim F(1, n-2)
\]
이라는 사실을 이용한다.

octave-3.4.0:75> MSR = SSR / 1
MSR =  332.26
octave-3.4.0:76> MSE = SSE / (n-2)
MSE =  0.46784
octave-3.4.0:77> F0 = MSR / MSE
ans =  710.19

유의확률(p-value)를 계산하면 다음과 같이 매우 작으므로 귀무가설을 기각한다.

octave-3.4.0:80> 1 - fcdf(F0, 1, n-2)
ans =  4.2286e-09

또는 $F(1, 8)$ 분포에서 유의수준 0.05에 해당하는 5.3177보다 $F_0=710.19$가 매우 크므로 귀무가설을 기각한다.

octave-3.4.0:81> finv(0.95, 1, n-2)
ans =  5.3177

추론

기울기에 대한 추론

기울기 $\beta_1$의 추정량 $b_1$은 $E(b_1)=\beta_1$이고 $\text{Var}(b_1) = \sigma^2/S_{XX}$인 정규분포를 따른다. 표준화하면
\[
\frac{b_1 - \beta_1}{\sigma/\sqrt{S_{XX}}} \sim N(0,1)
\]
이고 $\sigma^2$을 모르므로 대신에 추정량 $s^2$을 사용하면
\[
\frac{b_1 - \beta_1}{s/\sqrt{S_{XX}}} \sim t(n-2)
\]
이다. 여기서 분모는 표준오차 $\text{SE}(b_1)= s/\sqrt{S_{XX}}$이다.

기울기 추정량 $b_1$의 표준오차는

octave-3.4.0:88> sqrt(MSE/Sxx)
ans =  0.075305

귀무가설 $H_0: \beta_1 = 0$ 대 $H_1: \beta_1 \neq 0$을 검정하기 위하여 귀무가설 하에서 t-통계량을 계산하면

octave-3.4.0:89> t0 = b1 /sqrt(MSE/Sxx)
t0 =  26.649

이다. 유의확률이 매우 작으므로 또는 t0가 매우 크므로 귀무가설은 기각된다.

octave-3.4.0:99> 2 * (1 - tcdf(t0, n-2))
ans =  4.2286e-09
octave-3.4.0:100> tinv(0.975, n-2)
ans =  2.3060

[Octave] Matrix Algebra

Octave/MATLAB에서 행렬대수

기본 조작법

행렬 만들기

꺽쇠괄호로 묶고 행은 세미콜론, 열은 쉼표로 구분하여 입력한다. 열은 쉼표 없이 스페이스로 구분할 수도 있다. 행이나 열의 개수가 맞지 않으면 에러가 난다.

octave-3.4.0:67> A = [1, 2, 3; 4, 5, 6]
A =

   1   2   3
   4   5   6

일정한 간격의 숫자를 자동으로 생성할 때는 start:end 또는 start:step:end 형식의 문법을 이용할 수 있다.

octave-3.4.0:101> 1:10
ans =

    1    2    3    4    5    6    7    8    9   10

octave-3.4.0:102> 1:2:10
ans =

    1    3    5    7    9

0으로 채워진 행렬은 zeros(N, M), 1로 채워진 행렬은 ones(N, M)으로 생성할 수 있다.

octave-3.4.0:152> ones(4,1)
ans =

   1
   1
   1
   1

항등 행렬(identity matrix)는 eye(X)로 생성할 수 있고 대각원소 x를 가지는 대각행렬을 diag(x)로 생성할 수 있다.

octave-3.4.0:156> eye(3)
ans =

Diagonal Matrix

1 0 0
0 1 0
0 0 1

octave-3.4.0:157> diag([1 2 3])
ans =

Diagonal Matrix

1 0 0
0 2 0
0 0 3




인덱싱/슬라이싱

행렬 A에 대해 A(i,j)로 i행, j열의 원소를 가리킬 수 있다.

octave-3.4.0:92> A
A =

   1   2   3
   4   5   6

octave-3.4.0:90> A(2, 2)
ans =  5
octave-3.4.0:91> A(1, 2)
ans =  2

행렬 내의 일부 범위를 지정할 수도 있다.

octave-3.4.0:93> A(1:2, 1:2)
ans =

   1   2
   4   5

합치기: 행/열 묶기

행 벡터 x, y를 묶어 행렬을 만들 수 있다.

octave-3.4.0:116> x = 1:10
x =

    1    2    3    4    5    6    7    8    9   10

octave-3.4.0:117> y = 11:20
y =

   11   12   13   14   15   16   17   18   19   20

octave-3.4.0:118> [x; y]
ans =

    1    2    3    4    5    6    7    8    9   10
   11   12   13   14   15   16   17   18   19   20

열 벡터를 묶어 행렬을 만들 수도 있다.

octave-3.4.0:123> x = 1:3; y = 4:6; z = 7:9
z =

   7   8   9

octave-3.4.0:126> [x' y' z']
ans =

   1   4   7
   2   5   8
   3   6   9

두 행렬을 합칠 수도 있다.

octave-3.4.0:130> A
A =

   1   4   7
   2   5   8
   3   6   9

octave-3.4.0:131> B
B =

   11   12
   13   14
   15   16

octave-3.4.0:132> [A B]
ans =

    1    4    7   11   12
    2    5    8   13   14
    3    6    9   15   16

행렬의 크기

행렬 A에 대해 ndims(A) 명령을 내리면 2가 나온다. 2차원이기 때문이다. 컬럼의 수는 columns(A), 행의 수는 rows(A), 원소의 개수는 numel(A)로 알 수 있다. size(A)는 행과 열의 개수를 동시에 알려주고 length(A)는 행과 열 중 긴 것의 길이를 알려준다.

기본 연산

전치

실수 행렬 A의 전치(transpose)는 A' 또는 A.'로 가능하다.

octave-3.4.0:123> A = [1, 2, 3; 4, 5, 6]
A =

   1   2   3
   4   5   6

octave-3.4.0:124> A'
ans =

   1   4
   2   5
   3   6

참고로, 복소수 행렬 B의 경우에는 B.'는 전치, B'는 전치 켤레 행렬을 뜻한다.

행렬의 곱

두 행렬 A, B의 곱은 A*B로 계산할 수 있다.

octave-3.4.0:125> A = [1, 2, 3; 4, 5, 6]
A =

   1   2   3
   4   5   6

octave-3.4.0:126> B = [1, 1; 1, -1; -1, 1]
B =

   1   1
   1  -1
  -1   1

octave-3.4.0:127> A * B
ans =

   0   2
   3   5

행렬의 곱이 아니라 원소끼리의 곱(element by element multiplication)을 계산하려고 할 때는, A .* B와 같은 연산자를 이용한다.

octave-3.4.0:128> A = [1, 2, 3; 4, 5, 6]
A =

   1   2   3
   4   5   6

octave-3.4.0:129> B = [1, 0, 1; 0, 1, 0]
B =

   1   0   1
   0   1   0

octave-3.4.0:130> A .* B
ans =

   1   0   3
   0   5   0

트레이스

정방행렬 A의 대각원소의 합을 행렬 A의 트레이스(trace)라고 한다. trace(A)로 구할 수 있다.

octave-3.4.0:9> A = [1, 0, -1; 0, 2, 2; -1, 4, 5]
A =

   1   0  -1
   0   2   2
  -1   4   5

octave-3.4.0:10> trace(A)
ans =  8

행렬식

정방행렬 A의 행렬식(determinant)는 det(A) 또는 inverse(A)로 구할 수 있다.

octave-3.4.0:16> A = [1, 2; 3, 4]
A =

   1   2
   3   4

octave-3.4.0:17> det(A)
ans = -2

역행렬

정방행렬 A의 역행렬은 inv(A)로 구할 수 있다.

octave-3.4.0:16> A = [1, 2; 3, 4]
A =

   1   2
   3   4

octave-3.4.0:18> inv(A)
ans =

  -2.00000   1.00000
   1.50000  -0.50000

고유치와 고유벡터

정방행렬 A의 고유치(eigen value)와 고유벡터(eigenvalue)는 eig(A)로 구할 수 있다.

octave-3.4.0:30> A = [1, 4; 9, 1]
A =

   1   4
   9   1

octave-3.4.0:31> eig(A)
ans =

   7
  -5

octave-3.4.0:32> [v, lambda] = eig(A)
v =

   0.55470  -0.55470
   0.83205   0.83205

lambda =

Diagonal Matrix

   7   0
   0  -5

vec

행렬 $A$를 열벡터로 바꾸는 연산은 vec(A)로 할 수 있다.

octave-3.4.0:149> A
A =

   1   2   3
   4   5   7

octave-3.4.0:150> vec(A)
ans =

   1
   4
   2
   5
   3
   7

행렬의 유형

주어진 행렬 A의 유형을 matrix_type(A)를 이용하여 알 수 있다. unknown, full, positive definite, diagonal, permuted diagonal, upper, lower, banded, banded positive definite, singular로 구별해 준다.

octave-3.4.0:33> A = [2, 2, 1; 2, 5, 1; 1, 1, 2]
A =

   2   2   1
   2   5   1
   1   1   2

octave-3.4.0:34> matrix_type(A)
ans = Positive Definite

행렬의 분해

비정칙치 분해

행렬 $A$는
\[
A = U S V'
\]
로 분해될 수 있다. 이때 $U$, $V$는 직교행렬(orthogonal)이고 $S$는 대각행렬이다. $S$의 대각원소는 $A$의 비정칙치(singular values)라고 하며 이것은 $A'A$의 고유치의 양의 제곱근과 같다.

A =

   2   2   1
   2   5   1
   1   1   2

octave-3.4.0:180> [U, S, V] = svd(A)
U =

  -0.44910   0.29313  -0.84403
  -0.84403  -0.44910   0.29313
  -0.29313   0.84403   0.44910

S =

Diagonal Matrix

   6.41147         0         0
         0   1.81521         0
         0         0   0.77332

V =

  -0.44910   0.29313  -0.84403
  -0.84403  -0.44910   0.29313
  -0.29313   0.84403   0.44910

결과를 살펴보자. 우선, $U$와 $V$는 직교행렬이다.

octave-3.4.0:183> U' * U
ans =

   1.0000e+00  -1.9429e-16   1.6653e-16
  -1.9429e-16   1.0000e+00   2.2204e-16
   1.6653e-16   2.2204e-16   1.0000e+00

octave-3.4.0:184> V' * V
ans =

   1.00000  -0.00000   0.00000
  -0.00000   1.00000  -0.00000
   0.00000  -0.00000   1.00000

다음으로 $A=USV'$임을 확인해 보자.

octave-3.4.0:187> U * S * V'
ans =

   2.00000   2.00000   1.00000
   2.00000   5.00000   1.00000
   1.00000   1.00000   2.00000

다음으로 $A'A$의 고유치와 $S$의 대각원소의 제곱이 일치하는지 확인해 보자.

octave-3.4.0:193> eig(A' * A)
ans =

    0.59802
    3.29498
   41.10700

octave-3.4.0:192> diag(S) .** 2
ans =

   41.10700
    3.29498
    0.59802

촐레스키 분해

양정치 행렬(definite positive matrix) $A$는 항상
\[
A = R' R
\]
로 표현할 수 있으며 이때 $R$은 대각원소가 양수인 상삼각행렬(upper triangular matrix)이다. 이것을
촐레스키 분해(Cholesky decomposition)라고 하며 chol(A)로 계산할 수 있다.

octave-3.4.0:43> A
A =

   2   2   1
   2   5   1
   1   1   2

octave-3.4.0:46> matrix_type(A)
ans = Positive Definite
octave-3.4.0:44> R = chol(A)
R =

   1.41421   1.41421   0.70711
   0.00000   1.73205   0.00000
   0.00000   0.00000   1.22474

결과를 확인해 보자. 위에서 보듯이 R은 대각원소가 모두 양수이고 상삼각행렬이다. $R'R$을 계산해 보면 $A$를 얻을 수 있다.

octave-3.4.0:120> R' * R
ans =

   2.0000   2.0000   1.0000
   2.0000   5.0000   1.0000
   1.0000   1.0000   2.0000

QR 분해

임의의 실수 정방행렬 $A$의 경우에 항상
\[
A = QR
\]
로 나타낼 수 있다. 이때 $Q$는 직교행렬(orthogonal matrix)이다. 즉, $Q'Q=I$이다. $R$은 상삼각행렬(upper triangular matrix)이다. 이것을
QR decomposition이라고 하며 qr(A)로 계산할 수 있다.

octave-3.4.0:131> A = [2, 2, 1; 2, 5, 1; 1, 1, 2]
A =

   2   2   1
   2   5   1
   1   1   2

octave-3.4.0:132> [Q, R] = qr(A)
Q =

   0.66667  -0.59628  -0.44721
   0.66667   0.74536   0.00000
   0.33333  -0.29814   0.89443

R =

   3.00000   5.00000   2.00000
   0.00000   2.23607  -0.44721
   0.00000   0.00000   1.34164

결과를 확인해 보자. R은 위 결과에서 보듯이 상삼각행렬이고 다음에서 Q가 직교행렬인 것을 확인할 수 있다. 그리고 $QR =A$를 만족한다.

octave-3.4.0:134> Q' * Q
ans =

   1.00000  -0.00000   0.00000
  -0.00000   1.00000  -0.00000
   0.00000  -0.00000   1.00000

octave-3.4.0:140> Q * R
ans =

   2.00000   2.00000   1.00000
   2.00000   5.00000   1.00000
   1.00000   1.00000   2.00000

QR 분해는 정방행렬이 아닌 경우에도 가능하다.

Octave에서 확률 분포의 함수 이용하기

Octave는 MATLAB의 자유소프트웨어 판이라고 할 수 있다. 완전히는 아니지만 MATLAB의 기능을 대체할 수 있다. 문법도 MATLAB과 거의 똑같다.  그래서 Scientific Computing with MATLAB and Octave 같은 책도 있다.

Octave에서 표준정규분포의 확률밀도함수(pdf)를 그려보자. MATLAB에서도 동일하게 작동할 것이다.

octave> fplot('normpdf(x, 0, 1)', [-3,3])

누적분포함수(cdf)를 그리는 것도 마찬가지이다. normcdf를 사용하면 된다.

octave> fplot('normcdf(x, 0, 1)', [-3,3])

흔히 qunatile function이라고도 불리는 누적분포함수의 역함수도 널리 사용되는데 다음처럼 그릴 수 있다.

octave> fplot('norminv(x, 0, 1)', [0,1])

표준정규분포에서 난수를 생성하고 싶다면 normrnd(M, V, R, C)를 이용한다. 평균 M, 분산 V인 정규분포에서 R * C 난수 행렬을 만들어 준다.

octave> normrnd(0, 1, 3, 3)
ans =

     0.109137785229567    -0.643027153866707     0.718077536091356
     -1.14166438434455    -0.992687318234478     0.692413975077409
      0.18453718070133    -0.946799427044506    -0.521163772941814

함수에 대한 도움말을 보고 싶다면

octave> help('norminv')

이외에도 기본적인 확률 분포를 위한 함수들이 제공된다.

• 참조: http://www.network-theory.co.uk/docs/octave3/octave_235.html