Nans in complex eigendecomposition

Developers
1

I have a Stan model that uses an eigendecomposition of a complex matrix in the generated quantities block. I sometimes get nans in the orthogonal matrix of the eigendecomposition even though the original matrix has no nans. I’m using cmdstanpy: 1.2.0 and CmdStan 2.33.

The Stan program is below. The last ~20 lines of the generated quantities block prints out the number of nans in the original matrix and the orthogonal matrix outputted by the eigendecomposition.

functions {
  real khat(real xi, real alpha, real l) {
    real f = alpha * sqrt(2 * pi() * pow(l, 2)) * exp( (-2 * pow(pi(), 2) * pow(l, 2)) * pow(xi, 2));
    return f;
  }
}
data {
  int<lower=1> m;
  int<lower=1> N;
  array[N] real x;
  vector[N] y;
  complex_matrix[m, m] ftf;
  complex_vector[m] fty;
  vector[m] xis;
  complex_matrix[N, m] fmat;
}
transformed data {
  vector[N] mu = rep_vector(0, N);
  real alpha = 1;
  real sigma = 1;
}
parameters {
  real<lower=0> l;
}
model {
  real sq_sigma = square(sigma);

  matrix[N, N] L_K;
  matrix[N, N] K = gp_exp_quad_cov(x, alpha, l);

  // diagonal elements
  for (i in 1:N) {
    K[i, i] = K[i, i] + sq_sigma;
  }
  L_K = cholesky_decompose(K);

  l ~ normal(0, 1);
  alpha ~ std_normal();
  sigma ~ std_normal();

  y ~ multi_normal_cholesky(mu, L_K);
}
generated quantities {
  // diagonal matrix
  complex_vector[m] ds;
  for (i in 1:m) {
    ds[i] =  sqrt((xis[m] - xis[1]) / (m - 1) * khat(xis[i], alpha, l));
  }
  // construct matrix D * F*F * D
  complex_matrix[N, m] phi = diag_post_multiply(fmat, ds);
  complex_matrix[m, m] kmat_m = add_diag(conj(phi') * phi, rep_vector(pow(sigma, 2), m));

  // eigendecomposition
  tuple(complex_matrix[m, m], complex_vector[m]) ud;
  ud = eigendecompose(kmat_m);
  complex_matrix[m, m] u = ud.1;
  complex_vector[m] d = ud.2;

  // take inverse of diagonal matrix
  complex_vector[m] d_inv;
  real logdet = 0;
  for (i in 1:m) {
    d_inv[i] = 1 / d[i];
    logdet += log(get_real(d[i]));
  }

  // rand hand side, D * F^t y
  complex_vector[m] phi_y = ds .* fty;

  // solve linear system
  complex_vector[m] beta = u * diag_pre_multiply(d_inv, conj(u')) * phi_y;

  // construct phi for evaluation of posterior mean at data
  complex_vector[N] pos_mean = phi * beta;

  // check for nans in kmat_m
  int<lower=0> nnan = 0;
  for (i in 1:m) {
    for (j in 1:m) {
       if (is_nan(get_real(kmat_m[i, j])) || is_nan(get_imag(kmat_m[i, j])))
         nnan += 1;
    }
  }
  print("kmat nans:");
  print(nnan);

  // check for nans in u
  nnan = 0;
  for (i in 1:m) {
    for (j in 1:m) {
       if (is_nan(get_real(u[i, j])) || is_nan(get_imag(u[i, j])))
         nnan += 1;
    }
  }
  print("u nans:");
  print(nnan);

  if (nnan > 0)
    print(kmat_m);
}

Here’s the Python code:

import cmdstanpy as stan
import numpy as np

np.random.seed(1)

# data
n = 20
xs = np.random.rand(n)
xs.sort()
sig = 1.0
ys = np.sin(2 * np.pi * xs) + sig * np.random.randn(n)
m = 21
xis = np.linspace(-6, 6, num=m)

# precompute
fmat = np.outer(xs, xis)
fmat = np.exp(2 * np.pi * 1j * fmat)
fmat = np.matrix(fmat)
ftf = np.dot(fmat.H, fmat)
fty = np.dot(fmat.H, ys)
fty = np.ravel(fty)

# sample
data = {
    "m": np.size(xis),
    "N": np.size(ys),
    "x": xs,
    "y": ys,
    "ftf": ftf,
    "fty": fty,
    "xis": xis,
    "fmat": fmat
}
stan_file = 'program.stan'
inits = {'l': 0.5}
my_model = stan.CmdStanModel(stan_file=stan_file)
fit = my_model.sample(data=data, iter_sampling=50, iter_warmup=50, show_console=True, inits=inits)

# extract samples
samples = fit.stan_variables()['pos_mean'].astype(float)
pos_means = samples.mean(axis=0)
print(f'Nans?: {np.sum(np.isnan(pos_means)) > 0}')

Thanks in advance for any help.

2

I was able to re-create this error. It looks like it is coming from inside Eigen, the library we use for our matrix types and most of our linear algebra.

Here is a backtrace on my system of where the NANs are coming from :

Summary 
Program received signal SIGFPE, Arithmetic exception.
Eigen::internal::unary_evaluator<Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs2_op<std::complex<double> >, Eigen::Block<Eigen::Block<Eigen::Block<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false> const, -1, 1, false> const>, Eigen::internal::IndexBased, double>::coeff (col=0, row=0, this=0x7fffffffb660) at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/CoreEvaluators.h:583
583         return m_d.func()(m_d.argImpl.coeff(row, col));
(gdb) bt
#0  Eigen::internal::unary_evaluator<Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs2_op<std::complex<double> >, Eigen::Block<Eigen::Block<Eigen::Block<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false> const, -1, 1, false> const>, Eigen::internal::IndexBased, double>::coeff (col=0, row=0, this=0x7fffffffb660) at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/CoreEvaluators.h:583
#1  Eigen::internal::redux_evaluator<Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs2_op<std::complex<double> >, Eigen::Block<Eigen::Block<Eigen::Block<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false> const, -1, 1, false> const> >::coeffByOuterInner (inner=0, 
    outer=0, this=0x7fffffffb660) at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/Redux.h:381
#2  Eigen::internal::redux_impl<Eigen::internal::scalar_sum_op<double, double>, Eigen::internal::redux_evaluator<Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs2_op<std::complex<double> >, Eigen::Block<Eigen::Block<Eigen::Block<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false> const, -1, 1, false> const> >, 0, 0>::run<Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs2_op<std::complex<double> >, Eigen::Block<Eigen::Block<Eigen::Block<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false> const, -1, 1, false> const> > (xpr=..., func=..., 
    eval=...) at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/Redux.h:202
#3  Eigen::DenseBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs2_op<std::complex<double> >, Eigen::Block<Eigen::Block<Eigen::Block<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false> const, -1, 1, false> const> >::redux<Eigen::internal::scalar_sum_op<double, double> > (func=..., this=0x7fffffffb740) at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/Redux.h:418
#4  Eigen::DenseBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs2_op<std::complex<double> >, Eigen::Block<Eigen::Block<Eigen::Block<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false> const, -1, 1, false> const> >::sum (this=0x7fffffffb740)
    at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/Redux.h:463
#5  Eigen::MatrixBase<Eigen::Block<Eigen::Block<Eigen::Block<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false> const, -1, 1, false> >::squaredNorm (this=0x7fffffffb6a0) at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/Dot.h:98
#6  Eigen::MatrixBase<Eigen::Block<Eigen::Block<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false> >::makeHouseholder<Eigen::VectorBlock<Eigen::Block<Eigen::Block<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false>, -1> > (
    this=this@entry=0x7fffffffb9c0, essential=..., tau=..., beta=@0x7fffffffb910: -1.2914366647217514e-41)
    at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Householder/Householder.h:78
#7  0x00005555555b67f7 in Eigen::MatrixBase<Eigen::Block<Eigen::Block<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>, -1, 1, true>, -1, 1, false> >::makeHouseholderInPlace (this=this@entry=0x7fffffffb9c0, tau=..., beta=@0x7fffffffb910: -1.2914366647217514e-41)
    at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Householder/Householder.h:46
#8  0x00005555555e808a in Eigen::HessenbergDecomposition<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::_compute (
    matA=Eigen::Matrix<std::complex<double>,21,21,ColMajor> (data ptr: 0x555555896e50) = {...}, 
    hCoeffs=Eigen::Matrix<std::complex<double>,20,1,ColMajor> (data ptr: 0x555555884d40) = {...}, 
    temp=Eigen::Matrix<std::complex<double>,1,21,RowMajor> (data ptr: 0x555555884be0) = {...})
    at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/MapBase.h:94
