Cholesky frente a eigendecomposition para extraer muestras de una distribución normal multivariante

Question

Me gustaría extraer una muestra $\mathbf{x} \sim N\left(\mathbf{0}, \mathbf{\Sigma} \right)$ . Wikipedia sugiere utilizar un Cholesky o Eigendecomposition es decir $ \mathbf{\Sigma} = \mathbf{D}_1\mathbf{D}_1^T $ o $ \mathbf{\Sigma} = \mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^T $

Y, por tanto, la muestra puede extraerse a través de: $ \mathbf{x} = \mathbf{D}_1 \mathbf{v} $ o $ \mathbf{x} = \mathbf{Q}\sqrt{\mathbf{\Lambda}} \mathbf{v} $ donde $ \mathbf{v} \sim N\left(\mathbf{0}, \mathbf{I} \right) $

Wikipedia sugiere que ambos son igual de buenos para generar muestras, pero el método Cholesky tiene un tiempo de cálculo más rápido. ¿Es esto cierto? ¿Especialmente numéricamente cuando se utiliza un método monte-carlo, donde las varianzas a lo largo de las diagonales pueden diferir en varios órdenes de magnitud? ¿Existe algún análisis formal sobre este problema?

Answer 1

El problema fue estudiado por Straka et.al para la Filtro Kalman no perfeccionado que extrae muestras (deterministas) de una distribución Normal multivariante como parte del algoritmo. Con un poco de suerte, los resultados podrían ser aplicables al problema monte-carlo.

La descomposición de Cholesky (CD) y la descomposición de Eigen (ED), así como la propia Raíz cuadrada de matriz (MSR) son todas las formas en que se puede descomponer una matriz semidefinida positiva (PSD).

Considere la SVD de una matriz PSD, $P = USV^T$ . Dado que P es PSD, esto es en realidad lo mismo que la ED con $P = USU^T$ . Además, podemos dividir la matriz diagonal por su raíz cuadrada: $P = U\sqrt{S}\sqrt{S}^TU^T$ señalando que $\sqrt{S} = \sqrt{S}^T$ .

Ahora podemos introducir una matriz ortogonal arbitraria $O$ :

$P = U\sqrt{S}OO^T\sqrt{S}^TU^T = (U\sqrt{S}O)(U\sqrt{S}O)^T$ .

La elección de $O$ afecta realmente al rendimiento de la estimación, especialmente cuando hay elementos fuertemente no diagonales de la matriz de covarianza.

El documento estudiaba tres opciones de $O$ :

• $O = I$ que corresponde al DE;
• $O = Q$ del Descomposición QR de $U\sqrt{S} = QR$ que corresponde al CD; y
• $O = U^T$ lo que conduce a una matriz simétrica (es decir, MSR)

De lo cual se extrajeron las siguientes conclusiones en el documento después de mucho análisis (entrecomillado):

• Para una variable aleatoria a transformar con elementos no correlacionados, las tres MD consideradas proporcionan casi no hay diferencia en la calidad de la aproximación de la [Transformada no acentuada]. En En tal caso, se puede preferir la CD por su bajo coste.

• Si la variable aleatoria contiene elementos correlacionados, el uso de diferentes [descomposiciones] puede afectar significativamente a qu de la media o de la matriz de covarianza de la variable aleatoria transformada. Los dos casos anteriores muestran que debe preferirse la [ED].

• Si los elementos de la variable que se va a transformar presentan una fuerte correlación, de modo que la matriz de covarianza correspondiente es casi singular, hay que tener en cuenta otra cuestión, que es la estabilidad numérica. numérica del algoritmo que calcula la MD. La SVD es mucho más estable numéricamente para matrices de covarianza casi singulares que el algoritmo ChD.

Referencia:

• Straka, O.; Dunik, J.; Simandl, M. & Havlik, J. "Aspects and comparison of matrix decompositions in unscented Kalman filter", American Control Conference (ACC), 2013, 2013, 3075-3080.

Answer 2

He aquí una sencilla ilustración utilizando R para comparar el tiempo de cálculo de los dos métodos.

library(mvtnorm)
library(clusterGeneration)
set.seed(1234)
mean <- rnorm(1000, 0, 1)
sigma <- genPositiveDefMat(1000)
sigma <- sigma$Sigma

