Comment rendre une matrice positive définie?

J'essaie d'implémenter un algorithme EM pour le modèle d'analyse factorielle suivant;

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$

où est un vecteur aléatoire à p dimensions, est un vecteur à q dimensions de variables latentes et est une matrice pxq de paramètres. $W_j$ $a_j$ $B$

En raison d'autres hypothèses utilisées pour le modèle, je sais que où est la matrice de covariance de variance des termes d'erreur , = diag ( , , ..., ). $W_j\sim N(\mu, BB'+D)$ $D$ $e_j$ $D$ $\sigma_1^2$ $\sigma_2^2$ $\sigma_p^2$

Pour l'algorithme EM au travail, je fais des itérations de dôme impliquant une estimation de et matrices et pendant ces itérations je calculer l'inverse de à chaque itération en utilisant de nouvelles estimations de et . Malheureusement au cours des itérations, perd son caractère définitif positif (mais cela ne devrait pas parce qu'il s'agit d'une matrice de variance-covariance) et cette situation ruine la convergence de l'algorithme. Mes questions sont: $B$ $D$ $BB'+D$ $B$ $D$ $BB'+D$

Cette situation montre-t-elle qu'il y a un problème avec mon algorithme, car la probabilité devrait augmenter à chaque étape de l'EM?
Quels sont les moyens pratiques pour rendre une matrice positive définie?

Edit: Je calcule l'inverse en utilisant un lemme d'inversion de matrice qui déclare 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}$

où le côté droit n'implique que les inverses des matrices . $q\times q$

factor-analysis expectation-maximization Andy Amos
la source

Cela pourrait aider à mieux comprendre comment

«perd» son caractère définitif positif. Cela implique que

(ou les deux) deviennent définis non positifs. C'est difficile à faire lorsque

est calculé directement à partir de

et encore plus difficile lorsque

est calculé comme une matrice diagonale avec des carrés sur sa diagonale!

B B^{'} + D

$BB'+D$

B B^{'}

$BB'$

D

$D$

B B^{'}

$BB'$

B

$B$

D

$D$

whuber

@whuber Typiquement dans FA

, donc

n'est jamais défini définitivement positif. Mais (théoriquement)

devrait l'être, en supposant que les

sont tous supérieurs à zéro.

q < p

$q<p$

B B^{'}

$BB'$

B B^{'} + D

$BB' + D$

σ_{j}^{2}

$\sigma^2_j$

JMS

Ceci est lié à cette question: stats.stackexchange.com/questions/6364/…

Gilead

@JMS Merci. Je pense que mon commentaire est toujours pertinent:

peut être indéfini, mais ne doit toujours pas avoir de valeurs propres négatives. Cependant, des problèmes surviendront lorsque le plus petit des

est comparable à une erreur numérique dans l'algorithme d'inversion. Si tel est le cas, une solution consiste à appliquer SVD à

et zéro les valeurs propres vraiment petites (ou négatives), puis recalcule

et ajouter

B B^{'}

$BB'$

σ_{i}^{2}

$\sigma_i^2$

B B^{'}

$BB'$

B B^{'}

$BB'$

D

$D$

whuber

Ce doit être de petits éléments en

;

devrait être bien conditionné sinon puisque

D

$D$

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

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

q < p

$q<p$

JMS

Réponses:

OK, puisque vous faites FA, je suppose que est de rang de colonne complet et . Nous avons besoin de quelques détails supplémentaires. Cela peut être un problème numérique; cela peut également être un problème avec vos données. $B$ $q$ $q<p$

Comment calculez-vous l'inverse? Avez-vous besoin de l'inverse explicitement ou pouvez-vous ré-exprimer le calcul comme la solution d'un système linéaire? (c.-à-d. pour obtenir résoudre pour x, qui est généralement plus rapide et plus stable) $A^{-1}b$ $Ax=b$

Qu'arrive-t-il à ? Les estimations sont-elles vraiment petites / 0 / négatives? Dans un certain sens, c'est le lien critique, car est bien sûr déficient en rang et définit une matrice de covariance singulière avant d'ajouter , vous ne pouvez donc pas l'inverser. L'ajout de la matrice diagonale positive le rend techniquement complet, mais pourrait encore être horriblement mal conditionné si est petit. $D$ $BB'$ $D$ $D$ $BB'+D$ $D$

Souvent, l'estimation des variances idiosyncratiques (votre , les éléments diagonaux de ) est proche de zéro ou même négative; ce sont les cas Heywood. Voir par exemple http://www.technion.ac.il/docs/sas/stat/chap26/sect21.htm (tout texte FA devrait également en discuter, c'est un problème très ancien et bien connu). Cela peut résulter de la spécification incorrecte du modèle, des valeurs aberrantes, de la malchance, des éruptions solaires ... le MLE est particulièrement sujet à ce problème, donc si votre algorithme EM est conçu pour obtenir le MLE. $\sigma^2_i$ $D$

Si votre algorithme EM approche d'un mode avec de telles estimations, il est possible que perde sa définition positive, je pense. Il existe différentes solutions; Personnellement, je préférerais une approche bayésienne, mais même alors, vous devez être prudent avec vos priors (des priors incorrects ou même des priors appropriés avec trop de masse près de 0 peuvent avoir le même problème pour essentiellement la même raison) $BB'+D$

JMS
la source

Permettez-moi de souligner que dans la partie principale des algorithmes, vous ne voulez jamais inverser réellement une matrice. Vous devrez peut-être à la toute fin pour obtenir les estimations standard. Voir cet article de blog johndcook.com/blog/2010/01/19/dont-invert-that-matrix

Samsdram

Les valeurs de la matrice D deviennent de plus en plus petites à mesure que le nombre d'itérations augmente. C'est peut-être le problème comme vous l'avez souligné.

Andy Amos

@Andy Amos: Je parierais de l'argent dessus. Comme @whuber le fait remarquer, il est presque impossible que

ait des valeurs propres négatives si vous le calculez directement, et les zéros (de défaut de classement) devraient être pris en compte en ajoutant

depuis sa diagonale positive - à moins que certains d'entre eux les éléments sont vraiment petits. Essayez de générer des données à partir d'un modèle où

sont assez grandes et

. Plus il y a de données, mieux c'est pour que les estimations soient précises et stables. Cela vous indiquera au moins s'il y a un problème dans votre implémentation.

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