Имея дело с инверсией положительно определенной симметричной (ковариационной) матрицы?

27

В статистике и ее различных приложениях мы часто вычисляем ковариационную матрицу , которая является положительно определенной (в рассмотренных случаях) и симметричной для различных применений. Иногда нам понадобится инверсия этой матрицы для различных вычислений (например, квадратичные формы с этой инверсией в качестве (единственной) центральной матрицы). Учитывая качества этой матрицы и предполагаемое использование, мне интересно:

Каков наилучший, с точки зрения численной стабильности, способ вычисления или использования (скажем, для квадратичных форм или умножения матрицы на вектор в целом) этого обратного? Какая факторизация может пригодиться?

Бенджамин Аллевиус
источник

Ответы:

14

Факторизация Холецкого приводит к холеризоподобной факторизации обратного с верхней треугольной матрицей .C=RTRC1=SSTS=R1

На практике лучше придерживаться обратного фактора. Если невелика, то обычно лучше сохранить неявным, так как произведения матрицы-вектора могут быть вычислены путем решения двух треугольных систем и .RSy=C1xRTz=xRy=z

Арнольд Ноймайер
источник
25

Факторизация Холецкого имеет наибольшее значение для лучшей стабильности и скорости, когда вы работаете с ковариационной матрицей, поскольку ковариационная матрица будет положительной полуопределенной симметричной матрицей. Холески здесь естественный. НО...

Если вы намереваетесь вычислить факторизацию Холецкого, прежде чем вычислять ковариационную матрицу, сделайте себе одолжение. Сделайте задачу максимально устойчивой, вычислив QR-факторизацию вашей матрицы. (QR тоже быстрый.) То есть, если вы вычислите ковариационную матрицу как

C=ATA

где удалил столбец, а затем проследите, чтобы при формировании он вычеркивал число условия. Таким образом , лучше , чтобы сформировать QR - факторы , а не в явном виде вычисления Холецкого разложения .ACAATA

A=QR

Так как Q ортогональна,

C=(QR)TQR=RTQTQR=RTIR=RTR

Таким образом, мы получаем фактор Холецкого непосредственно из QR-факторизации в форме . Если A -Меньше QR факторизации доступна, это даже лучше , так как вам не нужно . A -Меньше QR быстрая вещь , чтобы вычислить, так не генерируется. Это становится просто последовательностью преобразований Домохозяина. (Столбец с поворотом, -less QR, по логике был бы еще более стабильным, за счет некоторой дополнительной работы по выбору стержней.)RTQQQQQ

Большим преимуществом использования QR здесь является то, что он очень численно устойчив к неприятным проблемам. Опять же, это потому, что нам никогда не приходилось формировать ковариационную матрицу непосредственно для вычисления фактора Холецкого. Как только вы сформируете произведение , вы возводите в квадрат номер условия матрицы. По сути, вы теряете информацию в тех частях этой матрицы, где у вас изначально было очень мало информации для начала.ATA

И наконец, как указывает другой ответ, вам даже не нужно вообще вычислять и хранить обратное, а использовать его неявно в форме обратных вычислений в треугольных системах.

pentavalentcarbon
источник
5
И если вам нужно оценить квадратичную форму на основе , вы можете сделать это стабильно, вычислив , т. Е. Выполняется одно прямое замещение и берется норма. C1x,C1x=x,(RTR)1x=RTx2
Кристиан Клэйсон
3

Я сделал это впервые за последнее время, используя предложения mathSE.

Я думаю, что SVD было рекомендовано большинством, но я остановился на простоте Cholesky:

Если матрица , то я разлагаю на треугольную матрицу используя Холецкого, так что . Затем я использую обратную замену или прямую замену (в зависимости от того, выбираю ли я L в качестве верхнего или нижнего треугольника), чтобы инвертировать , чтобы у меня было . Из этого я могу быстро вычислить .M=AAMLM=LLLL1M1=(LL)1=LL1


Начните с:

M=AA , где известно и неявно симметрично, а также положительно определено.M

Холесская факторизация:

MLL , где квадратное и неособоеL

Back-замена:

LL1 , вероятно, самый быстрый способ инвертировать (не цитируйте меня об этом, хотя)L

Умножение:

M1=(LL)1=LL1

Используемая запись: нижние индексы - это строки, верхние индексы - это столбцы, а - транспонированиеLL1


Мой алгоритм Холецкого (возможно, из Numeric Recipes или Wikipedia)

Lij=MijMiMjMiiMiMi

Это почти можно сделать на месте (вам нужно только временное хранилище для диагональных элементов, аккумулятор и несколько целочисленных итераторов).


Мой алгоритм обратной замены (из Числовых Рецептов, проверьте их версию, поскольку я, возможно, допустил ошибку с разметкой LaTeX)

(L1)ij={1/Liiif i=j(Li(LT)j)/Liiotherwise

Поскольку появляется в выражении, важен порядок, в котором вы выполняете итерацию по матрице (некоторые части матрицы результата зависят от других ее частей, которые должны быть вычислены заранее). Проверьте код Числовые рецепты для полного примера в коде. [Изменить]: На самом деле, просто посмотрите пример Числовые рецепты. Я слишком сильно упростил, используя точечные продукты, до такой степени, что вышеприведенное уравнение имеет циклическую зависимость, независимо от того, в каком порядке вы итерируете ...LT

Марк К Коуэн
источник
2

Если вы знаете, что матрица имеет обратную (т. Е. Если она действительно положительно определена) и если она не слишком велика, то разложение Холецкого дает подходящее средство для характеристики обратной матрицы.

Вольфганг Бангерт
источник