Cholesky versus eigendecomposition für die Entnahme von Stichproben aus einer multivariaten Normalverteilung

16

Ich würde gerne eine Probe ziehen xN(0,Σ) . Wikipedia schlägt vor, entweder Cholesky oder Eigendecomposition zu verwenden , dh Σ=D1D1T oder Σ=QΛQT

Und daher kann die Probe gezogen werden über: x=D1v oder x=Q.Λv wo vN(0,ich)

Wikipedia schlägt vor, dass beide gleich gut für die Generierung von Samples geeignet sind, aber die Cholesky-Methode hat die schnellere Rechenzeit. Ist das wahr? Besonders numerisch bei Verwendung einer Monte-Carlo-Methode, bei der sich die Varianzen entlang der Diagonalen um mehrere Größenordnungen unterscheiden können? Gibt es eine formale Analyse zu diesem Problem?

Damien
quelle
1
Damien, das beste Rezept, um sicherzustellen, dass ein Programm schneller ist, ist, es selbst in Ihrer Software zu überprüfen: Cholesky- und Eigen-Dekompositionsfunktionen können in verschiedenen Implementierungen in der Geschwindigkeit variieren. Die Cholesky-Methode ist populärer, AFAIK, aber die Eigenmethode ist möglicherweise flexibler.
TTNPHNS
1
Ich verstehe Cholesky schneller sein ( Wikipedia ) , während Eigendekomposition ist O ( N 3 ) ( Jacobi Eigenwert Algorithmus Ich habe jedoch zwei weitere Probleme:.? (1) Was bedeutet "potenziell flexiblere" Mittelwert und (2) Die Varianzen unterscheiden sich um mehrere Größenordnungen ( 10 - 4 vs 10 - 9 für die extremsten Elemente) - Hat dies einen Einfluss auf den ausgewählten Algorithmus?O(N3/3)O(N3)104109
Damien
@Damien Ein Aspekt von "flexibler" ist, dass die eigendecomposition, die für eine Kovarianzmatrix der SVD entspricht , abgeschnitten werden kann, um eine optimale niedrigrangige Approximation der vollständigen Matrix zu erhalten. Die abgeschnittene SVD kann direkt berechnet werden, anstatt das Ganze zu berechnen und dann die kleinen Eigenwerte herauszuwerfen.
GeoMatt22
Wie wäre es, wenn Sie meine Antwort unter Stapelüberlauf lesen : Ermitteln Sie die Scheitelpunkte der Ellipse auf einem Ellipsen-Kovarianz-Plot (erstellt von car::ellipse) . Obwohl die Frage in verschiedenen Anwendungen gestellt wird, ist die Theorie dahinter dieselbe. Dort sehen Sie schöne Figuren zur geometrischen Erklärung.
4.

Antworten:

12

Das Problem wurde von Straka et al. Für den Unscented Kalman Filter untersucht, der als Teil des Algorithmus (deterministische) Stichproben aus einer multivariaten Normalverteilung zieht. Mit etwas Glück könnten die Ergebnisse auf das Monte-Carlo-Problem zutreffen.

Die Cholesky-Dekomposition (CD) und die Eigen-Dekomposition (ED) - und im Übrigen die eigentliche Matrix-Quadratwurzel (MSR) - sind alle Möglichkeiten, wie eine positive semi-definite Matrix (PSD) zerlegt werden kann.

Betrachten sie den SVD eine PSD - Matrix . Da P PSD ist, ist dies tatsächlich das gleiche wie der ED mit P = U S U T . Außerdem können wir die Diagonalmatrix durch ihre Quadratwurzel teilen: P = U P=USVTP=USUTunter Hinweis darauf, dassP=USSTUT.S=ST

Wir können nun eine beliebige orthogonale Matrix einführen :O

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

Die Wahl von wirkt sich tatsächlich auf die Schätzleistung aus, insbesondere wenn die Kovarianzmatrix stark von der Diagonale abweicht.O

Die Arbeit untersuchte drei Möglichkeiten von :O

  • , was der ED entspricht;O=I
  • aus derQR-Zerlegungvon U O=Q, was der CD entspricht; undUS=QR
  • O=UT was zu einer symmetrischen Matrix führt (dh MSR)

Daraus wurden nach eingehender Analyse (Zitieren) die folgenden Schlussfolgerungen gezogen:

  • Für eine zu transformierende Zufallsvariable mit unkorrelierten Elementen liefern alle drei betrachteten MDs identische Sigma-Punkte und machen daher fast keinen Unterschied in der Qualität der [Unscented Transform] -Näherung. In einem solchen Fall kann die CD wegen ihrer geringen Kosten bevorzugt werden.

  • Wenn die Zufallsvariable korrelierte Elemente enthält, kann die Verwendung verschiedener [Zerlegungen] die Qualität der [Unscented Transform] -Näherung des Mittelwerts oder der Kovarianzmatrix der transformierten Zufallsvariablen erheblich beeinflussen. Die beiden oben genannten Fälle zeigten, dass die [ED] bevorzugt werden sollte.

  • Wenn die Elemente der zu transformierenden Variablen eine starke Korrelation aufweisen, so dass die entsprechende Kovarianzmatrix nahezu singulär ist, muss ein anderes Problem berücksichtigt werden, nämlich die numerische Stabilität des Algorithmus, der die MD berechnet. Die SVD ist für nahezu singuläre Kovarianzmatrizen numerisch wesentlich stabiler als die ChD.

Referenz:

  • Straka, O .; Dunik, J .; Simandl, M. & Havlik, J. "Aspekte und Vergleich von Matrixzerlegungen in nicht parfümierten Kalman-Filtern", American Control Conference (ACC), 2013, 2013, 3075-3080.
Damien
quelle
6

Hier ist eine einfache Abbildung, in der R verwendet wird, um die Berechnungszeit der beiden Methoden zu vergleichen.

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

Die Laufzeiten sind

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

Wenn Sie die Stichprobengröße auf 10000 erhöhen, sind die Laufzeiten gleich

> 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 

Hoffe das hilft.

Aaron Zeng
quelle
3

Hier ist die manuelle oder arme Demonstration:

> 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. 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
quelle
@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