Kronecker product in JAGS

  • rrankin

    rrankin - 2012-08-08

    Hello all,

    JAGS has some good matrix operations which other BUGS-esque programs lack, but
    so far all programs (that I know of) cannot compute the Kronecker product.
    Shouldn't it be theoretically possible, given JAGS's use of C++ matrix
    libraries? (if someone has the time and sophistication to build a module).

    My questions:

    1) can the Kronecker product be computed "manually" through a series of
    matrix-scalar multiplications? e.g.,

    Omega <- Sigma1 * Sigma2

    ... or something along those lines. I guess it boils down to whether the
    inheriting matrix Omega can be built-up from a series of sub-matrices.

    2) is there a resource out there to refer to how to compile such a module from
    C++ code?



  • rrankin

    rrankin - 2012-08-08

    ... also, I should say that the application of the above is for separable
    spatial-temporal processes, where either process has its own covariance
    matrix, and are subsequently combined via the Kronecker product. Thereafter,
    time-and-space indexed observations fall in a single 1-column vector and are
    drawn from a multivariate distribution with the Kronecker product as the
    covariance matrix. The following has a good review of its properties and

    Clark, J.S. and Gelfand, A.E. (eds.). 2006. Hierarchical Modelling for the
    Environmental Sciences: Statistical Methods and Applications. Oxford
    University Press, Oxford. i-x+205pp.

  • Martyn Plummer

    Martyn Plummer - 2012-08-31

    I'm sure you can define a Kronecker product element-wise.

  • rrankin

    rrankin - 2012-08-31


    Yes, I am able to define the kronecker product with element-wise
    multiplication. However, this triggers another error: the script fails to draw
    a sample for the separable covariance matrices which compose the kronecker
    product, with error "unable to find sampler".

    Delta ~ dwish(R, S)         # covariance matrix 1 (spatial)
    p ~ dbeta(1,1)                    # AR(1) parameter, to build correlation matrix 2
    Sigma <- p**distance_matrix    # correlation matrix 2, (temporal)
    # element-wise kronecker product
    for (i in 1:d){          # d = width of Delta
        for (j in 1:d){
            for (k in 1:s){     # s = width of Sigma 
                for (l in 1:s){
                     Kronecker[(i-1)*s+k,(j-1)*s+l] <- Delta[i,j]*Sigma[k,l]
    # use kronecker as covariance matrix 
    Omega <- inverse(Kronecker)
    Y[1:n] ~ dmnorm(zeros[1:n], Omega)

    E.g., it seems we cannot sample from a Wishart prior and then use the sample
    in the Kronecker production.

    Thank you,


  • tds151

    tds151 - 2012-08-31

    the issue is not the computation of the kronecker product. jags doesn't
    recognize the conjugacy of such a model (which is conjugate) and, thus, is not
    able to sample the posterior for a positive definite matrix from a non-
    conjugate formulation. this issue generalizes broadly to any conjugate
    specification for a covariance matrix, such as gelman's parameter expanded
    wishart that intends to reduce the correlations among the variance and
    covariance parameters. Dr. Plummer indicates it is possible to add such
    flexibility to JAGS, but it isn't a priority. You may consider Stan as an
    alternative, though it will be somewhat harder to work with,

  • Martyn Plummer

    Martyn Plummer - 2012-09-03

    I see the problem now. It is not just a question of defining a new function but getting the conjugate Wishart sampler to recognize it and do the appropriate calculations. This is a little more work but is probably worth doing.

    If you want me to work on this you need to send me a simple example.

  • tds151

    tds151 - 2012-10-03

    I have an example, now, including data. In any case, I appreciate your time to consider this and am more than happy to produce an example. I'd ask for a more guidance on the components that compose an example of value to you. For example, do you want me to .tex up a formulation for a specific inferential problem? Would you rather I conveyed something in pseudo-JAGS script? The example I have in my mind is my current project assessing the performance for teacher lesson delivery over mutliple classrooms of students, evaluation events and raters. The evaluation protocol includes multiple dimensions that compose lesson delivery performance, so it produces multivariate ordered responses. The evaluation events may be video-taped and later scored or scored "live". Our inferential task is to evaluate a whether these two different scoring modes express similar teacher-level correlations. So inside the GLM we define a set of multivariate teacher random effects as the mean of a latent continuous response (or the intrinsic traits of an ordered logit formulation). We have D scoring dimensions per mode and M = 2 modes. I define two covariance matrices, a D x D dimension covariance and an M x M mode covariance. The teacher effects are defined over now MD = 2D dimensions for each of T teachers with covariance composed as the kronecker product of the dimension and mode covariances. These lower dimensional precision matrices would receive wishart priors and the teacher random effects a multivariate Gaussian. An alternative formulation that avoids the need to hand-compose a kronecker product is to employ the matrix variate normal. (Dawid-1981). Lastly, an alternative to this formulation that would be very valuable for me would be to replace the inverse wishart prior on the dimension covariance matrix with a sparse, factor analytic contruction. In particular, S = LL' + Theta, where L would be restricted to lower triangular (which I do now in JAGS through a trick) and the elements would receive N(0,1) priors. Then the woodbury matrix identity could be used to reduce the dimension for computing the inverse of S, the dimension covariance.

    Last edit: tds151 2012-10-03
  • Maura Mezzetti

    Maura Mezzetti - 2013-07-04

    I am encountering the same problem as rrankin, were you able to solve it?

  • Martyn Plummer

    Martyn Plummer - 2013-07-08

    No, sorry no progress here.


Log in to post a comment.