Error when refining SPDE mesh in marked LGCP using inlabru with HPC

[Morgane Declerck]

 Dear INLA team,  

I am fitting a marked LGCP using inlabru, where the point process models observation locations and the mark model is a binomial response. 

The model includes cyclic random slopes across four tidal phases.

• The mark model includes both:

  • point_field * beta_link (shared spatial field linking the point process and mark model),

  • mark_field_tide (tide-specific spatial residual field).

• The model is run on an HPC cluster using Linux.

The model fits successfully with a coarser spatial mesh, but fails with a finer mesh (7–10 km resolution) with a singular matrix error.

The model runs reliably with the 15–20 km mesh but consistently fails with the 7–10 km mesh.

I am trying to understand whether this issue is due to numerical instability introduced by the finer mesh or a problem with the model specification.

The job runs on an HPC cluster (Linux) and fails with the following message:

processing file: SIM_lgcp_v4_7_2.Rmd

GitId: "devel"
Error:3 Reason: Matrix is (numerical) singular
Line:343 Function: GMRFLib_comp_posdef_inverse

Segmentation fault (core dumped)
'/uoa/home/.../INLA/bin/linux/64bit/inla.mkl.run'

The model runs correctly with the following mesh:

mesh <- fm_mesh_2d_inla(
  loc      = st_coordinates(gan_m),
  boundary = domain_m,
  max.edge = c(15000, 20000),
  cutoff   = 15000,
  offset   = c(10000, 40000),
  crs      = st_crs(gan_m)
)

However, when refining the mesh to approximately 7–10 km resolution, the model crashes:

mesh <- fm_mesh_2d_inla(
  loc      = st_coordinates(gan_m),
  boundary = domain_m,
  max.edge = c(7000, 10000),
  cutoff   = 5000,
  offset   = c(10000, 40000),
  crs      = st_crs(gan_m)
)

spde <- inla.spde2.pcmatern(
mesh = mesh,
prior.range = c(150000, 0.5),
prior.sigma = c(1, 0.05)
)

The model is a marked LGCP fitted with bru().

Point process

geometry ~ Inter_point + point_field

Mark model

mark ~ Inter_mark +
log_PEA_z * slope_PEA +
log_surface_chla_z * slope_CHLA +
sst_z * slope_SST +
log_MLD_z * slope_MLD +
sal_depth_int_z * slope_SAL +
log_dist_col_m_z +
point_field * beta_link +
mark_field_tide

Latent components

cmp <- ~ -1 +
point_field(geometry, model = spde) +

slope_PEA(
tide_group_id,
model = "generic0",
Cmatrix = Q_cyclic,
constr = TRUE,
rankdef = 1
) +

slope_CHLA(
tide_group_id,
model = "generic0",
Cmatrix = Q_cyclic,
constr = TRUE,
rankdef = 1
) +

slope_SST(
tide_group_id,
model = "generic0",
Cmatrix = Q_cyclic,
constr = TRUE,
rankdef = 1
) +

slope_MLD(
tide_group_id,
model = "generic0",
Cmatrix = Q_cyclic,
constr = TRUE,
rankdef = 1
) +

slope_SAL(
tide_group_id,
model = "generic0",
Cmatrix = Q_cyclic,
constr = TRUE,
rankdef = 1
) +

Inter_point(1) +
Inter_mark(1) +

beta_link(
main = 1,
model = "linear",
mean.linear = 0,
prec.linear = 1,
n = 1,
values = 1
) +

mark_field_tide(
geometry,
model = spde,
replicate = tide_group_id
)

Is this singular matrix error most likely due to:

1. the much finer mesh resolution,

2. the combination of point_field * beta_link and mark_field_tide.

3. the number of cyclic slope terms relative to only four tide groups,

4. or another numerical issue related to the SPDE discretisation?

Are there recommended strategies to stabilise the model when refining the mesh (for example changes to mesh construction, priors, or simplifying the latent model structure)?

Many thanks for any guidance.

Morgane

[Finn Lindgren]

Hi Morgane,

the first thing to try would be to change the modelling units from metres to kilometres. This should make the internal numerics more stable.

Also, don't use st_coordinates() to extract coordinates for input to fm_mesh_2d, as that loses the crs information. Just supply the actual sf object, so fm_mesh_2d can handle it appropriately.

I'd also make the outer mesh boundary explicit instead of the basic construction from "offset", and reduce the "cutoff" value to make it less likely to interfere with the max.edge setting.

crs_km <- fm_crs(gan_m, units = "km")
gan_km <- fm_transform(gan_m, crs_km)

domain_km <- fm_transform(domain_km, crs_km)

outer_boundary <- fm_nonconvex_hull(domain_km, 40)

mesh <- fm_mesh_2d_inla(
  loc      = gan_km,
  boundary = list(domain_km, outer_boundary),
  max.edge = c(15, 20),
  cutoff   = 10,
  offset   = c(10, 40), # Can likely be omitted.
  crs      = crs_km
)

mesh <- fm_mesh_2d_inla(
  loc      = gan_km,
  boundary = list(domain_km, outer_boundary),
  max.edge = c(7, 10),
  cutoff   = 5,
  offset   = c(10, 40),
  crs      = crs_km
)