#9  0x00005555555f8f68 in Eigen::HessenbergDecomposition<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::compute<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> > (matrix=..., this=0x7fffffffbd88)
    at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Eigenvalues/./HessenbergDecomposition.h:161
#10 Eigen::internal::complex_schur_reduce_to_hessenberg<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>, true>::run (_this=..., 
    matrix=Eigen::Matrix<std::complex<double>,21,21,ColMajor> (data ptr: 0x555555891b70) = {...}, computeU=computeU@entry=true)
    at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Eigenvalues/ComplexSchur.h:361
#11 0x00005555555f970d in Eigen::ComplexSchur<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::compute<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> > (this=this@entry=0x7fffffffbd58, matrix=..., computeU=computeU@entry=true)
    at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Eigenvalues/ComplexSchur.h:337
#12 0x00005555555f9996 in Eigen::ComplexEigenSolver<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::compute<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> > (this=this@entry=0x7fffffffbd30, matrix=..., computeEigenvectors=computeEigenvectors@entry=true)
    at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Eigenvalues/ComplexEigenSolver.h:270
#13 0x00005555555fa06d in Eigen::ComplexEigenSolver<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::ComplexEigenSolver<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> > (this=0x7fffffffbd30, matrix=..., computeEigenvectors=<optimized out>)
    at stan/lib/stan_math/lib/eigen_3.4.0/Eigen/src/Core/EigenBase.h:49
