Construct a row-wise kronecker product matrix

Modeling
1

Dear Stan users:

I have been having difficulty in constructing a row-wise Kronecker product in Stan.

Suppose my two matrices are

A<- matrix(c(1,2,3,4,0,0),nrow=2,byrow=TRUE)
B <- matrix(c(0,5,6,7),nrow=2,byrow=TRUE)

Then a row-wise Kronecker product in R can be coded as

A[, rep(seq(ncol(A)), each = ncol(B))] * B[, rep(seq(ncol(B)),ncol(A))] 

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    5    0   10    0   15
[2,]   24   28    0    0    0    0

Now, I would like to construct the same thing using Rstan, but I cannot get my head around the index. I know how to construct the normal Kronecker in Stan between two matrices

The following code credited a fellow Stan user Tanner Sorensen

data{
  int<lower=1> m; // in GAM model, m=p=num_data
  int<lower=1> n;
  int<lower=1> p; //in GAM model, m=p=num_data
  int<lower=1> q;
  matrix[m,n] A;
  matrix[p,q] B;
}
parameters{real<lower=0,upper=1> DUMMY;}
model{}
generated quantities{
  matrix[m*p,n*q] C;
  for (i in 1:m)
    for (j in 1:n)
      for (k in 1:p)
        for (l in 1:q)
          C[p*(i-1)+k,q*(j-1)+l] = A[i,j]*B[k,l];
}

As the row-wise Kron matrix in R is

B_true_int= matrix(nrow=dim(A)[1], ncol=dim(A)[2]*dim(B)[2])
for(i in 1:dim(A)[1]){
  B_true_int[i,]=kronecker(A[i,],B[i,])
}

Then my initial thought was to replicate Tanner’s code through first extracting the row of A and row of B.

I define

generated quantities{
  matrix[m,n*q] M;
  for (i in 1:m)
    for (j in 1:n)
      for (k in 1:p)
        for (l in 1:q)
          M[i, q*(j-1)+l] = to_matrix( row(A,i)) [i, j] * to_matrix( row(B, i))[k,l];
}

I have figured out a naive way of constructing a row-wise Kron product matrix.
In order to construct a row-wise Kron matrix, both matrices should have the same number of rows and thus through inspection, the relationship between a row-wise Kron and a full Kron matrix is as follows:

generated quantities{
//m is # of rows in matrix A
//n is # of cols in matrix A
//p is # of rows in matrix B
//q is # of rows in matrix B
  matrix[m*p,n*q] C; // full Kron product matrix
  matrix[m,n*q] M; //row-wise Kron product matrix 
  for (i in 1:m)
    for (j in 1:n)
      for (k in 1:p)
        for (l in 1:q)
          C[p*(i-1)+k,q*(j-1)+l] = A[i,j]*B[k,l];
     for(i in 1:m)
     M[i,] = C[m*(i-1)+i,];
}

However, currently I just experiment using a very small simulated data but my real dataset is quite large and I feel like it is quite redundant to first construct a full kronecker matrix and then extract desired rows from the matrix to construct the row-wise kronecker matrix.

Any ideas or advice that I can alternatively speed up the construction process?

Thanks in advance
Jaslene

2

I don’t really speak R or Kronecker well enough to put all these pieces together.

Is there a working Stan implementation of something you want to make faster? It looks like you have a function, but the way to optimize this is probably to make sure that however it’s used isn’t redundant. The point is to avoid building big intermediate structures if you don’t need them.

P.S. Can you edit the posts so that the code formats as code? I don’t know what’s going on with user posts, but a lot of them are coming through corrupted like this.

3

OK, I edited by removing empty vertical spaces and wrapping triple back-ticks around the code. There wasn’t any code markdown previously, which makes it really really hard for us to read code as otherwise this Discourse app just smashes everything together at the front of a line, except things that are exactly four spaces, which it treats as code, hence the original broken formatting.

4

Hi Bob, thanks for your reply. I have attached my R code and stan model below.
Essentially, I want to fit an additive model y=a0*x+f(x)+f(x,z)+iid errors.
where f(x) is constructed using 12 cubic b-spline basis functions
and f(x,z) is constructed through kronecker product of univariate b-splines of x and z, each consists of again 12 cubic b-spline basis functions. I want to construct b-spline basis functions within my stan code rather than using the functions outside. I just found the post https://github.com/milkha/Splines_in_Stan/blob/master/splines_in_stan.Rmd, where the author explains how to construct univariate b-splines and I try to mimic and construct a bivariate b-spline basis functions for my interaction term, as suggested in the GAM literature, that would be a Kron product of univariate b-splines.

Currently, I impose first-order difference penalty prior on f(x), I will need to impose the same penalty prior on f(x,z) but is yet to figure out how to do it properly.

Currently, using my method (quite redundant way), I verify the kroon product of univariate b-splines for f(x,z) to be the same as using the function kr in R. However, I would like to learn a more efficient way to construct it within stan code.
Any ideas or suggestions are greatly appreciated and welcome!
Thanks in advance,
Jaslene
b_spline_int_iid.stan (6.5 KB)

5

Just realised I can only upload one file at a time, here is the Rcode

gam model with penalty prior on univariate x and idd error.r (1.5 KB)

6

Hi Bob, I am wondering do you have time to take a look at the files I uploaded and perhaps give some suggestions?

Thank you so much for your time,
Jaslene

7

As I tried to say before, I’m not a linear algebra expert and not really who you want looking through these things. I’d have to start by looking up what a Kronecker product is!