Эффективное использование разложения Холецкого в R - PullRequest
Новая
       40

Эффективное использование разложения Холецкого в R

2 голосов
/ 28 декабря 2011

Этот вопрос относится к этому и этому одному

У меня есть две матрицы полного ранга A 1 , A 2 каждая с размером pxp и p-вектором y.

Эти матрицы тесно связаны в том смысле, что матрица A 2 является обновлением матрицы первого ранга A 1 .

Меня интересует вектор

β 2 |(β 1 , у, A 1 , A 2 , A 1 -1 })

где

β 2 = ( A 2 ' A 2 ) -1 ( A 2 ' y)

и

β 1 = ( A 1 ' A 1 ) -1 ( A 1 'y)

Теперь, в предыдущем вопрос здесь мне посоветовали оценить β 2 с помощью подхода Choleski, поскольку разложение Choleski легко обновить, используя функции R, такие как chud() в пакете SamplerCompare .

Ниже приведены две функции для решения линейных систем в R, первая использует функцию solve(), а вторая - подход Холецкого (вторая, которую я могу эффективно обновить). 

fx01 <- function(ll,A,y) chol2inv(chol(crossprod(A))) %*% crossprod(A,y)

fx03 <- function(ll,A,y) solve(A,y)

p <- 5
A <- matrix(rnorm(p^2),p,p)
y <- rnorm(p)

system.time(lapply(1:1000,fx01,A=A,y=y))
system.time(lapply(1:1000,fx03,A=A,y=y))

У меня вопрос: для p small обе функции кажутся сопоставимыми (на самом деле fx01еще быстрее).Но когда я увеличиваю p, fx01 становится все медленнее, поэтому при p = 100 fx03 в три раза быстрее, чем fx01.

Что вызывает ухудшение производительности fx01 и можно ли его улучшить / решить (может быть, моя реализация Choleski слишком наивна? Разве я не должен использовать функции созвездия Choleski, такие как backsolveи если да, то как?

  1. A %*% B является R-язычком для умножения матрицы A на B.
  2. crossprod(A,B) является R-язычком для A ' B (т.е. транспонировать матрицу A, умножающую матрицу / вектор B).
  3. solve(A,b) решает для x линейную систему A x = b.
  4. chol(A) - разложение Холецки матрицы PSD A .
  5. chol2inv вычислений ( X ' X ) -1 из (R части) QR-разложения X .

1 Ответ

3 голосов
/ 28 декабря 2011

Ваша реализация 'fx01', как вы упомянули, несколько наивна и выполняет гораздо больше работы, чем подход 'fx03'. В линейной алгебре (мои извинения за основной StackOverflow, не поддерживающий LaTeX!), 'Fx01' выполняет:

  • B: = A 'A примерно за n ^ 3 флопов.
  • L: = chol (B) примерно в 1/3 n ^ 3 флопа.
  • L: = inv (L) примерно в 1/3 n ^ 3 флопа.
  • B: = L 'L примерно в 1/3 n ^ 3 флопа.
  • z: = A y примерно за 2n ^ 2 флопа.
  • x: = B z примерно за 2n ^ 2 флопа.

Таким образом, стоимость выглядит очень похоже на 2n ^ 3 + 4n ^ 2, тогда как ваш подход 'fx03' использует стандартную процедуру 'решить', которая, вероятно, выполняет декомпозицию LU с частичным поворотом (2/3 n ^ 3 флопа) ) и два треугольника решает (плюс поворот) в 2n ^ 2 флопах. Таким образом, ваш подход «fx01» выполняет в три раза больше асимптотически, и это поразительно согласуется с вашими экспериментальными результатами. Обратите внимание, что если бы A был реальным симметричным или комплексным эрмитовым, то факторизация и решение LDL ^ T или LDL 'потребовали бы только вдвое меньше работы.

С учетом вышесказанного, я думаю, что вы должны заменить свое обновление Cholesky A 'A более стабильным QR-обновлением A, как я только что ответил в вашем предыдущем вопросе . Разложение QR стоит примерно 4/3 n ^ 3 флопов, и обновление ранга один для разложения QR составляет всего O (n ^ 2), поэтому этот подход имеет смысл только для общего A, когда существует более одного связанного решения, которое это просто модификация первого ранга.

...