본문으로 바로가기



R의 회귀분석계수 계산과정에 대하여 - gmp 패키지, 촐레스키(Cholesky) 분해, 그리고 QR 분해

category Statistics/Linear models 2017.08.08 14:11
R의 회귀분석계수 계산과정에 대하여 - gmp 패키지, 촐레스키 분해(Cholesky), 그리고 QR 분해

이 포스팅은 필자가 U of Iowa의 Intensive computing의 숙제를 한 내용을 정리한 포스팅이라는 점을 밝힌다. 포스팅에 사용된 패키지는 다음과 같다.

library(magrittr)
library(datasets)
library(microbenchmark)

오늘은 행렬의 분해를 사용하여 회귀분석의 계수를 구하는 과정에 대하여 알아보도록 하자. 회귀분석 계수를 구하는 방법은 여러가지가 있겠지만 그 중 촐레스키 분해(Cholesky decomposition)와 QR 분해를 살펴보도록 하겠다. 행렬의 분해와 회귀분석의 계수를 구하는 것이 관계가 있는 이유는 회귀 계수를 구하는 과정이 정규방정식(normal equation)을 푼는 문제와 같기 때문이다.

주어진 데이터 벡터를 y라고 할 때, y가 다음의 모델을 따른다고 가정하자.

\[ \begin{eqnarray*} \mathbf{y}_{n\times1} & = & X_{n\times p}\mathbf{\boldsymbol{\beta}}_{p\times1}+e_{n\times1}\\ & = & \left[\begin{array}{ccccc} \mathbf{1} & x_{1} & x_{2} & ... & x_{p}\end{array}\right]\left[\begin{array}{c} \beta_{0}\\ \beta_{1}\\ \vdots\\ \beta_{p} \end{array}\right]+e\\ & = & \beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+...+\beta_{p}x_{p}+e \end{eqnarray*} \]

이 가정하에서 선형모델의 오차항의 제곱의 합을 최소로 만드는 베타 벡터는 다음의 방정식을 만족하는 베타를 구함으로서 구할 수 있다.

\[ X^{T}X\mathbf{b}=X^{T}\mathbf{y} \]

이 내용에 대한 자세한 사항은 least square method 위키 링크와 필자의 예전 을 참고하면 좋을 것이다.

longley 데이터

실제 논문 appendix의 데이터

longley 데이터는 1967년에 J. W. Longley가 발표한 논문 An appraisal of least-squares programs from the point of view of the user의 부록에 담겨있는 데이터인데, 7개의 변수 데이터로 이루어져 있으며, 관측값은 16개 밖에 되지 않는 아주 좋은 toy example이다. 이 데이터가 유명한 이유는 변수 중 Employed를 종속변수로 회귀분석을 할 경우 변수들이 높은 상관관계를 갖는 것으로 유명하다. longley 데이터는 datasets 패키지에 담겨있으므로 다음과 같이 확인해볼 수 있다.

library(datasets)
longley
GNP.deflator GNP Unemployed Armed.Forces Population Year Employed
1947 83.00 234.29 235.60 159.00 107.61 1947 60.32
1948 88.50 259.43 232.50 145.60 108.63 1948 61.12
1949 88.20 258.05 368.20 161.60 109.77 1949 60.17
1950 89.50 284.60 335.10 165.00 110.93 1950 61.19
1951 96.20 328.98 209.90 309.90 112.08 1951 63.22
1952 98.10 347.00 193.20 359.40 113.27 1952 63.64
1953 99.00 365.38 187.00 354.70 115.09 1953 64.99
1954 100.00 363.11 357.80 335.00 116.22 1954 63.76
1955 101.20 397.47 290.40 304.80 117.39 1955 66.02
1956 104.60 419.18 282.20 285.70 118.73 1956 67.86
1957 108.40 442.77 293.60 279.80 120.44 1957 68.17
1958 110.80 444.55 468.10 263.70 121.95 1958 66.51
1959 112.60 482.70 381.30 255.20 123.37 1959 68.66
1960 114.20 502.60 393.10 251.40 125.37 1960 69.56
1961 115.70 518.17 480.60 257.20 127.85 1961 69.33
1962 116.90 554.89 400.70 282.70 130.08 1962 70.55