#14 0x00005555555fa6fd in stan::math::eigendecompose<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>, (void*)0> (
    m=Eigen::Matrix<std::complex<double>,21,21,ColMajor> (data ptr: 0x555555891b70) = {...})
    at stan/lib/stan_math/stan/math/prim/fun/eigendecompose.hpp:66
#15 0x00005555555fb775 in program_model_namespace::program_model::write_array_impl<boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014u, 0u, 2147483563u>, boost::random::linear_congruential_engine<unsigned int, 40692u, 0u, 2147483399u> >, std::vector<double, std::allocator<double> >, std::vector<int, std::allocator<int> >, std::vector<double, std::allocator<double> >, (void*)0, (void*)0, (void*)0> (
    this=this@entry=0x55555587e470, base_rng__=..., params_r__=std::vector of length 1, capacity 1 = {...}, 
    params_i__=std::vector of length 0, capacity 0, vars__=std::vector of length 3781, capacity 3781 = {...}, 
    emit_transformed_parameters__=emit_transformed_parameters__@entry=true, emit_generated_quantities__=true, pstream__=0x7fffffffc350)
    at ../../py/cmdstanpy/program.hpp:672
#16 0x00005555555fd26e in program_model_namespace::program_model::write_array<boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014u, 0u, 2147483563u>, boost::random::linear_congruential_engine<unsigned int, 40692u, 0u, 2147483399u> > > (
    pstream=0x7fffffffc350, emit_generated_quantities=true, emit_transformed_parameters=<optimized out>, 
    vars=std::vector of length 3781, capacity 3781 = {...}, params_i=std::vector of length 0, capacity 0, 
    params_r=std::vector of length 1, capacity 1 = {...}, base_rng=..., this=0x55555587e470) at /usr/include/c++/11/ext/new_allocator.h:89
