R crashes using the matern function

BertvanderVeen commented 3 years ago

Description:

While attempting to fit a spatial model in TMB, with matern covariance, R crashes when maximizing. When switching from matern to exp(-phi*H(i,j)) things work fine.

Reproducible Steps:

#include <TMB.hpp>
#include<math.h>

template<class Type>
Type objective_function<Type>::operator() ()
{
  using namespace Eigen;
  using namespace density;

  DATA_ARRAY(y); 
  DATA_MATRIX(H);
  DATA_ARRAY(NAidx);

  DATA_VECTOR(time);

  PARAMETER(lg_phi);
  PARAMETER(lg_tau);
  PARAMETER_VECTOR(lg_sigma);
  PARAMETER_VECTOR(beta);
  Type phi = exp(lg_phi);
  Type tau = exp(lg_tau);

  PARAMETER_VECTOR(b);
  PARAMETER_VECTOR(lambda);

  PARAMETER_ARRAY(u);

  parallel_accumulator<Type> nll(this);
  int Tt = u.cols();
  int num_lv = u.col(0).cols();
  int n = u.col(0).col(0).size();
  int p = y.col(0).col(0).size();

  matrix<Type> newlam(num_lv,p);

  for (int j=0; j<p; j++){
    for (int i=0; i<num_lv; i++){
      if (j < i){
        newlam(i,j) = 0;
      } else{
        newlam(i,j) = lambda(j);
        if (i > 0){
          newlam(i,j) = lambda(i+j+i*p-(i*(i-1))/2-2*i);
        }
      }
    }
  }

  matrix <Type> C(n,n);
  for (int i=0;i<n;i++)  {
    C(i,i)=Type(1);
    for (int j=0;j<i;j++)    {
      C(i,j) = matern(H(i,j), phi, tau);
      C(j,i) = C(i,j);
    }
  }

  MVNORM_t <Type> neg_den(C);
  for (int t=0; t<Tt; t++){
    for (int q=0; q<num_lv; q++){
        nll += neg_den(u.col(t).col(q));
    } 
    }

  array<Type> eta(n,p,Tt);
  eta.fill(0.0);
  for (int t=0; t<Tt; t++){
    for (int j=0; j<p; j++){
      for (int i=0; i<n; i++){
        eta(i,j,t) += b(j) + time(t)*beta(j);
      }
    }
  }

  for (int t=0; t<Tt; t++){
    for (int j=0; j<p; j++){
      for (int i=0; i<n; i++){
    for (int q=0; q<num_lv; q++){
      eta(i,j,t) += u(i,q,t)*newlam(q,j);
    }
      }
    }
  }

  vector <Type> sigma = exp(lg_sigma)/(Type(1.0)+exp(lg_sigma));
  for (int t=0; t<Tt; t++){
        for (int j=0; j<p; j++){
          for (int i=0; i<n; i++){
          if(NAidx(j,i,t)<1)nll -= dzipois(y(j,i,t),exp(eta(i,j,t)), sigma(j),true);
    }
        }
  }

  return nll;
}

Current Output:

Expected Output:

-

TMB Version:

1.7.16

R Version:

4.0.2,3.6.2, with microsoft R open with intel MKL

Operating System:

Windows 10, Ubuntu 18.04

kaskr commented 3 years ago

Likely duplicate of #189 (does it solve the problem to set openmp(1) ?).

matern calls the special functions lgamma and BesselK which both might throw warnings for extreme (large or small) parameters. Until #189 has been solved a possible workaround is to add a validpar implementation to your obj - see ?MakeADFun.

BertvanderVeen commented 3 years ago

That makes sense, exp(lg_phi) tends to infinity for my example. Setting openmp(1) does not solve the problem (in addition, it happens on vanilla R too). Thanks for the workaround, but do you happen to have a recommendation of a valid parameter range too?

kaskr commented 3 years ago

but do you happen to have a recommendation of a valid parameter range too?

No, it's problem specific

kaskr / adcomp