spde for R3 mesh

[Ran]

Dear INLA-helpers,

I calculated the \tilde{C}, G and B matrices for a tetrahedral mesh in R3 and would like to try INLA for my model.
When I create a spde object with boundary segments, I suppose that spde$internal$c0 is the \tilde{C} matrix, spde$internal$g1 is the B matrix and spde$internal$g2 is the G matrix due to its numbers of nonzeros in the sparse matrices, am I right?

I manipulated some of the parameters of the spde object like the following example:
spde$internal$c0 = Ct
spde$internal$g1 = B
spde$internal$g2 = G
spde$f$n = dim(Ct)[1]

Also I changed some of the graph attributes like
spde$mesh$graph$tv = facets (nx4, n = number of tetrahedra)
spde$mesh$graph$vt = vertices
...

Since I just blindly manipulate some inputs, I suppose, INLA would not work (certainly your codes are not simple). I am just wondering if there is a quick solution for using INLA with a R3 mesh (like changing some input parameters for S2 spde objects to make a R3 spde object) or you have to implement the whole thing..

Thank you very much for your help! I really appreciate it!
Ran.

[Finn Lindgren]

Hi,

the documentation for this advanced stuff is very much not written yet...
Some notes below:

On Tuesday, 11 September 2012 11:53:27 UTC+1, Ran wrote:

I calculated the \tilde{C}, G and B matrices for a tetrahedral mesh in R3 and would like to try INLA for my model.
When I create a spde object with boundary segments, I suppose that spde$internal$c0 is the \tilde{C} matrix, spde$internal$g1 is the B matrix and spde$internal$g2 is the G matrix due to its numbers of nonzeros in the sparse matrices, am I right?

No, not quite.  g1 is G and g2 is G*inv(C)*G, and we don't use B in the inla models at all, since we've only implemented Neumann boundary conditions.  Replacing G with G-B as in the theory unfortunately leads to strange models, so avoid that for now. I manipulated some of the parameters of the spde object like the following example:

 

spde$internal$c0 = Ct
spde$internal$g1 = B
spde$internal$g2 = G
spde$f$n = dim(Ct)[1]

Unfortunately, that won't work, due to the way these objects interact with the inla program. 

 

Since I just blindly manipulate some inputs, I suppose, INLA would not work (certainly your codes are not simple). I am just wondering if there is a quick solution for using INLA with a R3 mesh (like changing some input parameters for S2 spde objects to make a R3 spde object) or you have to implement the whole thing..

Yes, fortunately there is an alternative method for constructing more generic spde models that you can use instead, and which should be much easier to deal with, and which doesn't need to know anything about your mesh:
inla.spde2.generic

This functions constructs inla model objects with a specific precision matrix structure, where one can specify the different building blocks.  This is used by the inla.spde2.matern function to create Matern models on ordinary triangulation meshes, and what you need is a corresponding wrapper function for your 3D meshes.

When you have your computed C and G matrices, the following code should produce an SPDE model with alpha=2, i.e.
(kappa^2 - Laplacian ) (tau x(s)) = W(s)

spde = inla.spde2.generic(
 M0 = C,
 M1 = G,
 M2 = G %*% inla.qsolve(C, G),
 B0 = cbind(0, 1, 0), ## Make theta1 control log(tau)
 B1 = cbind(0, 0, 2), ## Make theta2 control log(kappa)
 B2 = cbind(1, 0, 0), transform="identity"), ## Unused part of the generic model; don't change
 theta.mu = c(0, 0), theta.Q = diag(2) )

Then you'll likely need some more details to use this in an inla model, but hopefully this can get you started.
For example, you can check if the models seems sensible by simulating a sample:

sample = inla.spde.sample(spde, theta=log(c(some.tau, some.kappa)))

The interpretation of kappa should in your 3D case roughly be sqrt(8*0.5) / spatial.range

/Finn

[Finn Lindgren]

Note: in my example code, it's important that C is diagonal, i.e. the matrix you called tilde{C}.

/Finn

[Ran]

Hi Finn,

thank you very much for the detailed and constructive feedback!

My comments:
1) INLA worked for some sampled fields using the funciton inla.spde2.generic:
I ran few simulations for spatial.range = 1, 100 and 1000, and with the real data. INLA worked for spatial.range = 1 and 100, however, for spatial.range = 1000, as well as when using real data, INLA failed to converge.