#17 stan::model::model_base_crtp<program_model_namespace::program_model>::write_array (this=0x55555587e470, rng=..., 
    theta=std::vector of length 1, capacity 1 = {...}, theta_i=std::vector of length 0, capacity 0, 
    vars=std::vector of length 3781, capacity 3781 = {...}, include_tparams=<optimized out>, include_gqs=true, msgs=0x7fffffffc350)
    at stan/src/stan/model/model_base_crtp.hpp:209
#18 0x00005555556801bf in stan::services::util::mcmc_writer::write_sample_params<stan::model::model_base, boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014u, 0u, 2147483563u>, boost::random::linear_congruential_engine<unsigned int, 40692u, 0u, 2147483399u> > > (this=this@entry=0x7fffffffc880, rng=..., sample=..., sampler=warning: RTTI symbol not found for class 'stan::mcmc::adapt_diag_e_nuts<stan::model::model_base, boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014u, 0u, 2147483563u>, boost::random::linear_congruential_engine<unsigned int, 40692u, 0u, 2147483399u> > >'
..., model=...) at /usr/include/c++/11/ext/new_allocator.h:89
#19 0x0000555555680c91 in stan::services::util::generate_transitions<stan::model::model_base, boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014u, 0u, 2147483563u>, boost::random::linear_congruential_engine<unsigned int, 40692u, 0u, 2147483399u> > > (
    sampler=warning: RTTI symbol not found for class 'stan::mcmc::adapt_diag_e_nuts<stan::model::model_base, boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014u, 0u, 2147483563u>, boost::random::linear_congruential_engine<unsigned int, 40692u, 0u, 2147483399u> > >'
..., num_iterations=num_iterations@entry=100, start=start@entry=100, finish=finish@entry=200, num_thin=num_thin@entry=1, 
    refresh=refresh@entry=100, save=true, warmup=false, mcmc_writer=..., init_s=..., model=..., base_rng=..., callback=..., logger=..., chain_id=1, 
    num_chains=1) at stan/src/stan/services/util/generate_transitions.hpp:73
#20 0x00005555556943e9 in stan::services::util::run_adaptive_sampler<stan::mcmc::adapt_diag_e_nuts<stan::model::model_base, boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014u, 0u, 2147483563u>, boost::random::linear_congruential_engine<unsigned int, 40692u, 0u, 2147483399u> > >, stan::model::model_base, boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014u, 0u, 2147483563u>, boost::random::linear_congruential_engine<unsigned int, 40692u, 0u, 2147483399u> > > (sampler=warning: RTTI symbol not found for class 'stan::mcmc::adapt_diag_e_nuts<stan::model::model_base, boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014u, 0u, 2147483563u>, boost::random::linear_congruential_engine<unsigned int, 40692u, 0u, 2147483399u> > >'
..., model=..., 
    cont_vector=std::vector of length 1, capacity 1 = {...}, num_warmup=num_warmup@entry=100, num_samples=num_samples@entry=100, 
    num_thin=num_thin@entry=1, refresh=100, save_warmup=false, rng=..., interrupt=..., logger=..., sample_writer=..., diagnostic_writer=..., 
    metric_writer=..., chain_id=1, num_chains=1) at stan/src/stan/services/util/run_adaptive_sampler.hpp:93
