# Eigenvectors of transition matrix is NaN

Hello,

I am trying to fit a Hidden Markov Model in CmdStan.
I specified a transition matrix `Gamma` that depends on the parameters of the model, and it is such that `Gamma[j, ]` is simplex for each state `j`.
Since I assume that my system is at steady state, I would like to define the initial state vector `rho` as a stationary probability vector, hence an eigenvector of `Gamma` with eigenvalue 1.
In the transformed parameter block I wrote the following:

``````complex_vector[n_states] eigvals = eigenvalues(Gamma);

complex_matrix[n_states, n_states] eigvecs = eigenvectors(Gamma);

simplex[n_states] rho;
for (v in 1 : n_states){
if (eigvals[v] == to_complex(1, 0)){
rho = get_real(eigvecs[1 : n_states, v]); // column vector
rho = rho / sum(rho); // normalize
}
}
``````

I do find a eigenvector with eigenvalue 1, but I repeatedly get the following error for `rho`:

``````Exception: test_new_func_v4_model_namespace::log_prob: rho is not a valid simplex. sum(rho) = nan, but should be 1 (in 'test_new_func_v4.stan', line 195, column 4 to column 26)
``````

and if I print the values of `rho` is indeed all `nan` entries.

// ------------------------------- //

In particular, what I found is that the evaluation `eigvals[v] == to_complex(1, 0)` is always failing, even if `eigvals[v]` appears to have real part 1 and imaginary part 0. If I print the output of

``````for (i in 1 : n_max){
print(get_real(eigvals[i]));
print(get_imag(eigvals[i]));
print(eigvals[i] == to_complex(1.0, 0.0));
print( get_real(eigvals[i]) == 1.0 );
print( get_imag(eigvals[i]) == 0.0 );
}
``````

I get for `i=0`

``````1
0
0
0
1
``````

I am not sure what I am not understandingâ€¦

CmdStan v2.34.1, CmdStanPy 1.2.2

@pgree experienced a similar issue which we tracked down to instability in one of our dependencies, which was fixed but not yet released by them: Nans in complex eigendecomposition

In his case the issue seemed to be tied to very small but not quite zero values - maybe something similar is happening here, which is why the equality test is failing? (`print()` will do some rounding)

1 Like

Thank you for your reply!! Is there a way I could get around this?

This seems to â€śfixâ€ť it:

``````for (v in 1 : n_states){
real r = get_real(eigvals[v]);
real im = get_imag(eigvals[v]);
if (im == 0 && r < 1.0001 && r > 0.9999){
rho = get_real(eigvecs[1 : n_states, v]); // column vector
rho = rho / sum(rho); // normalize
}
}
``````

So youâ€™re probably correct @WardBrian

Assuming it is the same issue, weâ€™re waiting on the Stan side for the Eigen package to release version 3.4.1, which includes the direct fix. I donâ€™t know if there is a timeline for when they are planning on doing that.

If you wanted to be a bit experimental, you could replace the copy of Eigen inside stan with their branch with the fix

There are other things you could try to improve the numerics, like the snippet you just posted (or the Eigen folks recommended matrix balancing [1401.5766] On matrix balancing and eigenvector computation to me) but they will be model-specific

Thank you very much!