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.236 0.4472
## [2,] 0.000 1.6733
t(cm) %*% cm
## [,1] [,2]
## [1,] 5 1
## [2,] 1 3
La función crossprod() calcula el productor \( 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.21429 -0.07143
## [2,] -0.07143 0.35714
m %*% mi
## [,1] [,2]
## [1,] 1.000e+00 5.551e-17
## [2,] 4.163e-17 1.000e+00
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
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
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.