#21 0x000055555569492f in stan::services::sample::hmc_nuts_diag_e_adapt<stan::model::model_base> (model=..., init=..., init_inv_metric=..., 
    random_seed=random_seed@entry=3522013070, chain=chain@entry=1, init_radius=init_radius@entry=2, num_warmup=num_warmup@entry=100, num_samples=100, 
    num_thin=1, save_warmup=false, refresh=100, stepsize=stepsize@entry=1, stepsize_jitter=stepsize_jitter@entry=0, max_depth=10, 
    delta=delta@entry=0.80000000000000004, gamma=gamma@entry=0.050000000000000003, kappa=kappa@entry=0.75, t0=t0@entry=10, init_buffer=75, 
    term_buffer=50, window=25, interrupt=..., logger=..., init_writer=..., sample_writer=..., diagnostic_writer=..., metric_writer=...)
    at stan/src/stan/services/sample/hmc_nuts_diag_e_adapt.hpp:102
#22 0x00005555556965fa in stan::services::sample::hmc_nuts_diag_e_adapt<stan::model::model_base, std::shared_ptr<stan::io::var_context>, stan::callbacks::writer, stan::callbacks::unique_stream_writer<std::basic_ofstream<char, std::char_traits<char> >, std::default_delete<std::basic_ofstream<char, std::char_traits<char> > > >, stan::callbacks::unique_stream_writer<std::basic_ofstream<char, std::char_traits<char> >, std::default_delete<std::basic_ofstream<char, std::char_traits<char> > > >, stan::callbacks::json_writer<std::basic_ofstream<char, std::char_traits<char> >, std::default_delete<std::basic_ofstream<char, std::char_traits<char> > > > > (model=..., num_chains=num_chains@entry=1, init=std::vector of length 1, capacity 1 = {...}, 
    random_seed=random_seed@entry=3522013070, init_chain_id=init_chain_id@entry=1, init_radius=init_radius@entry=2, num_warmup=num_warmup@entry=100, 
    num_samples=100, num_thin=1, save_warmup=false, refresh=100, stepsize=stepsize@entry=1, stepsize_jitter=stepsize_jitter@entry=0, max_depth=10, 
    delta=delta@entry=0.80000000000000004, gamma=gamma@entry=0.050000000000000003, kappa=kappa@entry=0.75, t0=t0@entry=10, init_buffer=75, 
    term_buffer=50, window=25, interrupt=..., logger=..., init_writer=std::vector of length 1, capacity 1 = {...}, 
    sample_writer=std::vector of length 1, capacity 1 = {...}, diagnostic_writer=std::vector of length 1, capacity 1 = {...}, 
    metric_writer=std::vector of length 1, capacity 1 = {...}) at stan/src/stan/services/sample/hmc_nuts_diag_e_adapt.hpp:553
#23 0x000055555560c99b in cmdstan::command (argc=<optimized out>, argv=<optimized out>) at src/cmdstan/command.hpp:654
#24 0x000055555560e68c in main (argc=<optimized out>, argv=<optimized out>) at src/cmdstan/main.cpp:6

@stevebronder @spinkney - is this an Eigen bug? It doesn’t look like it is using any of our autodiff types, just doubles

3

