Kompozycja Cholesky'ego i eigend do rysowania próbek z wielowymiarowego rozkładu normalnego

Chciałbym narysować próbkę . Wikipedia sugeruje albo stosując Cholesky'iego lub Eigendecomposition , tj i $\mathbf{x} \sim N\left(\mathbf{0}, \mathbf{\Sigma} \right)$$\mathbf{\Sigma} = \mathbf{D}_1\mathbf{D}_1^T$$\mathbf{\Sigma} = \mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^T$

Stąd przykład można pobrać za pomocą: lub gdzie \ mathbf {v} \ sim N \ left (\ mathbf {0}, \ mathbf {I} \ right) $\mathbf{x} = \mathbf{D}_1 \mathbf{v}$$\mathbf{x} = \mathbf{Q}\sqrt{\mathbf{\Lambda}} \mathbf{v}$$\mathbf{v} \sim N\left(\mathbf{0}, \mathbf{I} \right)$

Wikipedia sugeruje, że oba są równie dobre do generowania próbek, ale metoda Cholesky'ego ma szybszy czas obliczeń. Czy to prawda? Zwłaszcza numerycznie przy zastosowaniu metody Monte Carlo, gdzie wariancje wzdłuż przekątnych mogą się różnić o kilka rzędów wielkości? Czy jest jakaś formalna analiza tego problemu?

Problem został zbadany przez Straka i in. Dla Unscented Kalman Filter, który pobiera (deterministyczne) próbki z wielowymiarowego rozkładu normalnego jako części algorytmu. Przy odrobinie szczęścia wyniki mogą odnosić się do problemu Monte Carlo.

Rozkład Cholesky'ego (CD) i rozkład własny (ED) - i w tym przypadku faktyczny pierwiastek kwadratowy macierzy (MSR) to wszystkie sposoby, w jakie można rozbić dodatnią półokreśloną macierz (PSD).

Rozważmy SVD matrycy PSD, . Ponieważ P PSD, jest w rzeczywistości taka sama, jak z ED P = U S U T . Co więcej, możemy podzielić macierz diagonalną przez jej pierwiastek kwadratowy: P = U $P = USV^T$$P = USU^T$, zauważając, że$P = U\sqrt{S}\sqrt{S}^TU^T$.$\sqrt{S} = \sqrt{S}^T$

Możemy teraz wprowadzić dowolną macierz ortogonalną :$O$

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

Wybór faktycznie wpływa na wydajność estymacji, szczególnie gdy występują silne nie-diagonalne elementy macierzy kowariancji.$O$

W pracy przeanalizowano trzy opcje :$O$

• , co odpowiada ED;$O = I$
• zrozkładu QRo U $O = Q$, co odpowiada CD; i$U\sqrt{S} = QR$
• co prowadzi do macierzy symetrycznej (tj. MSR)$O = U^T$

Z których po wielu analizach (cytowaniu) wyciągnięto następujące wnioski:

• W przypadku zmiennej losowej, która ma zostać przekształcona z nieskorelowanymi elementami, wszystkie trzy rozważane MD zapewniają identyczne punkty sigma, a zatem nie mają prawie żadnej różnicy w jakości aproksymacji [Transformacja bezzapachowa]. W takim przypadku CD może być preferowane ze względu na niskie koszty.

• Jeżeli zmienna losowa zawiera skorelowane elementy, użycie różnych [rozkładów] może znacząco wpłynąć na jakość aproksymacji [Transformacja bezzapachowa] średniej lub macierzy kowariancji transformowanej zmiennej losowej. Dwa powyższe przypadki wykazały, że [ED] powinno być preferowane.

• Jeśli elementy zmiennej, która ma być transformowana, wykazują silną korelację, tak że odpowiadająca jej macierz kowariancji jest prawie pojedyncza, należy wziąć pod uwagę inny problem, którym jest stabilność liczbowa algorytmu obliczającego MD. SVD jest znacznie bardziej stabilny numerycznie dla prawie pojedynczych macierzy kowariancji niż ChD.

Odniesienie:

• Straka, O .; Dunik, J .; Simandl, M. i Havlik, J. „Aspekty i porównanie rozkładów macierzy w bezzapachowym filtrze Kalmana”, American Control Conference (ACC), 2013, 2013, 3075-3080.

Oto prosta ilustracja wykorzystująca R do porównania czasu obliczeń dwóch metod.

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")
  )

Czasy działania są następujące

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

Przy zwiększaniu wielkości próbki do 10000 czasy działania są następujące

> 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 

Mam nadzieję że to pomoże.

Oto instrukcja manualna lub pokaz biednego człowieka, który udowodni mi:

> 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. METODA 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. CHOLESKY METHOD:

$$\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 

---

Thank you to @userr11852 for pointing out to me that there is a better way to calculate the difference in performance between SVD and Cholesky, in favor of the latter, using the function microbenchmark. At his suggestion, here is the result:

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