위의 결과에서 확인할 수 있는 사항은 원본의 데이터를 변수끼리 스케일이 비슷해지도록 맞춰놓을 것을 확인할 수 있다.

회귀계수 정확한 값 구하기

R의 gmp 패키지는 숫자를 계산할 때 소수점을 이용하는 것이 아닌 정수 두 개의 나눗셈 형식으로 표현하는 방식을 택하므로써 정확한 계산을 할 수 있게 해준다. 이것을 이용하여 solve를 이용하여 정확한 회귀 계수를 구할 수 있다.

y <- as.vector(longley[,7])
X <- as.matrix(cbind(constant = 1, longley[,-7]))

# install.packages("gmp")
library(gmp)
## 
## Attaching package: 'gmp'
## The following objects are masked from 'package:base':
## 
##     %*%, apply, crossprod, matrix, tcrossprod
# X, y rationals
r_X <- as.bigq(round(1000 * X)) / as.bigq(1000)
r_y <- as.bigq(round(1000 * y)) / as.bigq(1000)
head(r_X)
## Big Rational ('bigq') 6 x 7 matrix:
##      [,1] [,2]   [,3]        [,4]    [,5]    [,6]        [,7]
## [1,] 1    83     234289/1000 1178/5  159     13451/125   1947
## [2,] 1    177/2  129713/500  465/2   728/5   13579/125   1948
## [3,] 1    441/5  129027/500  1841/5  808/5   109773/1000 1949
## [4,] 1    179/2  284599/1000 3351/10 165     110929/1000 1950
## [5,] 1    481/5  13159/40    2099/10 3099/10 4483/40     1951
## [6,] 1    981/10 346999/1000 966/5   1797/5  11327/100   1952
# coefficients
cef_exact <- as.double(solve(t(r_X) %*% r_X,
                             (t(r_X) %*% r_y)))
cef_exact
## [1] -3.482259e+03  1.506187e-02 -3.581918e-02 -2.020230e-02 -1.033227e-02
## [6] -5.110411e-02  1.829151e+00

촐레스키(Cholesky decomposition) 분해

어떤 행렬 \(A\)가 주어졌을때, 행렬 \(A\)가 모든 \(x\)에 대하여 \(x^T \mathbf{A} x > 0\)을 만족하면 우리는 \(A\)를 양정치(postive definite) 행렬이라고 부른다. 촐레스키 분해는 이러한 양정치 행렬에 대하여 적용할 수 있는데, 만약 행렬 \(A\)가 양정치 행렬이면, 다음의 식을 만족하는 하삼각(lower triangular) 행렬이 언제나 존재한다.

\[ \mathbf{A} = \mathbf{L} \mathbf{L}^T \]

따라서 우리에게 주어진 정규방정식의 \(\mathbf{X}^T \mathbf{X}\) 부분이 양정치 행렬을 만족한다면 정규방정식을 다음과 같이 나타낼 수 있을 것이다.

\[ \begin{align*} \mathbf{X}^{T}\mathbf{X}\beta & =\mathbf{X}^{T}y\\ \Rightarrow\mathbf{L}\mathbf{L}^{T}\beta & =\mathbf{X}^{T}y\\ \tag{1}\Rightarrow\mathbf{L}\left(\mathbf{L}^{T}\beta\right) & =\mathbf{X}^{T}y \end{align*} \]

따라서 위의 식에서 \(\beta\) 벡터는 다음의 두 단계를 통해서 구할 수 있다.

  1. Solve Lz=X’y for z (forward solve lower triangular system).
  2. Solve L’b=z for b (backward solve upper triangular system)

위의 방정식을 하삼각 행렬로 정규방정식을 표현하면 좋은 이유는 풀기가 훨씬 수월하기 때문이다.

R에서의 구현

위의 과정을 앞에서 살펴본 longley 데이터를 사용하여 알아보도록 하자.

library(matrixcalc)

