Cómo evitar que los vectores propios se diferencien de una computadora a otra

Siguiendo esta pregunta sobre cómo encontrar el estado estable de Markov, ahora me encuentro con el problema de que funciona perfectamente en mi computadora de laboratorio, pero no funciona en ninguna otra computadora. Específicamente, siempre encuentra el número correcto de valores propios cercanos a uno y, por lo tanto, qué nodos son nodos atractivos, pero no los encuentra de forma consistente y no están agrupados correctamente. Por ejemplo, al usar la matriz de transición de 64×64 a continuación, las computadoras en las que no funciona siempre producen una de tres colecciones diferentes de atractores al azar. En la matriz más pequeña M1 a continuación, todas las computadoras probadas obtienen el mismo resultado correcto para los grupos de atractores y la distribución estacionaria.

Todas las máquinas probadas ejecutan Win7x64 y WinPython-64bit-2.7.9.4. Una computadora siempre lo hace bien, otras tres siempre lo hacen de la misma manera. Según varias publicaciones que encontré como esto y esto , esto parece que podría deberse a diferencias en la precisión de punto flotante de los cálculos. Desafortunadamente no sé cómo arreglar esto; Quiero decir, no sé cómo modificar el código para extraer los valores propios de la matriz para forzar una precisión específica que todas las computadoras pueden manejar (y no creo que tenga que ser muy precisa para este propósito) .

Esa es mi mejor estimación actual de cómo podrían diferir los resultados. Si tiene una mejor idea de por qué sucede esto y cómo evitar que esto suceda, entonces eso también es genial.

Si hay una manera de hacer que scipy sea coherente de una ejecución a otra y de una computadora a otra, no creo que dependa de los detalles de mi método, sino porque se solicitó, aquí está. Ambas matrices tienen 3 atractores. En M1, el primero [1,2] es una órbita de dos estados, los otros dos [7] y [8] son ​​equilibrios. M2 es una matriz de transición de 64×64 con equilibrios en [2] y [26], así como una órbita que usa [7,8].

Pero en lugar de encontrar ese conjunto de atractores, a veces informa [[26], [2], [26]] y otras veces [[2,7,8,26], [2], [26]] y algunas veces … . no recibe la misma respuesta en cada ejecución, y nunca recibe [[2], [7,8], [26]] (en ningún orden).

import numpy as np import scipy.linalg M1 = np.array([[0.2, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.6, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.2, 0.0, 0.1, 0.1, 0.3, 0.3], [0.0, 0.0, 0.2, 0.2, 0.2, 0.2, 0.1, 0.1], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]]) M2 = np.genfromtxt('transitionMatrix1.csv', delimiter=',') # For easy switching M = M2 # Confirm the matrix is a valid Markov transition matrix #print np.sum(M,axis=1) 

El rest es el mismo código de la pregunta anterior, incluido aquí para su conveniencia.

 #create a list of the left eigenvalues and a separate array of the left eigenvectors theEigenvalues, leftEigenvectors = scipy.linalg.eig(M, right=False, left=True) # for stationary distribution the eigenvalues and vectors are always real, and this speeds it up a bit theEigenvalues = theEigenvalues.real #print theEigenvalues leftEigenvectors = leftEigenvectors.real #print leftEigenvectors # set how close to zero is acceptable as being zero...1e-15 was too low to find one of the actual eigenvalues tolerance = 1e-10 # create a filter to collect the eigenvalues that are near enough to zero mask = abs(theEigenvalues - 1) < tolerance # apply that filter theEigenvalues = theEigenvalues[mask] # filter out the eigenvectors with non-zero eigenvalues leftEigenvectors = leftEigenvectors[:, mask] # convert all the tiny and negative values to zero to isolate the actual stationary distributions leftEigenvectors[leftEigenvectors < tolerance] = 0 # normalize each distribution by the sum of the eigenvector columns attractorDistributions = leftEigenvectors / leftEigenvectors.sum(axis=0, keepdims=True) # this checks that the vectors are actually the left eigenvectors attractorDistributions = np.dot(attractorDistributions.T, M).T # convert the column vectors into row vectors (lists) for each attractor attractorDistributions = attractorDistributions.T print attractorDistributions # a list of the states in any attractor with the stationary distribution within THAT attractor #theSteadyStates = np.sum(attractorDistributions, axis=1) #print theSteadyStates 

La respuesta desafortunada es que no hay manera de arreglar la semilla para scipy, y por lo tanto no hay forma de forzarla a generar valores consistentes. Esto también significa que no hay forma de que pueda producir respuestas correctas de manera confiable porque solo una respuesta es correcta. Mis bashs de obtener una respuesta definitiva o una solución de la gente escéptica fueron completamente descartados , pero alguien puede encontrar algo de sabiduría en esas palabras al enfrentar este problema.

Como ejemplo concreto del problema, cuando ejecuta el código anterior, a veces puede obtener el siguiente conjunto de vectores propios que supuestamente representan los estados estacionarios de cada uno de los atractores del sistema. La computadora de mi hogar siempre produce este resultado (que es diferente de mi computadora portátil y de laboratorio). Como se indica en la pregunta, los atractores correctos son [[2],[7,8],[26]] . Los equilibrios de [2] y [6] se identifican correctamente, pero la distribución para [7,8] cambio devuelve una distribución de probabilidad no válida sobre [2,26] . La respuesta correcta es [0.19835, 0.80164] sobre [7,8] respectivamente. Mi computadora de laboratorio encuentra esa solución correctamente, pero hasta ahora otras seis computadoras no lo han logrado.

Lo que esto significa es que (a menos que haya algún otro error no identificado en mi código) scipy.linalg no sirve para encontrar estados estables de modelos de Markov. Aunque funciona parte del tiempo, no se puede confiar en que proporcione la respuesta correcta y, por lo tanto, debe evitarse por completo … al menos para los estados estacionarios del modelo de Markov y probablemente para todo lo que tenga que ver con los vectores propios. Simplemente no funciona.

Publicaré un código sobre cómo generar de manera confiable la distribución estacionaria de un modelo de Markov sin usar scipy si alguien hace una pregunta al respecto. Funciona un poco más lento, pero siempre es igual y siempre es correcto.

 [[ 0. 0. 0. ] [ 0. 0. 0. ] [ 0.25707958 1. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0.06867772 0. 1. ] [ 0. 0. 0. ] [ 0. 0. 0. ] [ 0. 0. 0. ] ... [ 0. 0. 0. ]]