¿Cómo hacer que una matriz sea definida positiva?

Question

Estoy intentando implementar un algoritmo EM para el siguiente modelo de análisis factorial;

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

donde $W_j$ es un vector aleatorio p-dimensional, $a_j$ es un vector q-dimensional de variables latentes y $B$ es una matriz pxq de parámetros.

Como resultado de otros supuestos utilizados para el modelo, sé que $W_j\sim N(\mu, BB'+D)$ donde $D$ es la matriz de covarianza varianza de los términos de error $e_j$ , $D$ = diag( $\sigma_1^2$ , $\sigma_2^2$ ,..., $\sigma_p^2$ ).

Para que el algoritmo EM funcione, estoy haciendo iteraciones de cúpula que implican la estimación de $B$ y $D$ matrices y durante estas iteraciones estoy calculando la inversa de $BB'+D$ en cada iteración utilizando nuevas estimaciones de $B$ y $D$ . Desafortunadamente, durante el curso de las iteraciones, $BB'+D$ pierde su definición positiva (pero no debería porque es una matriz de varianza-covarianza) y esta situación arruina la convergencia del algoritmo. Mis preguntas son:

1. ¿Demuestra esta situación que hay algo mal en mi algoritmo, ya que la probabilidad debería aumentar en cada paso de EM?

2. ¿Cuáles son las formas prácticas de hacer que una matriz sea definida positiva?

Edición: Estoy calculando la inversa utilizando un lema de inversión de matrices que establece que:

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

donde el lado derecho implica sólo los inversos de $q\times q$ matrices.

Answer 1

OK, ya que estás haciendo FA asumo que $B$ tiene rango de columna completo $q$ y $q<p$ . Aunque necesitamos algunos detalles más. Puede tratarse de un problema numérico; también puede ser un problema con tus datos.

¿Cómo se calcula la inversa? ¿Necesitas la inversa explícitamente, o puedes reexpresar el cálculo como la solución de un sistema lineal? (es decir, para obtener $A^{-1}b$ resolver $Ax=b$ para x, que suele ser más rápido y estable)

¿Qué le ocurre a $D$ ? ¿Son las estimaciones realmente pequeñas/0/negativas? En cierto sentido es el eslabón crítico, porque $BB'$ es, por supuesto, de rango deficiente y define una matriz de covarianza singular antes de añadir $D$ por lo que no se puede invertir. Sumando la matriz diagonal positiva $D$ técnicamente lo hace de rango completo pero $BB'+D$ todavía podría estar horriblemente mal acondicionado si $D$ es pequeño.

A menudo, la estimación de las varianzas idiosincrásicas (su $\sigma^2_i$ los elementos diagonales de $D$ ) es cercano a cero o incluso negativo; estos casos se denominan casos de Heywood. Véase, por ejemplo http://www.technion.ac.il/docs/sas/stat/chap26/sect21.htm (cualquier texto de FA debería hablar de esto también, es un problema muy antiguo y conocido). Esto puede ser el resultado de una mala especificación del modelo, valores atípicos, mala suerte, erupciones solares... el MLE es particularmente propenso a este problema, así que si su algoritmo EM está diseñado para obtener el MLE, tenga cuidado.

Si su algoritmo EM se aproxima a un modo con tales estimaciones es posible que $BB'+D$ perder su definición positiva, creo. Hay varias soluciones; personalmente preferiría un enfoque bayesiano, pero incluso en ese caso hay que tener cuidado con las variables a priori (las variables a priori inadecuadas o incluso las adecuadas con demasiada masa cerca de 0 pueden tener el mismo problema básicamente por la misma razón).