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

Этот вопрос связан с этим и этот один

У меня есть две матрицы полного ранга A1, A2, каждая из которых имеет размерность p x p и p-вектор y.

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

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

β2 | (β1, д, A1, A2, А1-1})

куда

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

и

β1 = (A1' A1)-1 (A1' г)

Теперь, в предыдущем вопросе здесь мне посоветовали оценить β2 с помощью подхода Холецкого, поскольку разложение Холецкого легко обновить с помощью функций 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.

r linear-algebra numerical-methods
person user189035    schedule 28.12.2011    source источник
comment
Та да! Юникод и немного разметки. Не так удобно, как LaTeX, но мы обходимся тем, что есть.   -  person joran    schedule 28.12.2011

Ответы (1)


arrow_upward
3
arrow_downward

Ваша реализация «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, когда есть более чем одно связанное решение, которое это просто модификация первого ранга.

person Jack Poulson    schedule 28.12.2011
comment
К сожалению, да, власть определила, что рендеринг LaTeX будет НАСТОЛЬКО сильно замедлять работу. Вы можете посмотреть на мою редакцию вопроса, чтобы получить несколько советов о том, как обойтись без него. Неудобно, но работает нормально. - person joran; 28.12.2011
comment
Джек Поулсон, я не нашел реализации обновления QR в R. Учитывая то, что у нас есть (реализовано обновление декомпозиции Холецкого), можете ли вы порекомендовать более эффективный подход? - person user189035; 28.12.2011
comment
Джек Поулсон:› на самом деле существует гораздо больше, чем одно обновление первого ранга A. A является общей (действительной) матрицей (не симметричной). - person user189035; 28.12.2011
comment
Я бы порекомендовал реализовать его самостоятельно или переключиться на Octave или Matlab, у которых есть рутина. Подход A'A не только менее эффективен, но и менее стабилен, и это обновление ранга 2 для A'A, несмотря на то, что кто-то сказал в вашем предыдущем вопросе. - person Jack Poulson; 29.12.2011
comment
Случайно это оказалось повторяющейся ошибкой среди статистиков (хотя очевидно, что основное утверждение неверно, учитывая мой временной эксперимент, см., например, страницу 15 здесь: stat.uiowa.edu/~luke/classes/248/homework.pdf). Кажется, предпринимаются некоторые попытки связать R и Octave. Я попробую это и дам вам знать. В любом случае спасибо за помощь и подсказки. - person user189035; 29.12.2011