Como tornar uma matriz positiva definitiva?

10

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

Wj=μ+Baj+ejforj=1,,n

onde é um vetor aleatório p-dimensional, é um vetor q-dimensional de variáveis ​​latentes e é uma matriz de parâmetros pxq.a j BWjajB

Como resultado de outras suposições usadas para o modelo, eu sei que que D é a matriz de covariância de variância dos termos de erro e_j , D = diag ( \ sigma_1 ^ 2 , \ sigma_2 ^ 2 , ..., \ sigma_p ^ 2 ).DWjN(μ,BB+D)DejDσ12σ22σp2

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:

  1. Essa situação mostra que há algo errado com meu algoritmo, pois a probabilidade deve aumentar a cada etapa do EM?

  2. 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:

(BB+D)1=D1D1B(Iq+BD1B)1BD1

onde o lado direito envolve apenas o inverso de matrizes q×q .

Andy Amos
fonte
11
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! BB+DBBDBBBD
whuber
@whuber Normalmente em FA , então nem sempre é positivo. Mas (teoricamente) deveria ser, assumindo que os são todos maiores que zero. B B B B + D σ 2 jq<pBBBB+Dσj2
JMS
Isso está relacionado a esta pergunta: stats.stackexchange.com/questions/6364/…
Gilead
11
@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 . σ 2 i B B B B DBBσi2BBBBD
whuber
11
Tem que haver pequenos elementos em ; deve estar bem condicionado, caso contrário, uma vez queI q + B D - 1 B q < pDIq+BD1Bq<p
JMS

Respostas:

3

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.q q < pBqq<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 x = bA1bAx=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.B B D D B B + D DDBBDDBB+DD

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. Dσi2D

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
11
@ 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. D σ 2 iq B 2 i qσ 2 iBBDσi2qBiq2σi2
JMS