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.