Composição de Cholesky versus eigend para extrair amostras de uma distribuição normal multivariada

16

Gostaria de desenhar uma amostra . A Wikipedia sugere usar uma composição de Cholesky ou Eigend , por exemplo , ou xN(0,Σ)Σ=D1D1TΣ=QΛQT

E, portanto, o exemplo pode ser obtido via: ou onde x=D1vx=QΛvvN(0,I)

A Wikipedia sugere que ambos são igualmente bons para gerar amostras, mas o método de Cholesky tem o tempo de computação mais rápido. Isso é verdade? Especialmente numericamente ao usar o método de Monte Carlo, onde as variações ao longo das diagonais podem diferir em várias ordens de magnitude? Existe alguma análise formal sobre esse problema?

Damien
fonte
1
Damien, a melhor receita para garantir que programa é mais rápido é verificar você mesmo no seu software: as funções de decomposição de Cholesky e Eigen podem variar em velocidade em diferentes implementações. A maneira Cholesky é mais popular, AFAIK, mas a maneira própria pode ser potencialmente mais flexível.
ttnphns
1
Eu entendo Cholesky para ser mais rápido ( Wikipedia ), enquanto eigendecomposition é O ( N 3 ) ( Jacobi Eigenvalue Algoritmo No entanto, tenho mais dois problemas:.? (1) O que significa "potencialmente mais flexível" média e (2) os desvios diferem por várias ordens de grandeza ( 10 - 4 vs 10 - 9 para a maioria dos elementos extremos) - isto tem uma influência sobre o algoritmo seleccionado?O(N3/3)O(N3)104109
Damien
@ Damien, um aspecto de "mais flexível" é que a composição do eigend, que para uma matriz de covariância corresponde ao SVD , pode ser truncada para obter uma aproximação ideal de baixo escalão da matriz completa. O SVD truncado pode ser calculado diretamente, em vez de calcular a coisa completa e, em seguida, jogar fora os pequenos autovalores.
GeoMatt22
Que tal ler minha resposta no Stack Overflow: obtenha vértices da elipse em um gráfico de covariância da elipse (criado por car::ellipse) . Embora a pergunta seja feita em uma aplicação diferente, a teoria por trás é a mesma. Você verá figuras agradáveis ​​para explicação geométrica lá.
李哲源

Respostas:

12

O problema foi estudado por Straka et.al para o Filtro Kalman Unscented, que extrai amostras (determinísticas) de uma distribuição Normal multivariada como parte do algoritmo. Com alguma sorte, os resultados podem ser aplicáveis ​​ao problema de monte-carlo.

A decomposição de Cholesky (CD) e a decomposição de Eigen (ED) - e, nesse caso, a raiz quadrada de matriz (MSR) real são todas as maneiras pelas quais uma matriz semi-definida positiva (PSD) pode ser decomposta.

Considere o SVD de uma matriz de PSD, . Uma vez que P é PSD, isto é, na verdade, o mesmo que o ED com P = L S L T . Além disso, podemos dividir a matriz diagonal por sua raiz quadrada: P = U P=USVTP=USUT, observando queP=USSTUT.S=ST

Agora podemos introduzir uma matriz ortogonal arbitrária :O

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

A escolha de realmente afeta o desempenho da estimativa, especialmente quando há fortes elementos fora da diagonal da matriz de covariância.O

O artigo estudou três opções de :O

  • , que corresponde ao DE;O=I
  • dadecomposição QRde U O=Q, o qual corresponde ao CD; eUS=QR
  • que leva a uma matriz simétrica (isto é, MSR)O=UT

Das quais foram tiradas as seguintes conclusões no trabalho após muita análise (citação):

  • Para uma variável aleatória a ser transformada com elementos não correlacionados, todos os três MDs considerados fornecem pontos sigma idênticos e, portanto, eles quase não fazem diferença na qualidade da aproximação [Transformação sem cheiro]. Nesse caso, o CD pode ser preferido por seus baixos custos.

  • Se a variável aleatória contiver elementos correlacionados, o uso de diferentes [decomposições] poderá afetar significativamente a qualidade da aproximação [Transformada Não Centrada] da matriz de média ou covariância da variável aleatória transformada. Os dois casos acima mostraram que o [DE] deve ser preferido.

  • Se os elementos da variável a ser transformada exibem forte correlação para que a matriz de covariância correspondente seja quase singular, outra questão deve ser levada em consideração, que é a estabilidade numérica do algoritmo que calcula o MD. O SVD é muito mais numericamente estável para matrizes de covariância quase singulares que o ChD.

Referência:

  • Straka, O .; Dunik, J .; Simandl, M. & Havlik, J. "Aspectos e comparação de decomposições matriciais em filtro de Kalman sem cheiro", American Control Conference (ACC), 2013, 2013, 3075-3080.
Damien
fonte
6

Aqui está uma ilustração simples usando R para comparar o tempo de computação dos dois métodos.

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

Os tempos de execução são

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

Ao aumentar o tamanho da amostra para 10000, os tempos de execução são

> 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 

Espero que isto ajude.

Aaron Zeng
fonte
3

Aqui está a demonstração do manual, ou pobre homem, prove-me-yourself:

> 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. MÉTODO 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
fonte
@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