Is the fail of convergence related to the issue with the intrinsic field, since kappa in those cases may be close to zero?.. In my case, kappa needs to be small since the entries of my C and G matrices are very large. They are calculated based on volumes of the tetrahedra/triangles of the mesh in which the minimum distance between two nodes are 70km. 
Which unit do you usually use to calculated the matrices from a mesh , meter^2, kilometer^2? (For example, in Fig 5 in your paper Lindgren et al (2011) you plot the covariance function against the great circle distance. is the unit in km?)

Although INLA worked for spatial.range = 1 and 100, the estimation of the fields is not accurate..the values are too small. I think that this issue is related to strong regularization since the values of the C and G matrices are very large.

2) A small question: while using the function inla.qsolve(C,G), I noticed that it ran out of memory on my own laptop (RAM 4GB). Does the function convert the input matrices into full matrices? (C is diagonal and G is sparse)

Thank you very much for the help!
Ran

[Ran]

I think that I need to normalized the distance between two nodes by min.dist/6378137 while calculating my C and G matrices, is it correct?

[Finn Lindgren]

On Saturday, 15 September 2012 14:30:55 UTC+1, Ran wrote:

Is the fail of convergence related to the issue with the intrinsic field, since kappa in those cases may be close to zero?.. In my case, kappa needs to be small since the entries of my C and G matrices are very large. They are calculated based on volumes of the tetrahedra/triangles of the mesh in which the minimum distance between two nodes are 70km. 

Yes, the distance units used for the model can have a large impact on the numerical properties, and this seems to exactly the case in your setting.
 

Which unit do you usually use to calculated the matrices from a mesh , meter^2, kilometer^2? (For example, in Fig 5 in your paper Lindgren et al (2011) you plot the covariance function against the great circle distance. is the unit in km?)

I almost always use a unit radius sphere for my spherical surface meshes, and then rescale the output to desired plotting scale; The x-axis in Fig 5 has units in radians [0, pi], whereas the y-axis in Fig 8(c) is in km, with values obtained by rescaling radians into km by multiplying with the earth radius.
 

Although INLA worked for spatial.range = 1 and 100, the estimation of the fields is not accurate..the values are too small. I think that this issue is related to strong regularization since the values of the C and G matrices are very large.

Yes. 

 

2) A small question: while using the function inla.qsolve(C,G), I noticed that it ran out of memory on my own laptop (RAM 4GB). Does the function convert the input matrices into full matrices? (C is diagonal and G is sparse)

Ah, sorry about that; hadn't tested that code on large matrices, and forgot that this isn't necessarily indended to work (the second parameter is intended to typically be a vector, not a matrix...). INstead,
  solve(C, G)
should likely do the right thing, as will
  Diagonal(nrow(C), 1/diag(C)) %*% G
since C must be diagonal.

/Finn

[Finn Lindgren]

On Saturday, 15 September 2012 18:42:31 UTC+1, Ran wrote:

I think that I need to normalized the distance between two nodes by min.dist/6378137 while calculating my C and G matrices, is it correct?

Any rescaling factor that allows the numerical calculations to work are ok, so if that factor helps you should be ok.  As I wrote in my other reply, I usually rescale my distances to a unit sphere, i.e. rescale by the radius of the earth, so 1/6378.137 for the equatorial radius in km.  I don't know if also multiplying with the min.dist helps; do let me know if that is the case!

/Finn

[Ran]

Hi Finn,

I hope you are good, Happy new years!

Quick questions on creating non-stationary SPDE model using the inla.spde2.generic function: From the examples I found in diverse tutorials and forum posts you first create basis functions via inla.mesh.basis which requires a mesh and then use the function inla.spde2.matern to create the SPDE object. 

My questions:

0) Is it possible to generate the basis functions for tau and kappa using some INLA functions (maybe with inla.fmesher.smorg?) when I only have the coordinates of the vertices (=mesh$loc?)  and facets (= mesh$graph$tv?) of the mesh (I do not have the mesh in the inla-class)? 

1) If one got the basis functions for tau and kappa, is it possible to put those values in inla.spde2.generic function, specified in the arguments B0, B1 or B2?

Thank you very much for your help!

Best,Ran

On Wednesday, September 12, 2012 4:10:58 PM UTC+2, Finn Lindgren wrote: