scipy.linalg.eig return autovalores complejos para la matriz de covarianza?

Question

Los valores propios de una matriz de covarianza deben ser reales y no negativos porque las matrices de covarianza son simétricas y semi positivas definidas.

Sin embargo, echar un vistazo al siguiente experimento con scipy:

>>> a=np.random.random(5) 
>>> b=np.random.random(5) 
>>> ab = np.vstack((a,b)).T 
>>> C=np.cov(ab) 
>>> eig(C) 
7.90174997e-01 +0.00000000e+00j, 
2.38344473e-17 +6.15983679e-17j, 
2.38344473e-17 -6.15983679e-17j, 
-1.76100435e-17 +0.00000000e+00j, 
5.42658040e-33 +0.00000000e+00j 

Sin embargo, reproduciendo el ejemplo anterior en Matlab funciona correctamente:

a = [0.6271, 0.4314, 0.3453, 0.8073, 0.9739] 
b = [0.1924, 0.3680, 0.0568, 0.1831, 0.0176] 
C=cov([a;b]) 
eig(C) 
-0.0000 
-0.0000 
0.0000 
0.0000 
0.7902 

Answer 1

Se han planteado dos cuestiones:

1. Los valores propios devueltos por scipy.linalg.eig no son reales.
2. Algunos de los autovalores son negativos.

Ambos problemas son el resultado de errores introducidos por errores de truncamiento y redondeo, que siempre ocurren con algoritmos iterativos que usan aritmética de coma flotante. Tenga en cuenta que los resultados de Matlab también produjeron valores propios negativos.

Ahora, para un aspecto más interesante de la cuestión: ¿por qué es real el resultado de Matlab, mientras que el resultado de SciPy tiene algunos componentes complejos?

eig de Matlab detecta si la matriz de entrada es real simétrica o hermitiana y utiliza la factorización de Cholesky cuando lo es. Consulte la descripción del argumento chol en el eig documentation. Esto no se hace automáticamente en SciPy.

Si desea utilizar un algoritmo que explote la estructura de una matriz real simétrica o hermitiana, use scipy.linalg.eigh. Para el ejemplo de la pregunta:

>>> eigh(C, eigvals_only=True) 
array([ -3.73825923e-17, -1.60154836e-17, 8.11704449e-19, 
     3.65055777e-17, 7.90175615e-01]) 

Este resultado es el mismo que el de Matlab, si se redondearán con el mismo número de dígitos de precisión que Matlab impresa.

Answer 2

Lo que estás experimentando es la inestabilidad numérica debido a las limitaciones en la precisión del punto flotante.

Tenga en cuenta que:

(1) MATLAB también arrojaron valores negativos, pero el formato de impresión se establece en short y que no ven la precisión completa de la doble almacenado en la memoria. Use format long g para imprimir más decimales

(2) Todas las partes imaginarias devueltas por numpy's linalg.eig están cerca de la precisión de la máquina. Por lo tanto, deberías considerarlos cero.