crossprod(X) %>% is.positive.definite
## [1] TRUE
crossprod(X) %>% eigen(symmetric = TRUE) %$% values > 0
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE

위의 \(\mathbf{X}^T \mathbf{X}\) 행렬이 먼저 양정치 행렬인지 판단하기 위해서 matrixcalc 패키지에 들어있는 is.positiv.definite 함수를 사용하던지, 아니면 그냥 기본 함수들을 이용하여 행렬(\(\mathbf{X}^T \mathbf{X}\)이 대칭행렬이므로)의 eigen value들이 모두 양수인지를 확인하면 된다. 위의 결과로 보아 우리는 촐레스키 분해를 만족하는 하삼각행렬이 존재함을 알 수 있다.

# X always includes a constant vector
regCef_chol <- function(X, y, center = FALSE){
    if (center){
        cm <- colMeans(X[, -1])
        X <- cbind(1, sweep(X[, -1], 2, cm))
    }
    
    a <- crossprod(X)
    b <- crossprod(X, y)
    
    # cholesky decomposition
    upper <- chol(a)
    
    # we need to do the following
    # z <- forwardsolve( t(upper), b )
    # but transpose can be avoid by using
    # backsolve & transpose option = TRUE
    z <- backsolve( upper, b, transpose = TRUE)
    cef <- as.vector(backsolve(upper, z))
    
    if (center) {
        cef[1] <- cef[1] - cm %*% cef[-1]
    }
    cef
}

R에서 제공하는 chol() 함수를 사용하여 촐레스키 분해를 이용한 회귀계수 구하는 함수를 정의하였다. 위의 코드에 대하여 몇가지 짚고 넘어가자.

  1. R에서 t(X) %*% X 의 코드를 쓰기보다는 crossprod(X) 함수를 써야하는데, 그 이유는 transpose를 피할 수 있기 때문이다. 이에 대한 내용은 다음의 stackoverflow를 참조하자.
  2. chol() 함수는 결과값으로 상삼각(upper triangle) 행렬을 사용하고 있었다. 따라서 코드는 위에서 살펴본 식과 살짝 다르다.
  3. regCef_chol() 함수는 centering을 이용하여 구하는 방법과 그렇지 않는 방법 두 가지를 구현하고 있다. center 옵션을 FALSE로 설정한 경우, 데이터 행렬을 그대로 사용하고, TRUE로 설정한 경우는 데이터 행렬을 각 변수별 평균값을 뺀, 즉, centered된 행렬을 사용하여 구한 후 다시 원래 데이터를 사용하여 구한 계수의 형식으로 복원시킨다. 이것과 관련한 내용은 필자의 이전 능선회귀분석의 coefficient path 그리기를 참고하면 좋을 것이다.
  4. 위의 식 (1)을 그대로 구현하려면 forwardsolve와 backsolve를 이용해야 하지만, transpose 옵션을 이용한 트릭을 사용하면 backsolve함수를 두 번 이용하여 transpose를 피할 수 있다. 아래의 코드를 내용을 참고하도록 하자.
a <- crossprod(X)
b <- crossprod(X, y)

# cholesky decomposition
upper <- chol(a)
microbenchmark(
    forwardsolve( t(upper), b ),
    backsolve( upper, b, transpose = TRUE)
)
## Unit: microseconds
##                                   expr   min    lq     mean median    uq
##              forwardsolve(t(upper), b) 9.331 9.797 10.42754  9.798 9.798
##  backsolve(upper, b, transpose = TRUE) 5.132 5.599  6.15877  5.599 6.065
##     max neval
##  47.120   100
##  34.525   100
z1 <- forwardsolve( t(upper), b )
z2 <- backsolve( upper, b, transpose = TRUE)
z1 == z2
##      [,1]
## [1,] TRUE
## [2,] TRUE
## [3,] TRUE
## [4,] TRUE
## [5,] TRUE
## [6,] TRUE
## [7,] TRUE

자, 이제 우리가 정의한 함수를 사용하여 계수를 구해보도록 하자.

