Como tornar uma matriz positiva definitiva?

Estou tentando implementar um algoritmo EM para o seguinte modelo de análise fatorial;

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$

onde é um vetor aleatório p-dimensional, é um vetor q-dimensional de variáveis latentes e é uma matriz de parâmetros pxq. $W_j$ $a_j$ $B$

Como resultado de outras suposições usadas para o modelo, eu sei que que é a matriz de covariância de variância dos termos de erro , = diag ( , , ..., ). $W_j\sim N(\mu, BB'+D)$ $D$ $e_j$ $D$ $\sigma_1^2$ $\sigma_2^2$ $\sigma_p^2$

Para o algoritmo EM para o trabalho, eu estou fazendo iterações cúpula envolvendo estimativa de $B$ e $D$ matrizes e durante estes iterações eu estou calculando o inverso do $BB'+D$ em cada iteração usando novas estimativas de $B$ e $D$ . Infelizmente, durante o curso das iterações, $BB'+D$ perde sua definição positiva (mas não deve ser porque é uma matriz de variância-covariância) e essa situação arruina a convergência do algoritmo. Minhas perguntas são:

Essa situação mostra que há algo errado com meu algoritmo, pois a probabilidade deve aumentar a cada etapa do EM?
Quais são as formas práticas de tornar uma matriz positiva definitiva?

Edit: Estou computando o inverso usando um lema de inversão de matriz que afirma que:

(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}$

onde o lado direito envolve apenas o inverso de matrizes $q\times q$ .

factor-analysis expectation-maximization Andy Amos
fonte

Pode ajudar a entender melhor como "perde" sua definição positiva. Isso implica que ou (ou ambos) estão se tornando definitivos não positivos. Isso é difícil de fazer quando é calculado diretamente de e ainda mais difícil quando é calculado como uma matriz diagonal com quadrados na diagonal!

B B^{'} + D

$BB'+D$

B B^{'}

$BB'$

D

$D$

B B^{'}

$BB'$

B

$B$

D

$D$

whuber

@whuber Normalmente em FA , então nem sempre é positivo. Mas (teoricamente) deveria ser, assumindo que os são todos maiores que zero.

q < p

$q<p$

B B^{'}

$BB'$

B B^{'} + D

$BB' + D$

σ_{j}^{2}

$\sigma^2_j$

JMS

Isso está relacionado a esta pergunta: stats.stackexchange.com/questions/6364/…

Gilead

@JMS Obrigado. Acho que meu comentário ainda é pertinente: pode ser indefinido, mas ainda não deve ter autovalores negativos. Problemas surgirão quando o menor dos for comparável ao erro numérico no algoritmo de inversão. Se este for o caso, uma solução é aplicar SVD para e zero, os realmente pequenos (ou negativos) valores próprios, em seguida, recalcular e adicione .

B B^{'}

$BB'$

σ_{i}^{2}

$\sigma_i^2$

B B^{'}

$BB'$

B B^{'}

$BB'$

D

$D$

whuber

Tem que haver pequenos elementos em ; deve estar bem condicionado, caso contrário, uma vez que

D

$D$

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

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

q < p

$q<p$

JMS

Respostas:

OK, já que você está fazendo FA, estou assumindo que é da coluna completa e . Precisamos de mais alguns detalhes. Isso pode ser um problema numérico; também pode ser um problema com seus dados. $B$ $q$ $q<p$

Como você está computando o inverso? Você precisa do inverso explicitamente ou pode reexprimir o cálculo como a solução para um sistema linear? (ou seja, para obter resolva para x, que normalmente é mais rápido e mais estável) $A^{-1}b$ $Ax=b$

O que está acontecendo com ? As estimativas são realmente pequenas / 0 / negativas? Em certo sentido, é o elo crítico, porque é obviamente deficiente na classificação e define uma matriz de covariância singular antes de adicionar , portanto você não pode invertê-lo. A adição da matriz diagonal positiva tecnicamente a torna completa, mas ainda pode estar terrivelmente mal condicionada se for pequeno. $D$ $BB'$ $D$ $D$ $BB'+D$ $D$

Muitas vezes, a estimativa para as variações idiossincráticas (your , os elementos diagonais de ) é próxima de zero ou até negativa; esses são chamados casos de Heywood. Veja, por exemplo, http://www.technion.ac.il/docs/sas/stat/chap26/sect21.htm (qualquer texto da FA também deve discutir isso, é um problema muito antigo e bem conhecido). Isso pode resultar de erros de especificação do modelo, outliers, má sorte, erupções solares ... o MLE é particularmente propenso a esse problema; portanto, se o seu algoritmo EM for projetado para obter a aparência do MLE. $\sigma^2_i$ $D$

Se o seu algoritmo EM estiver se aproximando de um modo com essas estimativas, é possível que o perca sua definição positiva, eu acho. Existem várias soluções; pessoalmente, eu preferiria uma abordagem bayesiana, mas mesmo assim você precisa ter cuidado com seus priores (priores impróprios ou priores apropriados com muita massa perto de 0 podem ter o mesmo problema basicamente pelo mesmo motivo) $BB'+D$

JMS
fonte

Deixe-me dizer que na parte principal dos algoritmos, você nunca deseja inverter uma matriz. Você pode precisar no final, para obter as estimativas padrão. Veja este post no blog johndcook.com/blog/2010/01/19/dont-invert-that-matrix

Samsdram

Os valores da matriz D estão diminuindo à medida que o número de iterações aumenta. Talvez este seja o problema como você apontou.

Andy Amos

@ Andy Amos: Eu apostaria dinheiro nisso. Como o @whuber aponta, é quase impossível que tenha autovalores negativos se você o computar diretamente, e os zeros (por serem deficientes em classificação) devem ser resolvidos adicionando desde a diagonal positiva - a menos que alguns desses elementos são realmente pequenos. Tente gerar alguns dados a partir de um modelo em que seja bem grande e . Quanto mais dados, melhor para que as estimativas sejam precisas e estáveis. Isso dirá pelo menos se há algum problema na sua implementação.

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