(Or alternatively, fine_mesh <- fm_subdivide(mesh, 1), that will yield a similar mesh with edges aligned with the original mesh)

I would also recommend considering not specifying "loc" _at all_, or to specify a regular triangle grid, like

mesh <- fm_mesh_2d_inla(

  loc = fm_hexagon_lattice(domain_km, edge_len = 15),

  boundary = list(domain_km, outer_boundary),

  max.edge = c(16, 20), # The max.edge value(s) need to be larger than the hexagon lattice edge_len.
  cutoff   = 10,
  offset   = c(10, 40), # Can likely be omitted.
  crs      = crs_km
)

This both improves the spde approximation and the point process likelihood construction, and improves the numerics.

If there are still numerical singularity issues after that change, one can revisit the other potential issues.

I believe you can simplify your component definitions, such as using the weights argument (always the second unnamed argument):

log_PEA(
  tide_group_id,

  log_PEA_z,

  model = "generic0",
  Cmatrix = Q_cyclic,
  constr = TRUE,
  rankdef = 1
)

and then replacing log_PEA_z * slope_PEA with log_PEA in the predictor formula expression.

If this is done for all the components, the pre-processing in inlabru will be a bit faster.

Finn

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/c341b7f1-4046-4e10-b141-b7bb23ac2711n%40googlegroups.com.

--

Finn Lindgren
email: finn.l...@gmail.com

[Morgane Declerck]

Hi Finn,

Thank you again for your suggestions. I implemented the changes you recommended and the model now runs to completion, which is a clear improvement compared to before.

In particular, I:

1. changed the modelling units from metres to kilometres,
2. avoided using st_coordinates() when constructing the mesh,
3. made the outer mesh boundary explicit using fm_nonconvex_hull(),
4. reduced the cutoff,
5. and simplified the tide-varying effects using the weighted component formulation.

The relevant parts of the workflow now look like this:

crs_km    <- fm_crs(gan_m, units = "km")
gan_km    <- fm_transform(gan_m, crs_km)

domain_km <- fm_transform(domain_m, crs_km)

outer_boundary <- fm_nonconvex_hull(domain_km2, 40)

mesh <- fm_mesh_2d_inla(
  loc      = gan_km,

  boundary = list(domain_km2, outer_boundary),

  max.edge = c(7, 10),
  cutoff   = 5,

  crs      = crs_km

)

spde <- inla.spde2.pcmatern(
  mesh = mesh,

  prior.range = c(150, 0.5),

  prior.sigma = c(1, 0.05)
)

  The component structure now follows the weighted formulation:

  cmp <- ~ -1 +

  point_field(geometry, model = spde) +

  log_PEA(tide_group_id, log_PEA_z,
          model="generic0", Cmatrix=Q_cyclic, constr=TRUE, rankdef=1) +

  log_surface_chla(tide_group_id, log_surface_chla_z,
          model="generic0", Cmatrix=Q_cyclic, constr=TRUE, rankdef=1) +

  sst_effect(tide_group_id, sst_z,
          model="generic0", Cmatrix=Q_cyclic, constr=TRUE, rankdef=1) +

  log_MLD(tide_group_id, log_MLD_z,
          model="generic0", Cmatrix=Q_cyclic, constr=TRUE, rankdef=1) +

  sal_effect(tide_group_id, sal_depth_int_z,
          model="generic0", Cmatrix=Q_cyclic, constr=TRUE, rankdef=1) +

  Inter_point(1) +
  Inter_mark(1) +

  beta_link(main = 1, model = "linear",

            mean.linear = 0, prec.linear = 1,
            n = 1, values = 1) +

  log_dist_col_m_z +

  mark_field_tide(geometry, model = spde, replicate = tide_group_id)

  The model is then fitted using:

   INLA::inla.setOption(num.threads = "1:1")

fit <- bru(
  components = cmp,
  obs_points,
  obs_marks,
  options = list(
    control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE),
    control.predictor = list(compute = TRUE),
    control.inla = list(
      strategy = "simplified.laplace",
      int.strategy = "auto"
    ),
    bru_verbose = TRUE
  )
)

I also tested the alternative mesh construction using a regular lattice (fm_hexagon_lattice(...)) instead of loc = gan_km.

In both mesh versions the model runs, but during fitting I repeatedly obtain the warning:

*** warning *** iterative process seems to diverge, 'vb.correction' is aborted
*** Please (re-)consider your model, priors, confounding, etc.

Given that the warning appears with both mesh constructions, I was wondering how you would interpret it. In particular:

1. How would you recommend addressing this warning?

2. Does it mainly indicate a numerical issue with the variational Bayes correction step, or could it suggest a more fundamental issue with the model specification?

3. If the model otherwise runs to completion, to what extent might this warning affect the reliability of the outputs (fixed effects, spatial fields, predictions, WAIC/DIC, etc.)?

4. Since the warning occurs with both loc = gan_km and the hexagonal lattice mesh, would your next step be to simplify the model structure (for example removing the tide-phase replication on mark_field_tide), or would you first reconsider priors or integration settings?

Many thanks again for your help — your previous suggestions already allowed the model to run successfully.

Best wishes,

Morgane

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/CAHC4gzY_ubs9HoFYQrZGE2Z7i27DiMqDznzCPPFiQV2GqANVPg%40mail.gmail.com.