cef1 <- regCef_chol( X, y )
cef2 <- regCef_chol( X, y, center = TRUE)

# result test
cef1
## [1] -3.482259e+03  1.506187e-02 -3.581918e-02 -2.020230e-02 -1.033227e-02
## [6] -5.110411e-02  1.829151e+00
cef2
## [1] -3.482259e+03  1.506187e-02 -3.581918e-02 -2.020230e-02 -1.033227e-02
## [6] -5.110411e-02  1.829151e+00
cef1 == cef2
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE

재미있는 사실은 cef1과 cef2가 콘솔상에는 똑같은 숫자로 표시되어 같은 벡터처럼 보이지만, 실제 두값을 비교해보았을 때에는 다르다는 결과가 나온다.

앞에서 gmp구한 cc_exact와는 얼마나 차이가 나는지 유클리디안 놈을 사용하여 구해보자.

(cef_exact - cef1)^2 %>% sum()
## [1] 1.219947e-09
(cef_exact - cef2)^2 %>% sum()
## [1] 6.784948e-18

센터링을 사용하여 구한 계수가 gmp를 사용하여 구한 정확한 OLS 회귀계수와 상당히 가까움을 알 수 있다.

QR 분해

\(n \times p\)의 실수 행렬 \(A\)가 주어졌을때(\(n \geq p\)), 우리는 \(A\)를 항상 다음과 같이 분해할 수 있다.

\[ A = QR \] 위의 식에서 \(Q\)\(n \times p\)의 행렬이 되고, \(Q^TQ=I_p\)을 만족하는 행렬이 된다. 또한 \(R\)은 상삼각행렬(upper triangular matrix)가 된다. 만약, 행렬 \(R\)의 대각 원소들이 모두 0이 아니라면, 상삼각행렬 \(R\)은 역행렬을 가지게 된다.

QR 분해를 이용한 회귀계수를 구하는 방법은 앞에서 살펴본 촐레스키 분해보다 좀 더 선호되는데 그 이유는 다음과 같다.

  1. QR분해의 경우 일반적인 행렬에 대하여 언제나 존재하고, 회귀 계수를 구하는 과정에 있어서 좀 더 안정적이다.
  2. 결정적으로 \(\mathbf{X}^T \mathbf{X}\) 행렬을 구하지 않아도 된다는 것에 있다. 그 이유을 살펴보도록 하자.

\[ \begin{eqnarray*} X^{T}X\mathbf{\beta} & = & X^{T}\mathbf{y}\\ \left(QR\right)^{T}\left(QR\right)\mathbf{\beta} & = & \left(QR\right)^{T}\mathbf{y}\\ R^{T}Q^{T}QR\mathbf{\beta} & = & R^{T}Q^{T}\mathbf{y}\\ R^{T}R\mathbf{\beta} & = & R^{T}Q^{T}\mathbf{y} \end{eqnarray*} \]

여기서 만약 우리가 QR 분해를 했을 경우, 상삼각행렬의 대각원소가 0이 아니라면, 역행렬이 존재하므로, 위의 식은 다음과 같이 된다는 것을 알 수 있다.

\[ R\mathbf{\beta}=Q^{T}\mathbf{y} \]

따라서 베타 벡터를 구하는 데에 있어서 데이터 매트릭스를 QR분해를 한번 적용한 후, 분해된 두 매트릭스를 사용하여 backsolve를 적용하면 되기 때문에 촐레스키 분해를 사용한 것보다 훨씬 계산이 간편해지게 되는 것이다.

R에서의 QR 분해

QR_X <- qr(X)
QR_X$qr %>% dim
## [1] 16  7

R에서 제공하는 qr 함수를 사용하여 X를 분해해보면 좀 당황할 수 있다. 왜냐하면 Q 행렬과 R 행렬이 아닌 그냥 사각행렬이 나오기 때문이다. 하지만 설명서를 읽어보면 결과값으로 나온 사각행렬에 Q행렬과 R행렬의 정보가 같이 들어있다는 것을 알 수 있다.

    The upper triangle contains the $\bold{R}$ of the decomposition and the lower triangle contains information on the $\bold{Q}$ of the decomposition (stored in compact form)

