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

16

Chciałbym narysować próbkę . Wikipedia sugeruje albo stosując Cholesky'iego lub Eigendecomposition , tj i xN.(0,Σ)Σ=D1re1T.Σ=QΛQT.

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

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?

Damien
źródło
1
Damien, najlepszą receptą, aby upewnić się, który program jest szybszy, jest sprawdzenie go samemu w swoim oprogramowaniu: Funkcje rozkładu Cholesky'ego i Eigen'a mogą różnić się prędkością w różnych implementacjach. Droga Choleskiego jest bardziej popularna, AFAIK, ale droga własna może być potencjalnie bardziej elastyczna.
ttnphns
1
Rozumiem Cholesky być szybsze O(N3/3) ( Wikipedia ), natomiast eigendecomposition jest O(N3) ( Jacobi Eigenvalue Algorytm Mam jednak dwa dalsze problemy.? (1) Co to znaczy "potencjalnie bardziej elastyczny" oznacza i (2) odchylenia różni się o kilka rzędów wielkości ( 104 vs 109 w większości elementów krańcowych) - to ma wpływ na wybranego algorytmu?
Damien
@Damien Jednym z aspektów „bardziej elastycznego” jest to, że skład eigend, który dla macierzy kowariancji odpowiada SVD , może zostać obcięty, aby uzyskać optymalne przybliżenie pełnej macierzy niskiego rzędu. Skrócone SVD można obliczyć bezpośrednio, zamiast obliczać całą rzecz, a następnie wyrzucać małe wartości własne.
GeoMatt22
Co powiesz na przeczytanie mojej odpowiedzi w Stack Overflow: Uzyskaj wierzchołki elipsy na wykresie kowariancji elipsy (utworzonym przez car::ellipse) . Chociaż pytanie jest zadawane w różnych zastosowaniach, leżąca u podstaw teoria jest taka sama. Zobaczysz tam ładne figury do geometrycznego wyjaśnienia.
李哲源

Odpowiedzi:

12

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=USVTP=USUT, zauważając, żeP=USSTUT.S=ST

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

.P=USOOTSTUT=(USO)(USO)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; iUS=QR
  • co prowadzi do macierzy symetrycznej (tj. MSR)O=UT

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.
Damien
źródło
6

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.

Aaron Zeng
źródło
3

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)   

Corr=[1.8.2.81.7.2.71]

N=[[,1][,2][,3][1,]1.08063380.65639130.8400443[2,]1.14342410.17297380.9884772[999999,]0.48618270.035630062.1176976[1000000,]0.43945511.692655171.9534729]

1. METODA SVD:

[U[3×3]Σ0.5[d1000d2000d3]NT[3×106]]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:

[Ch[c1100c21c220c31c32c33]NT[3×106]]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
Antoni Parellada
źródło
@user11852 Thank you. I read cursorily the entry on microbenchmark and it really makes a difference.
Antoni Parellada
Sure, but does it have a difference in estimation performance?
Damien
Good point. I haven't had time to explore the package.
Antoni Parellada