Как сделать матрицу положительно определенной?

Я пытаюсь реализовать алгоритм EM для следующей модели факторного анализа;

W_{j} = μ + B a_{j} + e_{j} for j = 1, \dots, n

$W_j = \mu+B a_j+e_j \quad\text{for}\quad j=1,\ldots,n$

где - p-мерный случайный вектор, $W_j$ $a_j$ - это q-мерный вектор скрытых переменных, а - матрица параметров pxq. $B$

В результате других предположений, использованных для модели, я знаю, что где - ковариационная матрица дисперсии ошибочных членов , = diag ( , , ..., $W_j\sim N(\mu, BB'+D)$ $D$ $e_j$ $D$ $\sigma_1^2$ $\sigma_2^2$ $\sigma_p^2$ ).

Для алгоритма EM к работе, я делаю итерации купольные с участием оценки и матриц и в течение этих итераций я вычисляя обратное на каждой итерации , используя новые оценки и . К сожалению, в ходе итераций $B$ $D$ $BB'+D$ $B$ $D$ $BB'+D$ теряет свою положительную определенность (но не должно, потому что это матрица дисперсии-ковариации), и эта ситуация разрушает сходимость алгоритма. Мои вопросы:

Показывает ли эта ситуация, что с моим алгоритмом что-то не так, поскольку вероятность должна увеличиваться на каждом шаге ЭМ?
Каковы практические способы сделать матрицу положительно определенной?

Редактировать: я вычисляю обратное с помощью леммы обращения матрицы, которая утверждает, что:

(B B^{'} + D)^{- 1} = D^{- 1} - D^{- 1} B (I_{q} + B^{'} D^{- 1} B)^{- 1} B^{'} D^{- 1}

$(BB'+D)^{-1}=D^{-1}-D^{-1}B (I_q+B'D^{-1}B)^{-1} B'D^{-1}$

где правая часть включает в себя только обратные матрицы . $q\times q$

factor-analysis expectation-maximization Энди Амос
источник

Это может помочь лучше понять, как

«теряет» свою положительную определенность. Это подразумевает, что либо

либо

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

вычисляется непосредственно из

B B^{'} + D

$BB'+D$

B B^{'}

$BB'$

D

$D$

B B^{'}

$BB'$

B

$B$ и еще сложнее, когда

вычисляется как диагональная матрица с квадратами на диагонали!

D

$D$

whuber

@whuber Обычно в FA

, поэтому

не всегда определенно положительно. Но (теоретически)

должно быть, предполагая, что

q < p

$q<p$

B B^{'}

$BB'$

B B^{'} + D

$BB' + D$

больше нуля.

σ_{j}^{2}

$\sigma^2_j$

JMS

Это связано с этим вопросом: stats.stackexchange.com/questions/6364/…

Gilead

@JMS Спасибо. Я думаю, что мой комментарий все еще уместен:

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

сравнимо с числовой ошибкой в алгоритме инверсии. Если это так, то одним из решений является применение SVD к

и обнуление действительно малых (или отрицательных) собственных значений, а затем пересчитать

B B^{'}

$BB'$

σ_{i}^{2}

$\sigma_i^2$

B B^{'}

$BB'$

и добавить

B B^{'}

$BB'$

D

$D$

whuber

Это должны быть маленькие элементы в

;

противном случае

должно быть хорошо обусловлено, поскольку

D

$D$

I_{q} + B^{'} D^{- 1} B

$I_q + B'D^{-1}B$

q < p

$q<p$

JMS

Ответы:

Хорошо, так как вы делаете FA, я предполагаю, что имеет полный ранг столбца и $B$ $q$ $q<p$ . Нам нужно еще несколько деталей, хотя. Это может быть численная проблема; это также может быть проблемой с вашими данными.

Как вы вычисляете обратное? Вам нужно обратное явно, или вы можете повторно выразить расчет как решение для линейной системы? (то есть, чтобы получить решите для x, который обычно быстрее и более стабилен) $A^{-1}b$ $Ax=b$

Что происходит с ? Действительно ли оценки малы / 0 / отрицательны? В некотором смысле это критическое звено, потому что конечно, имеет недостаток ранга и определяет сингулярную ковариационную матрицу перед добавлением , поэтому вы не можете ее инвертировать. Добавление положительной диагональной матрицы технически делает ее полным рангом, но может быть ужасно плохо обусловлен, если $D$ $BB'$ $D$ $D$ $BB'+D$ $D$ мало.

Часто оценка для идиосинкразических дисперсий (ваши , диагональные элементы ) близка к нулю или даже отрицательна; это так называемые чехлы Heywood. См., Например, http://www.technion.ac.il/docs/sas/stat/chap26/sect21.htm (любой текст FA должен обсуждать это, это очень старая и хорошо известная проблема). Это может произойти из-за неправильной спецификации модели, выбросов, неудачи, солнечных вспышек ... MLE особенно подвержен этой проблеме, поэтому, если ваш EM-алгоритм предназначен для того, чтобы вывести MLE наружу. $\sigma^2_i$ $D$

Если ваш EM-алгоритм приближается к режиму с такими оценками, это возможно для , чтобы потерять свою положительную определенность, я думаю. Есть различные решения; лично я предпочел бы байесовский подход, но даже тогда вы должны быть осторожны с вашими приорами (неправильные приоры или даже правильные приоры со слишком большой массой около 0 могут иметь ту же проблему в основном по той же причине) $BB'+D$

JMS
источник

Позвольте мне повторить, что в основной части алгоритмов вы никогда не захотите инвертировать матрицу. Вам может понадобиться в самом конце, чтобы получить стандартные оценки, хотя. Смотрите этот пост в блоге johndcook.com/blog/2010/01/19/dont-invert-that-matrix

Samsdram,

Значения матрицы D становятся все меньше и меньше с увеличением числа итераций. Может быть, это проблема, как вы указали.

Энди Амос

@ Энди Амос: Я бы поставил на это деньги. Как @whuber указывает, почти невозможно, чтобы

имел отрицательные собственные значения, если вы вычисляете его напрямую, а о нулях (из-за недостатка ранга) следует заботиться, добавляя

с его положительной диагонали - если только некоторые из них элементы действительно маленькие. Попробуйте сгенерировать некоторые данные из модели, где

B B^{'}

$BB'$

D

$D$

σ_{i}^{2}

$\sigma_i^2$

\sum_{q} B_{i q}^{2} \approx σ_{i}^{2}

$\sum_q B_{iq}^2 \approx \sigma_i^2$ . Чем больше данных, тем лучше, чтобы оценки были точными и стабильными. Это по крайней мере скажет вам, если есть проблема в вашей реализации.

JMS