¿Cómo hacer una matriz positiva definida?

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

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$

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

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

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

¿Esta situación muestra que hay algo mal con mi algoritmo ya que la probabilidad debería aumentar en cada paso de EM?
¿Cuáles son las formas prácticas de hacer una matriz positiva definida?

Editar: estoy calculando el inverso usando un lema de inversión de matriz que establece 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}$

donde el lado derecho involucra solo las inversas de matrices. $q\times q$

factor-analysis expectation-maximization Andy Amos
fuente

Podría ayudar a comprender mejor cómo

"pierde" su definición positiva. Esto implica que

(o ambos) se están convirtiendo en definitivos no positivos. ¡Eso es difícil de hacer cuando

se calcula directamente desde

y aún más difícil cuando

se calcula como una matriz diagonal con cuadrados en su diagonal!

B B^{'} + D

$BB'+D$

B B^{'}

$BB'$

D

$D$

B B^{'}

$BB'$

B

$B$

D

$D$

whuber

@whuber Normalmente en FA

, por lo que

no es nunca definitivo positivo. Pero (en teoría)

debería ser, suponiendo que las

's sean mayores que cero.

q < p

$q<p$

B B^{'}

$BB'$

B B^{'} + D

$BB' + D$

σ_{j}^{2}

$\sigma^2_j$

JMS

Esto está relacionado con esta pregunta: stats.stackexchange.com/questions/6364/…

Gilead

@JMS Gracias. Creo que mi comentario sigue siendo pertinente:

puede ser indefinido, pero aún así no debe tener ningún valor propio negativo. Sin embargo, surgirán problemas cuando el más pequeño de

sea comparable al error numérico en el algoritmo de inversión. Si este es el caso, una solución es aplicar a SVD

y de cero a cabo las muy pequeñas (o negativas) valores propios, a continuación, volver a calcular

y añadir

B B^{'}

$BB'$

σ_{i}^{2}

$\sigma_i^2$

B B^{'}

$BB'$

B B^{'}

$BB'$

D

$D$

whuber

Tiene que ser elementos pequeños en

;

debería estar bien acondicionado de lo contrario ya que

D

$D$

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

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

q < p

$q<p$

JMS

Respuestas:

OK, ya que estás haciendo FA, supongo que es de rango de columna completo y . Sin embargo, necesitamos algunos detalles más. Esto puede ser un problema numérico; También puede ser un problema con sus datos. $B$ $q$ $q<p$

¿Cómo estás calculando el inverso? ¿Necesita el inverso explícitamente o puede volver a expresar el cálculo como la solución a un sistema lineal? (es decir, para obtener resuelva para x, que generalmente es más rápido y más estable) $A^{-1}b$ $Ax=b$

¿Qué le está pasando a ? ¿Son las estimaciones realmente pequeñas / 0 / negativas? En cierto sentido, es el enlace crítico, porque es, por supuesto, deficiente en el rango y define una matriz de covarianza singular antes de agregar , por lo que no puede invertirla. Agregar la matriz diagonal positiva técnicamente lo hace rango completo, pero aún podría estar terriblemente mal condicionado si es pequeño. $D$ $BB'$ $D$ $D$ $BB'+D$ $D$

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

Si su algoritmo EM se acerca a un modo con tales estimaciones, es posible que pierda su definición positiva, creo. Hay varias soluciones; Personalmente, preferiría un enfoque bayesiano, pero incluso entonces debe tener cuidado con sus antecedentes (los antecedentes impropios o incluso los anteriores apropiados con demasiada masa cerca de 0 pueden tener el mismo problema básicamente por la misma razón) $BB'+D$

JMS
fuente

Permítanme decir que en la parte principal de los algoritmos, nunca querrán invertir realmente una matriz. Sin embargo, es posible que necesite al final para obtener las estimaciones estándar. Ver esta publicación de blog johndcook.com/blog/2010/01/19/dont-invert-that-matrix

Samsdram

Los valores de la matriz D se hacen cada vez más pequeños a medida que aumenta el número de iteraciones. Tal vez este sea el problema como usted señaló.

Andy Amos

@Andy Amos: Apostaría dinero por ello. Como @whuber señala que es casi imposible que

tenga valores propios negativos si lo está calculando directamente, y los ceros (por ser deficientes en el rango) deben ser atendidos agregando

desde su diagonal positiva, a menos que algunos de esos Los elementos son realmente pequeños. Intente generar algunos datos de un modelo donde

es bastante grande y

. Cuantos más datos, mejor para que las estimaciones sean precisas y estables. Eso al menos le dirá si hay un problema en su implementación.

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