eigen.time <- system.time(
  rmvnorm(n=1000, mean=mean, sigma = sigma, method = "eigen")
  )

chol.time <- system.time(
  rmvnorm(n=1000, mean=mean, sigma = sigma, method = "chol")
  )

Los tiempos de ejecución son

> eigen.time
   user  system elapsed 
   5.16    0.06    5.33 
> chol.time
   user  system elapsed 
   1.74    0.15    1.90

Al aumentar el tamaño de la muestra a 10000, los tiempos de ejecución son

> eigen.time <- system.time(
+   rmvnorm(n=10000, mean=mean, sigma = sigma, method = "eigen")
+   )
> 
> chol.time <- system.time(
+   rmvnorm(n=10000, mean=mean, sigma = sigma, method = "chol")
+   )
> eigen.time
   user  system elapsed 
   15.74    0.28   16.19 
> chol.time
   user  system elapsed 
   11.61    0.19   11.89 

Espero que esto ayude.

Answer 3

Este es el manual, o la pobre demostración de "pruébamelo a mí mismo":

> set.seed(0)
> # The correlation matrix
> corr_matrix = matrix(cbind(1, .80, .2, .80, 1, .7, .2, .7, 1), nrow=3)
> nvar = 3 # Three columns of correlated data points
> nobs = 1e6 # One million observations for each column
> std_norm = matrix(rnorm(nvar * nobs),nrow=nobs, ncol=nvar) # N(0,1)   

$$\text{Corr}=\small \begin{bmatrix} 1 & .8 & .2\\ .8& 1 & .7 \\ .2&.7&1 \end{bmatrix}$$

$$\text{N}=\tiny \begin{bmatrix} & [,1] & [,2] & [,3] \\ [1,] & -1.0806338 & 0.6563913 & 0.8400443 \\ [2,] & -1.1434241 & -0.1729738 & -0.9884772 \\ \vdots & \vdots & \vdots & \vdots \\ \vdots & \vdots & \vdots & \vdots \\ [999999,] & 0.4861827 & 0.03563006 & -2.1176976 \\ [1000000,] & -0.4394551 & 1.69265517 & -1.9534729\\ \end{bmatrix}$$

1. MÉTODO SVD:

$$\left[ \bf \underset{[3 \times 3]}{\color{blue}{\Large\,U}}\,\,\,\,\,\underset{\tiny \begin{bmatrix}\sqrt{d_1}&0&0\\0&\sqrt{d_2}&0\\0&0&\sqrt{d_3}\end{bmatrix}}{\Large\color{blue}{\Sigma^{0.5}}} \, \underset{[3\times 10^6]}{\Large\color{blue}{N^T}} \right]^T$$

> ptm <- proc.time()
> # Singular Value Decomposition method:
> svd = svd(corr_matrix)   
> rand_data_svd = t(svd$u %*% (diag(3) * sqrt(svd$d)) %*% t(std_norm))
> proc.time() - ptm
   user  system elapsed 
   0.29    0.05    0.34 
> 
> ptm <- proc.time()

2. MÉTODO CHOLESKY:

$$\bf \left[ \underset{\begin{bmatrix}c_{11}&0&0\\c_{21}&c_{22}&0\\c_{31}&c_{32}&c_{33}\end{bmatrix}}{\Large\color{blue}{\text{Ch}}}\,\,\underset{[3\times 10^6]}{\Large\color{blue}{N^T}} \right]^T$$

> # Cholesky method:
> chole = t(chol(corr_matrix))
> rand_data_chole = t(chole %*% t(std_norm))
> proc.time() - ptm
   user  system elapsed 
   0.25    0.03    0.31 

---

Gracias a @userr11852 por indicarme que hay una forma mejor de calcular la diferencia de rendimiento entre SVD y Cholesky, a favor de esta última, utilizando la función microbenchmark . A sugerencia suya, he aquí el resultado:

microbenchmark(chol(corr_matrix), svd(corr_matrix))
Unit: microseconds
              expr     min     lq      mean  median      uq     max neval cld
 chol(corr_matrix)  24.104  25.05  28.74036  25.995  26.467  95.469   100  a 
  svd(corr_matrix) 108.701 110.12 116.27794 111.065 112.719 223.074   100   b