Разложение Холецкого против собственного для рисования образцов из многомерного нормального распределения

16

Я хотел бы нарисовать образец xN(0,Σ) . Википедия предлагает использовать либо разложение Холецкого, либо Собственное , то есть или Σ=D1D1TΣзнак равноQΛQT

И, следовательно, образец может быть получен через: или где Иксзнак равноD1vИксзнак равноQΛvv~N(0,я)

Википедия предполагает, что они оба одинаково хороши для генерации выборок, но метод Холецкого имеет меньшее время вычислений. Это правда? Особенно численно при использовании метода Монте-Карло, где отклонения по диагонали могут отличаться на несколько порядков? Есть ли формальный анализ по этой проблеме?

Damien
источник
1
Дэмиен, лучший рецепт, чтобы убедиться в том, какая программа быстрее, это проверить ее самостоятельно в своем программном обеспечении: функции разложения Холески и Эйгена могут различаться по скорости в разных реализациях. Путь Холецкого более популярен, AFAIK, но собственный путь может быть потенциально более гибким.
ttnphns
1
Я понимаю , Cholesky быть быстрее О(N3/3) ( Википедия ) , тогда как eigendecomposition является О(N3) ( Jacobi собственных значений Алгоритм Однако, у меня есть еще две проблемы:.? (1) Что делает "потенциально более гибким" средним и (2) Различия отличаются на несколько порядков ( 10-4 против 10-9 для самых экстремальных элементов) - имеет ли это отношение к выбранному алгоритму?
Дэмиен
@Damien Одним из аспектов «более гибкого» является то, что собственное разложение, которое для ковариационной матрицы соответствует SVD , может быть усечено для получения оптимального аппроксимации низкого ранга полной матрицы. Усеченный SVD может быть вычислен напрямую, вместо того, чтобы вычислять всю вещь, а затем выбрасывать небольшие собственные значения.
GeoMatt22
Как насчет прочтения моего ответа в Stack Overflow: получить вершины эллипса на графике ковариации эллипса (созданного car::ellipse) . Хотя вопрос задается в разных приложениях, теория остается той же. Там вы увидите красивые фигуры для геометрического объяснения.
李哲源

Ответы:

12

Проблема была изучена с помощью Straka et.al для фильтра Калмана без запаха , который рисует (детерминированные) образцы из многомерного нормального распределения , как часть алгоритма. Если повезет, результаты могут быть применимы к проблеме Монте-Карло.

Разложение Холецкого (CD) и Собственное разложение (ED) - и в этом отношении фактический корень квадратной матрицы (MSR) - это все способы, в которых положительная полуопределенная матрица (PSD) может быть разбита.

пзнак равноUSВTпзнак равноUSUTпзнак равноUSSTUTSзнак равноST

О

пзнак равноUSООTSTUTзнак равно(USО)(USО)T

О

О

  • Ознак равноя
  • Ознак равноQUSзнак равноQр
  • Ознак равноUT

Из которого после большого анализа (цитирования) в статье были сделаны следующие выводы:

  • Для преобразуемой случайной величины с некоррелированными элементами все три рассматриваемых MD обеспечивают идентичные сигма-точки и, следовательно, они практически не влияют на качество аппроксимации [Unscented Transform]. В таком случае CD может быть предпочтительным из-за его низкой стоимости.

  • Если случайная переменная содержит коррелированные элементы, использование различных [разложений] может существенно повлиять на качество аппроксимации [несцентированного преобразования] среднего значения или ковариационной матрицы преобразованной случайной величины. Два вышеупомянутых случая показали, что [ED] должен быть предпочтительным.

  • Если элементы переменной, подлежащей преобразованию, демонстрируют сильную корреляцию, так что соответствующая ковариационная матрица является почти сингулярной, необходимо принять во внимание еще одну проблему, а именно числовую устойчивость алгоритма, вычисляющего MD. SVD гораздо более численно устойчив для почти сингулярных ковариационных матриц, чем ChD.

Ссылка:

  • Straka, O .; Дуник, Дж .; Симандл М. и Хавлик Дж. «Аспекты и сравнение матричных разложений в бесценовом фильтре Калмана», American Control Conference (ACC), 2013, 2013, 3075-3080.
Damien
источник
6

Вот простая иллюстрация использования R для сравнения времени вычисления двух методов.

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

Время работы

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

При увеличении размера выборки до 10000, время выполнения

> 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 

Надеюсь это поможет.

Аарон Зенг
источник
3

Вот руководство, или демонстрация «Докажи сам себе» для бедняков:

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

Коррзнак равно[1+0,80,2+0,81+0,70,2+0,71]

Nзнак равно[[,1][,2][,3][1,]-1.08063380.65639130.8400443[2,]-1.1434241-0.1729738-0.9884772[999999,]0.48618270.03563006-2.1176976[1000000,]-0.43945511.69265517-1.9534729]

1. СВД МЕТОД:

[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. ЧОЛЕСКИЙ МЕТОД:

[Ch[с1100с21с220с+31с32с33]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 

Спасибо @ userr11852 за указание на то, что есть лучший способ вычислить разницу в производительности между SVD и Cholesky, в пользу последнего, используя функцию microbenchmark. По его предложению, вот результат:

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
Антони Пареллада
источник
@ user11852 Спасибо. Я внимательно читаю запись, microbenchmarkи это действительно имеет значение.
Антони Пареллада
Конечно, но есть ли разница в оценке производительности?
Дэмиен
Хорошая точка зрения. У меня не было времени изучить пакет.
Антони Пареллада