I’m having a hard time replicating this outside of Stan, but just grabbing one iteration of k_mat it looks like you have many very small values (like 1e-200 small). I’ll look into this more but my guess is that dividing by any of these small numbers could produce +/- infinity values and then doing an operation on a combination of +/- infinity causes a NaN. Or multiplying them together would truncate to 0 and then a division would also return a NaN

           (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)             (1,0)  (4.10228e-278,0)  (5.13064e-250,0)  (4.99368e-225,0)  (5.81158e-204,0)  (1.76197e-186,0) (-5.50726e-173,0) (-3.26903e-163,0)  (7.06947e-158,0)  (1.24648e-155,0)  (2.71149e-157,0) (-2.59112e-163,0)  (1.30494e-172,0)  (2.46975e-186,0) (-1.91872e-203,0) (-7.73057e-226,0)  (9.19957e-250,0)  (1.88814e-278,0)  (5.82033e-312,0)             (0,0)
            (0,0) (-2.47284e-278,0)             (1,0)  (1.38744e-216,0)   (2.4762e-192,0)  (3.43923e-171,0)  (5.71164e-154,0)   (2.4711e-140,0) (-1.10218e-130,0) (-9.33605e-125,0)  (2.88109e-123,0)  (7.24901e-125,0)  (2.25024e-130,0) (-3.06856e-140,0)  (2.20528e-153,0)  (5.95593e-171,0) (-6.60291e-192,0)  (-3.7963e-218,0)  (6.44677e-246,0)  (1.88814e-278,0)             (0,0)
            (0,0) (-6.04933e-250,0) (-8.36347e-217,0)             (1,0)  (9.55555e-163,0)  (2.43362e-142,0)   (4.8234e-125,0)  (1.14309e-111,0)  (7.05723e-102,0)  (-4.49183e-96,0)  (-5.42948e-94,0)   (2.39098e-96,0)   (8.5847e-102,0)  (3.80277e-111,0)     (-7.4e-125,0)  (7.58904e-142,0)  (2.92481e-163,0) (-4.62711e-188,0)  (-3.7963e-218,0)  (9.19957e-250,0)             (0,0)
            (0,0) (-5.47394e-225,0) (-2.91959e-192,0) (-5.76007e-163,0)             (1,0)  (1.34013e-116,0)  (4.87047e-100,0)   (1.37752e-86,0)   (4.65852e-77,0)   (4.10421e-71,0)  (-3.72772e-69,0)   (-6.4299e-71,0)   (4.04062e-77,0)   (2.07025e-86,0)   (1.30865e-99,0) (-3.63396e-117,0)  (5.31816e-138,0)  (2.92481e-163,0) (-6.60291e-192,0) (-7.73057e-226,0)             (0,0)
            (0,0) (-1.83954e-204,0)    (-3.77e-171,0) (-2.86939e-142,0)  (-8.0783e-117,0)             (1,0)   (3.82729e-78,0)   (1.98491e-65,0)    (8.0111e-56,0)   (3.86606e-50,0)   (4.86044e-48,0)  (-6.29963e-50,0)  (-1.55061e-55,0)    (1.3905e-65,0)   (1.01665e-78,0)   (9.17059e-96,0) (-3.63396e-117,0)  (7.58904e-142,0)  (5.95593e-171,0) (-1.91872e-203,0)             (0,0)
            (0,0) (-1.66831e-186,0)  (-1.8079e-154,0)  (-5.2873e-125,0) (-5.74258e-100,0)  (-2.30708e-78,0)             (1,0)    (2.2258e-47,0)   (1.64726e-38,0)   (9.48721e-33,0)   (6.53342e-31,0)   (1.17212e-32,0)   (-2.1679e-38,0)  (-7.61466e-48,0)   (9.74419e-62,0)   (1.01665e-78,0)   (1.30865e-99,0)     (-7.4e-125,0)  (2.20528e-153,0)  (2.46975e-186,0)             (0,0)
            (0,0) (-8.48269e-173,0) (-2.33975e-140,0)  (-3.6182e-112,0)     (-1.51e-86,0)  (-2.34033e-65,0)  (-1.34171e-47,0)             (1,0)   (2.63592e-24,0)   (2.78377e-19,0)   (2.28789e-17,0)   (2.24835e-19,0)   (5.75601e-25,0)  (-1.51919e-34,0)  (-7.61466e-48,0)    (1.3905e-65,0)   (2.07025e-86,0)  (3.80277e-111,0) (-3.06856e-140,0)  (1.30494e-172,0)             (0,0)
            (0,0)  (3.65797e-163,0) (-1.69766e-130,0)  (-6.6821e-102,0)  (-1.47456e-77,0)  (-8.78157e-56,0)  (-1.94221e-38,0)  (-1.58893e-24,0)             (1,0)   (6.35667e-09,0)   (9.57979e-08,0)   (1.12353e-09,0)   (1.57557e-15,0)   (5.75601e-25,0)   (-2.1679e-38,0)  (-1.55061e-55,0)   (4.04062e-77,0)   (8.5847e-102,0)  (2.25024e-130,0) (-2.59112e-163,0)             (0,0)
            (0,0)  (1.44871e-157,0)  (1.04468e-124,0)  (-6.91864e-96,0)  (-3.88604e-71,0)  (-1.22372e-50,0)  (-1.03996e-32,0)  (-3.28223e-19,0)  (-3.83179e-09,0)       (1.00479,0)      (0.312161,0)   (0.000671321,0)   (1.12353e-09,0)   (2.24835e-19,0)   (1.17212e-32,0)  (-6.29963e-50,0)   (-6.4299e-71,0)   (2.39098e-96,0)  (7.24901e-125,0)  (2.71149e-157,0)             (0,0)
            (0,0)   (2.1493e-155,0)  (5.90405e-123,0)   (6.07546e-94,0)  (-5.74171e-69,0)  (-4.60208e-48,0)  (-2.06802e-31,0)  (-2.50793e-17,0)  (-1.12951e-07,0)      (-0.18817,0)       (34.5771,0)      (0.312161,0)   (9.57979e-08,0)   (2.28789e-17,0)   (6.53342e-31,0)   (4.86044e-48,0)  (-3.72772e-69,0)  (-5.42948e-94,0)  (2.88109e-123,0)  (1.24648e-155,0)             (0,0)
            (0,0) (-1.71726e-157,0)  (1.24995e-124,0)   (4.89971e-96,0)   (7.19491e-71,0)  (-9.70316e-50,0)  (-1.10982e-32,0)  (-7.11668e-20,0)  (-1.23159e-09,0)  (-0.000791528,0)      (-0.18817,0)       (1.00479,0)   (6.35667e-09,0)   (2.78377e-19,0)   (9.48721e-33,0)   (3.86606e-50,0)   (4.10421e-71,0)  (-4.49183e-96,0) (-9.33605e-125,0)  (7.06947e-158,0)             (0,0)
            (0,0)  (1.41319e-163,0) (-1.42514e-130,0)  (1.48026e-101,0)   (8.28022e-77,0)   (1.73509e-55,0)  (-3.33915e-38,0)  (-5.45005e-25,0)  (-4.98714e-16,0)  (-1.23159e-09,0)  (-1.12951e-07,0)  (-3.83179e-09,0)             (1,0)   (2.63592e-24,0)   (1.64726e-38,0)    (8.0111e-56,0)   (4.65852e-77,0)  (7.05723e-102,0) (-1.10218e-130,0) (-3.26903e-163,0)             (0,0)
            (0,0)   (6.1819e-173,0)  (1.67358e-140,0)  (-2.4084e-111,0)   (3.56972e-86,0)   (2.84948e-65,0)   (8.52063e-48,0)  (-2.33997e-34,0)  (-5.45005e-25,0)  (-7.11668e-20,0)  (-2.50793e-17,0)  (-3.28223e-19,0)  (-1.58893e-24,0)             (1,0)    (2.2258e-47,0)   (1.98491e-65,0)   (1.37752e-86,0)  (1.14309e-111,0)   (2.4711e-140,0) (-5.50726e-173,0)             (0,0)
            (0,0) (-4.01464e-186,0)  (1.04471e-153,0)  (4.03593e-125,0) (-8.28805e-100,0)   (1.75301e-78,0)   (1.99682e-61,0)   (8.52063e-48,0)  (-3.33915e-38,0)  (-1.10982e-32,0)  (-2.06802e-31,0)  (-1.03996e-32,0)  (-1.94221e-38,0)  (-1.34171e-47,0)             (1,0)   (3.82729e-78,0)  (4.87047e-100,0)   (4.8234e-125,0)  (5.71164e-154,0)  (1.76197e-186,0)             (0,0)
            (0,0) (-7.46538e-204,0) (-9.68154e-171,0)  (3.59515e-142,0)  (1.98195e-117,0)    (-5.808e-96,0)   (1.75301e-78,0)   (2.84948e-65,0)   (1.73509e-55,0)  (-9.70316e-50,0)  (-4.60208e-48,0)  (-1.22372e-50,0)  (-8.78157e-56,0)  (-2.34033e-65,0)  (-2.30708e-78,0)             (1,0)  (1.34013e-116,0)  (2.43362e-142,0)  (3.43923e-171,0)  (5.81158e-204,0)             (0,0)
            (0,0)  (1.54577e-224,0) (-2.56906e-192,0) (-4.75437e-163,0)  (2.51937e-138,0)  (1.98195e-117,0) (-8.28805e-100,0)   (3.56972e-86,0)   (8.28022e-77,0)   (7.19491e-71,0)  (-5.74171e-69,0)  (-3.88604e-71,0)  (-1.47456e-77,0)     (-1.51e-86,0) (-5.74258e-100,0)  (-8.0783e-117,0)             (1,0)  (9.55555e-163,0)   (2.4762e-192,0)  (4.99368e-225,0)             (0,0)
            (0,0)  (5.47223e-250,0)  (7.59093e-217,0) (-1.80032e-188,0) (-4.75437e-163,0)  (3.59515e-142,0)  (4.03593e-125,0)  (-2.4084e-111,0)  (1.48026e-101,0)   (4.89971e-96,0)   (6.07546e-94,0)  (-6.91864e-96,0)  (-6.6821e-102,0)  (-3.6182e-112,0)  (-5.2873e-125,0) (-2.86939e-142,0) (-5.76007e-163,0)             (1,0)  (1.38744e-216,0)  (5.13064e-250,0)             (0,0)
            (0,0)  (6.64721e-280,0)  (3.83476e-246,0)  (7.59093e-217,0) (-2.56906e-192,0) (-9.68154e-171,0)  (1.04471e-153,0)  (1.67358e-140,0) (-1.42514e-130,0)  (1.24995e-124,0)  (5.90405e-123,0)  (1.04468e-124,0) (-1.69766e-130,0) (-2.33975e-140,0)  (-1.8079e-154,0)    (-3.77e-171,0) (-2.91959e-192,0) (-8.36347e-217,0)             (1,0)  (4.10228e-278,0)             (0,0)
            (0,0) (-2.55093e-311,0)  (6.64721e-280,0)  (5.47223e-250,0)  (1.54577e-224,0) (-7.46538e-204,0) (-4.01464e-186,0)   (6.1819e-173,0)  (1.41319e-163,0) (-1.71726e-157,0)   (2.1493e-155,0)  (1.44871e-157,0)  (3.65797e-163,0) (-8.48269e-173,0) (-1.66831e-186,0) (-1.83954e-204,0) (-5.47394e-225,0) (-6.04933e-250,0) (-2.47284e-278,0)             (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)             (1,0)
4

Thanks @WardBrian @stevebronder, any other ideas on this?

FYI – taking only the real component of the matrix and computing an eigendecomposition does not produce NaNs.

Also, @stevebronder, the complex component of that matrix is machine zero, but that shouldn’t in general be true for kmat_m, the matrix causing problems.

5

I finally got a reproducing example in pure Eigen, so I’ve opened ComplexEigenSolver produces NaNs when original matrix has none (#2773) · Issues · libeigen / eigen · GitLab

1 Like
6

@pgree - just to keep you up to date, the Eigen maintainers have made two changes which combined seem to have resolved this issue. One is already scheduled for inclusion in 3.4.1, and I’ve requested that they consider including the other as well, so hopefully this will be resolved in the next Eigen update.

In the mean time, another maintainer suggested balancing as a solution (and most likely as the reason numpy etc were getting the right answer before): c++ - Eigen - Balancing matrix for eigenvalue - Stack Overflow

1 Like
7

thanks for your help and the update @WardBrian