이러한 결과값을 우리가 확인하기 위해서는 추가적인 함수를 사용해줘야 한다.

qr.Q(QR_X) %>% dim
## [1] 16  7
qr.R(QR_X)
##      constant GNP.deflator        GNP  Unemployed Armed.Forces
## 1947       -4   -406.72500 -1550.7938 -1277.32500  -1042.67500
## 1948        0     41.79551   381.7172   224.61743    125.26181
## 1949        0      0.00000    49.8229   -31.18589    -29.98495
## 1950        0      0.00000     0.0000  -282.06021    164.42555
## 1951        0      0.00000     0.0000     0.00000   -170.35326
## 1952        0      0.00000     0.0000     0.00000      0.00000
## 1953        0      0.00000     0.0000     0.00000      0.00000
##         Population          Year
## 1947 -469.69600000 -7818.0000000
## 1948   26.37951035    18.2758880
## 1949    4.19691804     1.7752596
## 1950   -3.18974111    -1.4529877
## 1951    0.04624846    -0.4491762
## 1952    1.46320173    -0.2819006
## 1953    0.00000000    -0.6693051

따라서 우리가 위에서 구한 회귀계수를 구하는 방법은 다음과 같이 구할 수 있을 것이다.

backsolve(qr.R(QR_X), crossprod(qr.Q(QR_X), y)) %>% as.vector()
## [1] -3.482259e+03  1.506187e-02 -3.581918e-02 -2.020230e-02 -1.033227e-02
## [6] -5.110411e-02  1.829151e+00

하지만 R에서는 다음과 같은 간단한 함수를 제공한다.

cef_qr <- qr.solve(X, y)
cef_qr
##      constant  GNP.deflator           GNP    Unemployed  Armed.Forces 
## -3.482259e+03  1.506187e-02 -3.581918e-02 -2.020230e-02 -1.033227e-02 
##    Population          Year 
## -5.110411e-02  1.829151e+00

gmp를 사용한 계수와의 차를 이용하여 앞에서 구한 촐레스키 분해를 사용한 방법과 정확도를 비교해 보도록 하자. 센터링을 사용한 촐레스키 분해를 이용한 방법보다 더 정확한 값을 갖는 것을 알 수 있다. 이러한 이유로 R에서 제공하는 lm함수는 QR 분해를 이용하여 회귀계수를 구하고 있다.

(cef_exact - cef1)
## [1] -3.492773e-05  2.954203e-10 -9.547554e-10 -1.460293e-10 -4.604652e-11
## [6]  2.386142e-09  1.793108e-08
(cef_exact - cef2)
## [1] -2.604793e-09  2.344444e-13 -1.148734e-13 -1.602191e-14 -4.739265e-15
## [6]  5.411088e-13  1.314060e-12
(cef_exact - cef_qr)
##      constant  GNP.deflator           GNP    Unemployed  Armed.Forces 
## -5.911716e-12  6.002143e-16 -9.714451e-17 -1.734723e-17 -3.469447e-17 
##    Population          Year 
## -1.422473e-15  3.108624e-15
cef_lm <- lm(y~0+X)$coefficients
cef_lm == cef_qr
##     Xconstant XGNP.deflator          XGNP   XUnemployed XArmed.Forces 
##          TRUE          TRUE          TRUE          TRUE          TRUE 
##   XPopulation         XYear 
##          TRUE          TRUE

Reference

[1] Luke Tierney, STAT:7400 Computer Intensive Statistics Note: http://homepage.divms.uiowa.edu/~luke/classes/STAT7400/notes.pdf

[2] J. W. Longley (1967) An appraisal of least-squares programs from the point of view of the user. Journal of the American Statistical Association 62, 819–841.

[3] The QR Decomposition of a Matrix in R: http://astrostatistics.psu.edu/su07/R/html/base/html/qr.html

[4] Monahan, John F. A primer on linear models. CRC Press, 2008.


SHARE TO

신고


티스토리 툴바