Geometric anisotropy

[James Thorson]

Hi all,

I'm writing to inquire whether there is documentation for how to implement "geometric anisotropy" in R-INLA and/or the SPDE approach?  Geometric anisotropy is obviously a very frequent model to explore when isotropy fails to explain the data well and, while I see some discussion of general anisotropy and nonstationarity, it would be convenient to have this intermediate level to explore.

Cheers,

Jim

[Daniel Simpson]

Hi James,

Right now we do not have that implemented. 

We've written a few papers on it, but the challenge is coming up with a useful parameterisation that has few enough hyperparameters to be useful in INLA.

Conceivably, we could have a model of stationary anisotropy, but we haven't implemented that either.

These papers are listed here: http://arxiv.org/find/stat/1/au:+Fuglstad_G/0/1/0/all/0/1

Best wishes,

Dan

--
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 post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at http://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/groups/opt_out.

---

Daniel Simpson

Researcher,

Department of Mathematical Sciences,

NTNU, Trondheim, Norway.

Room: SB-II 1222

Phone: +47 941 74 803

Email: dp.si...@gmail.com

Email: daniel....@math.ntnu.no

Skype: dp.simpson

Website: http://math.ntnu.no/~daniesi/

[James Thorson]

Dan,

At the risk of trying your patience:

do you have an example of how to calculate the precision matrix arising from Eq. 6 of Fuglstad et al. (arXiv) that you linked?  

For isotropy, I've previously computed:

Q = kappa^4*G0 + 2*kappa^2*G1 + G2

where:

G0=spde$param.inla$M0

G1=spde$param.inla$M1

G2=spde$param.inla$M2

and feeding the precision matrix Q into other software (i.e. Template Model Builder).  Is there an analogous way to compute Q given the parameterization of geometric anistropy H in that paper and the mesh computed by INLA?

Cheers,

Jim

[Daniel Simpson]

Not really. You have to program it up seperatly (G1 is very different in this case).  The appendices of the first paper (chronologically) give the details of how to code it up yourself.  Alternatively, you could use some finite elements software, although I don't have any specific suggestions about that.

D

[James Thorson]

Dan,

Thanks for your previous assistance.

I'm still trying to code up anisotropy in TMB using the SPDE approximation (as obviously inspired by INLA).  I couldn't find your suggested approach in Fuglstad arXiv articles, but have made some progress via the Appendix A in Lindgren et al. 2011 (Proc R Soc Stat).  Specifically, I'm trying to re-code the computation of M0, M1, and M2 matrices given a particular triangularization, so that I can then include the H-matrix with estimated parameters (Eq. 20 of Lindgren et al. 2011).  

However, I can't work out the notation for the boundary conditions, i.e., how the B-matrix is used when computing M0, M1, and M2.  The re-coding in R (shown for easy debugging) is attached.

Is anyone able to modify this as necessary to include the B-matrix?  Or provide a reference that explains the use of the B-matrix more explicitly when computing M0/M1/M2?

I hope this R-code would be a useful reference for others trying to understand the details of the SPDE approximation, and therefore hope it is worth the effort.  

Jim

D

Dan

To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion-group+unsub...@googlegroups.com.

To post to this group, send email to r-inla-disc...@googlegroups.com.

Visit this group at http://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/groups/opt_out.

[James Thorson]

Hi all,

In case anyone is following this thread, I thought I'd post a final version of the R-code that replicates the computation of M0, M1, and M2 from R-INLA (i.e., the "inla.spde2.matern" function) given a particular triangular mesh.  This code *appears* to allow inclusion of the H-matrix to estimate geometric anisotropy.  However, I still welcome feedback on:

1.  Whether there's some mistake in my interpretation of the H-matrix

2.  How the B-matrix from Lindgren et al. (2011) enters into the computation (I calculate the B-matrix but it doesn't appear necessary to replicate the output from "inla.spde2.matern").

Cheers,

Jim

[Finn Lindgren]

On 21/05/14 23:34, James Thorson wrote:
> Hi all,
>
> In case anyone is following this thread, I thought I'd post a final
> version of the R-code that replicates the computation of M0, M1, and M2
> from R-INLA (i.e., the "inla.spde2.matern" function) given a particular
> triangular mesh. This code *appears* to allow inclusion of the H-matrix
> to estimate geometric anisotropy. However, I still welcome feedback on:
>
> 1. Whether there's some mistake in my interpretation of the H-matrix

For stationary anisotropy, take a look at equation (20) of Appendix A of
the Lindgren et al (2011) JRRSB SPDE/GMRF paper. When the SPDE is
defined as

(\kappa^2 - \nabla \cdot H \nabla) x(s) = W(s)

the matrix in the G-element calculation should be the "adjugate" of H,
in the non-singular case (the only one relevant here) defined as

adj(H) = det(H) H^{-1}

With that modification your code looks correct.

> 2. How the B-matrix from Lindgren et al. (2011) enters into the
> computation (I calculate the B-matrix but it doesn't appear necessary to
> replicate the output from "inla.spde2.matern").

The B-matrix is only used for for special constructions when no using
Neumann conditions, so it's correct that is doesn't and shouldn't affect
the inla.spde2.matern() output in the normal cases.

One important note on your code however is that you're not using sparse
matrices; while it's possible to convert non-sparse storage into sparse,
you still use a lot more memory while building the matrices than is
necessary. It's much more efficient to first build index-value triplets
(i,j,value), and build the sparse matrix using sparseMatrix()
as a post-processing step. This will also automatically accumulate the
values for duplicated elements, so no special code for that is needed.
I've attached Matlab code that uses this approach; the syntax is very
similar to what one would need in R; the Matlab sparse() function does
exactly the same thing as R's sparseMatrix() from the Matrix package.

Also note that the internal storage in spde$param.inla might change, so
for a safer comparison of the matrices, compare with the output from

inla.fem(mesh, order=2)

instead, as that is the more direct interface to the C++-code used
internally by R-INLA.


Experimental non-stationary anisotropic code is in the fmesher program:

https://code.google.com/p/inla/source/browse/fmesher/src/mesh.cc#2370

and can be accessed through the internal "smorgasbord" function:

gamma = c(...) ## length n vector of baseline diffusion
vec = matrix(...) ## n-by-3 matrix of anisotropy vectors
smorg = (INLA:::inla.fmesher.smorg(mesh.globe$loc, mesh.globe$graph$tv,
fem=2,
aniso=list(gamma,vec),
output=c("g1aniso","g2aniso")))
g1 = smorg$g1aniso
g2 = smorg$g2aniso

In the stationary case, gamma should be rep(something, mesh$n)
and vec should have mesh$n equal rows, with third column equal to zero
(for planar meshes; on the sphere they should be tangential).

Finn L