Cholesky versus eigendecomposition pour le prélèvement d'échantillons à partir d'une distribution normale multivariée

16

Je voudrais dessiner un échantillon . Wikipedia suggère d'utiliser soit une composition de Cholesky soit une composition Eigend , c'est-à-dire Σ = D 1 D T 1 ou Σ = Q Λ Q TxN(0,Σ)Σ=D1D1TΣ=QΛQT

Et donc l'échantillon peut être tiré via: ou x = Q x=D1vvN(0,I)x=QΛvvN(0,I)

Wikipedia suggère qu'ils sont tous deux également bons pour générer des échantillons, mais la méthode Cholesky a le temps de calcul le plus rapide. Est-ce vrai? Surtout numériquement lors de l'utilisation d'une méthode de monte-carlo, où les variances le long des diagonales peuvent différer de plusieurs ordres de grandeur? Existe-t-il une analyse formelle de ce problème?

Damien
la source
1
Damien, la meilleure recette pour vous assurer que le programme est plus rapide est de le vérifier vous-même sur votre logiciel: les fonctions de décomposition de Cholesky et Eigen peuvent varier en vitesse dans différentes implémentations. La voie Cholesky est plus populaire, AFAIK, mais la voie propre peut être potentiellement plus flexible.
ttnphns
1
Je comprends Cholesky pour être plus rapide ( Wikipedia ) alors que eigendecomposition est O ( N 3 ) ( Jacobi Eigenvalue algorithme Cependant, j'ai deux autres problèmes:.? (1) Qu'est - ce que "potentiellement plus flexible" moyenne et (2) les écarts différer de plusieurs ordres de grandeur ( 10 - quatre vs 10 - 9 pour la plupart des éléments extrêmes) - ce que cela a une influence sur l'algorithme sélectionné?O(N3/3)O(N3)104109
Damien
@Damien un aspect de "plus flexible" est que la composition de l'eigend, qui pour une matrice de covariance correspond à la SVD , peut être tronquée pour obtenir une approximation optimale de bas rang de la matrice complète. Le SVD tronqué peut être calculé directement, plutôt que de calculer l'intégralité de la chose et de rejeter ensuite les petites valeurs propres.
GeoMatt22
Que diriez-vous de lire ma réponse sur Stack Overflow: Obtenez des sommets de l'ellipse sur un tracé de covariance d'ellipse (créé par car::ellipse) . Bien que la question soit posée dans des applications différentes, la théorie sous-jacente est la même. Vous y verrez de belles figures pour une explication géométrique.
李哲源

Réponses:

12

Le problème a été étudié par Straka et al pour le filtre de Kalman non parfumé qui tire des échantillons (déterministes) d'une distribution normale multivariée dans le cadre de l'algorithme. Avec un peu de chance, les résultats pourraient être applicables au problème de monte-carlo.

La décomposition de Cholesky (CD) et la décomposition propre (ED) - et d'ailleurs la racine carrée de matrice réelle (MSR) sont toutes des façons de décomposer une matrice semi-définie positive (PSD).

Considérons la SVD d'une matrice de PSD, . C'est depuis P est PSD, en fait le même que l'ED avec P = U S U T . De plus, nous pouvons diviser la matrice diagonale par sa racine carrée: P = U P=USVTP=USUT, notant queP=USSTUT.S=ST

Nous pouvons maintenant introduire une matrice orthogonale arbitraire :O

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

Le choix de affecte en fait les performances d'estimation, en particulier lorsqu'il existe de forts éléments hors diagonale de la matrice de covariance.O

L'article a étudié trois choix d' :O

  • , ce qui correspond à l'ED;O=I
  • partir de ladécomposition QRde U O=Q, qui correspond au CD; etUS=QR
  • O=UT

D'où les conclusions suivantes ont été tirées dans le document après de nombreuses analyses (citant):

  • Pour une variable aléatoire à transformer avec des éléments non corrélés, les trois MD considérés fournissent des points sigma identiques et, par conséquent, ils ne font pratiquement aucune différence sur la qualité de l'approximation [Transformation non parfumée]. Dans un tel cas, le CD peut être préféré pour ses faibles coûts.

  • Si la variable aléatoire contient des éléments corrélés, l'utilisation de différentes [décompositions] peut affecter de manière significative la qualité de l'approximation [Transformation non parfumée] de la matrice moyenne ou de covariance de la variable aléatoire transformée. Les deux cas ci-dessus ont montré que le [DE] devrait être préféré.

  • Si les éléments de la variable à transformer présentent une forte corrélation de sorte que la matrice de covariance correspondante est presque singulière, un autre problème doit être pris en compte, qui est la stabilité numérique de l'algorithme de calcul du MD. La SVD est beaucoup plus stable numériquement pour des matrices de covariance presque singulières que la ChD.

Référence:

  • Straka, O .; Dunik, J .; Simandl, M. & Havlik, J. «Aspects and comparaison of matrix decompositions in unscented Kalman filter», American Control Conference (ACC), 2013, 2013, 3075-3080.
Damien
la source
6

Voici une illustration simple utilisant R pour comparer le temps de calcul des deux méthodes.

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

Les temps de course sont

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

Lors de l'augmentation de la taille de l'échantillon à 10000, les durées de fonctionnement sont

> 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 

J'espère que cela t'aides.

Aaron Zeng
la source
3

Voici la démonstration manuelle ou la démonstration du pauvre homme:

> 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.sept.2.sept1]

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
la source
@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