Eigenvectors of transition matrix is NaN

Modeling
1

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…

Thank you in advance for your help!

CmdStan v2.34.1, CmdStanPy 1.2.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
3

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

4

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

5

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

6

Thank you very much!