Komposisi Cholesky versus eigend untuk menggambar sampel dari distribusi normal multivariat

16

Saya ingin menggambar sampel . Wikipedia menyarankan untuk menggunakan komposisi Cholesky atau Eigend , yaitu atau xN(0,Σ)Σ=D1D1TΣ=QΛQT

Dan karenanya sampel dapat diambil melalui: atau mana x=D1vx=QΛvvN(0,saya)

Wikipedia menyatakan bahwa keduanya sama-sama bagus untuk menghasilkan sampel, tetapi metode Cholesky memiliki waktu komputasi yang lebih cepat. Apakah ini benar? Khususnya secara numerik ketika menggunakan metode monte-carlo, di mana varians sepanjang diagonal mungkin berbeda dengan beberapa urutan besarnya? Apakah ada analisis formal tentang masalah ini?

Damien
sumber
1
Damien, resep terbaik untuk memastikan program apa yang lebih cepat adalah dengan memeriksanya sendiri di perangkat lunak Anda: Fungsi dekomposisi Cholesky dan Eigen dapat bervariasi dalam kecepatan dalam implementasi yang berbeda. Cara Cholesky lebih populer, AFAIK, tetapi cara eigen berpotensi lebih fleksibel.
ttnphns
1
Saya mengerti Cholesky menjadi lebih cepat ( Wikipedia ) sedangkan eigendecomposition adalah O ( N 3 ) ( Jacobi Eigenvalue Algoritma Namun, saya punya dua masalah lebih lanjut:.? (1) Apa yang "berpotensi lebih fleksibel" mean dan (2) Varians berbeda dengan beberapa urutan besarnya ( 10 - 4 vs 10 - 9 untuk elemen yang paling ekstrim) - apakah ini memiliki pengaruh pada algoritma yang dipilih?O(N3/3)O(N3)104109
Damien
@ Damen salah satu aspek dari "lebih fleksibel" adalah bahwa komposisi eigend, yang untuk matriks kovarians sesuai dengan SVD , dapat dipotong untuk mendapatkan pendekatan tingkat rendah optimal dari matriks penuh. SVD terpotong dapat dihitung secara langsung, daripada menghitung semuanya dan kemudian membuang nilai eigen kecil.
GeoMatt22
Bagaimana dengan membaca jawaban saya di Stack Overflow: Dapatkan simpul dari elips pada plot kovarians elips (dibuat oleh car::ellipse) . Meskipun pertanyaannya diajukan dalam aplikasi yang berbeda, teori di baliknya adalah sama. Anda akan melihat angka-angka bagus untuk penjelasan geometris di sana.
李哲源

Jawaban:

12

Masalahnya dipelajari oleh Straka et.al untuk Filter Kalman Unscented yang mengambil (deterministik) sampel dari distribusi Normal multivariat sebagai bagian dari algoritma. Dengan sedikit keberuntungan, hasilnya mungkin berlaku untuk masalah monte-carlo.

Cholesky Decomposition (CD) dan Eigen Decomposition (ED) - dan dalam hal ini Matrix Square Root (MSR) yang sebenarnya adalah semua cara di mana matriks semi-definite positif (PSD) positif dapat dipecah.

Pertimbangkan SVD dari matriks PSD, . Karena P adalah PSD, ini sebenarnya sama dengan ED dengan P = U S U T . Selain itu, kita dapat membagi matriks diagonal dengan akar kuadratnya: P = U P=USVTP=USUT, mencatat ituP=USSTUT.S=ST

Kami sekarang dapat memperkenalkan matriks ortogonal sewenang-wenang :O

.P=USOOTSTUT=(USO)(USO)T

Pilihan sebenarnya mempengaruhi kinerja estimasi, terutama ketika ada elemen off-diagonal yang kuat dari matriks kovarians.O

Makalah ini mempelajari tiga pilihan :O

  • , yang sesuai dengan UGD;O=I
  • daridekomposisi QRdari U O=Q, yang sesuai dengan CD; danUS=QR
  • yang mengarah ke matriks simetris (yaitu MSR)O=UT

Dari mana kesimpulan berikut diambil di koran setelah banyak analisis (mengutip):

  • Untuk variabel acak yang akan diubah dengan elemen yang tidak berkorelasi, ketiga MD yang dianggap memberikan poin sigma identik dan karenanya mereka hampir tidak membuat perbedaan pada kualitas perkiraan [Unscented Transform]. Dalam kasus seperti itu CD lebih disukai karena biayanya yang rendah.

  • Jika variabel acak berisi elemen berkorelasi, penggunaan [dekomposisi] yang berbeda dapat secara signifikan mempengaruhi kualitas [Unscented Transform] perkiraan dari rata-rata atau matriks kovarians dari variabel acak yang diubah. Dua kasus di atas menunjukkan bahwa [ED] harus lebih disukai.

  • Jika elemen-elemen dari variabel yang akan diubah menunjukkan korelasi yang kuat sehingga matriks kovarian yang sesuai hampir tunggal, masalah lain harus diperhitungkan, yaitu stabilitas numerik dari algoritma yang menghitung MD. SVD jauh lebih stabil secara numerik untuk matriks kovarians yang hampir tunggal daripada CHD.

Referensi:

  • Straka, O .; Dunik, J .; Simandl, M. & Havlik, J. "Aspek dan perbandingan dekomposisi matriks dalam filter Kalman tanpa wewangian", American Control Conference (ACC), 2013, 2013, 3075-3080.
Damien
sumber
6

Berikut ini adalah ilustrasi sederhana menggunakan R untuk membandingkan waktu komputasi dari kedua metode.

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

Waktu berjalan adalah

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

Ketika meningkatkan ukuran sampel menjadi 10.000, waktu berjalan adalah

> 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 

Semoga ini membantu.

Aaron Zeng
sumber
3

Inilah manual, atau demonstrasi orang miskin, buktikan-kepada-saya sendiri:

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

Kor=[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. SVD METHOD:

[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
sumber
@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