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