Быстрое вычисление / оценка линейной системы низкого ранга

10

Линейные системы уравнений распространены в вычислительной статистике. Одна особая система, с которой я столкнулся (например, в факторном анализе), это система

Ax=b

где Здесь D - диагональная матрица n × n со строго положительной диагональю, Ω - симметричная положительная полуопределенная матрица m × mm n ), а B - произвольная n × м матрица. Нас просят решить диагональную линейную систему (легкую), которая возмущена матрицей низкого ранга. Наивный способ решить эту проблему выше, чтобы инвертировать A , используя формулу Woodbury в

A=D+BΩBT
Dn×nΩm×mmnBn×mA, Однако это не совсем правильно, поскольку факторизации Холецкого и QR обычно могут значительно ускорить решение линейных систем (и нормальных уравнений). Недавно я подошел к следующей статье , в которой, похоже, используется подход Холецкого, и упоминается численная нестабильность инверсии Вудбери. Тем не менее, документ представляется в виде черновика, и я не смог найти численных экспериментов или подтверждающих исследований. Каков уровень техники для решения описанной мной проблемы?
с промежутками
источник
1
@gappy, вы рассматривали возможность использования QR (или Cholesky) разложения для матрицы (средний член в формуле Вудберри)? Остальные операции представляют собой простые умножения матриц. Основным источником нестабильности является вычисление Ω - 1 . Поскольку т < < п Я подозреваю , это применение QR или Cholesky в сочетании с Вудбери будет быстрее , чем QR по всей матрице A . Это, конечно, не состояние дел, а общие замечания. Ω-1+ВD-1ВTΩ-1м<<NA
mpiktas
Я подозреваю, что то, что защищает Матиас Сигер, находится в от уровня техники, он очень яркий парень, и такого рода проблемы неоднократно возникают в моделях, которые он исследует. Я использую методы, основанные на Холецком, по тем же причинам. Я подозреваю, что в «Матричных вычислениях» Голуба и Ван Лоана есть обсуждение, которое является стандартным справочником для такого рода вещей (хотя у меня нет моей копии в руках). ε
Дикран Marsupial
Следует отметить , что, принимая ваша задача эквивалентна решению системы ( я + ˉ Б Ом ˉ Б Т ) х = ˉ б , где ˉ Ь = Д - 1 / 2 б . Это немного упрощает проблему. Теперь, если Σ = ˉ B Ω ˉ B T , мы знаем, что Σ положительно полуопределена с не болееВ¯знак равноD-1/2В(я+В¯ΩВ¯T)Иксзнак равноб¯б¯знак равноD-1/2бΣзнак равноВ¯ΩВ¯TΣ положительных собственных значений. Поскольку m n , нахождение m наибольших собственных значений и соответствующих собственных векторов может быть выполнено различными способами. Тогда решение имеет вид x = Q ( I + Λ ) - 1 Q T ˉ b, где Σ = Q Λ Q T дает собственное разложение Σ . мм«NмИксзнак равноQ(я+Λ)-1QTб¯Σзнак равноQΛQTΣ
кардинал
(я+В¯ΩВ¯T)D1/2Иксзнак равноб¯Иксзнак равноD-1/2Q(я+Λ)-1QTD-1/2бD1/2Иксв обоих случаях.) Обратите внимание, что все инверсии имеют диагональные матрицы и поэтому тривиальны.
кардинал
Ω-1+ВTD-1В

Ответы:

2

«Матричные вычисления» Голуба и ван Лоан подробно обсуждаются в главе 12.5.1 об обновлении факторизаций QR и Холецкого после обновлений ранга p.

Фабиан Педрегоса
источник
Я знаю, и соответствующие функции Lapack упоминаются как в статье, на которую я ссылаюсь, так и в книге. Интересно, однако, что является лучшей практикой для рассматриваемой проблемы, а не для общей проблемы обновления.
gappy