jueves, 10 de diciembre de 2015

La descomposición de Choleski

Si \(A\) es una matriz cuadrada simétrica y definida positiva (todos sus valores propios son estrictamente positivos), entonces existen matrices \(B\) tales que \(B^2=A\).

La descomposición de Choleski de la matriz \(A\) es una matriz \(U\) triangular superior tal que \(A = U'U\).

En R, esta descomposición se obtiene con la función chol().

( m <- matrix(c(5,1,1,3),2,2) )
##      [,1] [,2]
## [1,]    5    1
## [2,]    1    3
( cm <- chol(m) )
##          [,1]      [,2]
## [1,] 2.236068 0.4472136
## [2,] 0.000000 1.6733201
t(cm) %*% cm
##      [,1] [,2]
## [1,]    5    1
## [2,]    1    3

La función crossprod() calcula el producto \(U'U\) de forma mucho más eficiente

crossprod(cm)
##      [,1] [,2]
## [1,]    5    1
## [2,]    1    3

La inversa de la matriz \(A\) se puede calcular como \(A^{-1} = U^{-1} (U')^{-1}\) o con la función chol2inv()

( mi <- chol2inv(cm) )
##             [,1]        [,2]
## [1,]  0.21428571 -0.07142857
## [2,] -0.07142857  0.35714286
m %*% mi
##              [,1]         [,2]
## [1,] 1.000000e+00 5.551115e-17
## [2,] 2.775558e-17 1.000000e+00

Resolver un sistema con la descomposición de Choleski

Supongamos un sistema lineal \(Ax=b\) donde la matriz \(A\) es simétrica y definida positiva. Entonces \[ U'Ux=b \rightarrow Ux = (U')^{-1}b\] De modo que podemos resolver el sistema en dos pasos

  1. Resolver \(U'y=b\), de modo que \(y=(U')^{-1}b\)
  2. Resolver \(Ux=y\)

En R, podemos utilizar la función solve() para resolver los sitemas lineales, pero cuando las matrices son triangulares es mejor utilizar las funciones backsolve() y forwardsolve() para matrices triangulares superiores e inferiores, respectivamente.

b <- c(11,5)
y <- forwardsolve(t(cm), b)
( x <- backsolve(cm, y) )
## [1] 2 1

Mínimos cuadrados generalizados

Supongamos que la matriz de varianzas-covarianzas de los errores es \(\text{var}(\epsilon)=\sigma^2\Sigma\), donde \(\sigma^2\) es desconocida pero \(\Sigma\) es conocida. Es decir, conocemos las varianzas y las correlaciones entre los errores, pero no su escala. Con la ayuda de la descomposición de Choleski podemos escribir \(\Sigma = U'U = SS'\), donde \(S=U'\) es triangular inferior. De este modo podemos transformar la regresión así: \[ \begin{eqnarray} y &=& X\beta + \epsilon \\ S^{-1}y &=& S^{-1}X\beta + S^{-1}\epsilon \\ \tilde{y} &=& \tilde{X}\beta + \tilde{\epsilon} \end{eqnarray} \] En este último modelo se verifica que \[ \text{var}(\tilde{\epsilon}) = \text{var}(S^{-1}\epsilon) = S^{-1}\text{var}(\epsilon)(S')^{-1} = (U')^{-1}\sigma^2 U'U U^{-1} = \sigma^2 I \] De esta forma hemos transformado un problema GLS en un caso ordinario OLS entre \(\tilde{y}\) y la matriz \(\tilde{X}\) donde los errores \(\tilde{\epsilon}\) son i.i.d. En el modelo transformado la suma de cuadrados es \[ (S^{-1}y - S^{-1}X\beta)'(S^{-1}y - S^{-1}X\beta) = (y - X\beta)'(S')^{-1}S^{-1}(y - X\beta) = (y - X\beta)'\Sigma^{-1}(y - X\beta) \] que se minimiza con \[ \hat{\beta}=(X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1}y \] de forma que \[ \text{var}(\hat{\beta}) = \sigma^2 (X'\Sigma^{-1}X)^{-1} \] Como \(\tilde{\epsilon}=S^{-1}\epsilon\), el análisis de los residuos se debe aplicar a \(S^{-1}e\). Si disponemos de la matriz \(\Sigma\) correcta, entonces estos residuos deben ser aproximadamente i.i.d. El principal problema para aplicar este método es que en la práctica desconocemos \(\Sigma\) y hay que estimarla.