construct a row-wise kronecker product matrix

182 views
Skip to first unread message

Jaslene Lin

unread,
May 25, 2017, 1:47:20 AM5/25/17
to Stan users mailing list

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];
          
         

}

However, it does not work. I would like some advice or any ideas would be greatly appreciated.

Thanks in advance.

Jaslene

Jaslene Lin

unread,
May 25, 2017, 1:56:28 AM5/25/17
to Stan users mailing list

In MGLM package, there is a function kr( A, B) that would automatically produce a row-wise Kron matrix however, for my model, I do need to build in the row-wise Kron matrix within my transformed data block so I guess I have to do it manually as Stan does not have built-in function for Kron matrix at the moment. 


Jaslene Lin

unread,
May 25, 2017, 7:21:53 AM5/25/17
to Stan users mailing list
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,];

}


This is probably not the smartest and the most ideal way, I open to any suggestions and advice.
Thanks in advance:)

Jaslene

Bob Carpenter

unread,
May 25, 2017, 4:39:03 PM5/25/17
to stan-...@googlegroups.com
I'm not an expert here, but there's usually so much duplication in
these Kronecker products that you don't want to just blow them out
into a dense matrix, but rather do specialized operations with the
full matrix remaining implicit in the calculations.

- Bob
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

Jaslene Lin

unread,
May 25, 2017, 5:00:10 PM5/25/17
to Stan users mailing list
Hi Bob:

Thanks for your reply. Yes, currently I am just experimenting with a simulated example, which is a very small dataset. Later, I will need to use this in a much larger dataset and I am afraid of the running time. Are there any examples that I can look at in terms of doing specialised operations with the full matrix remaining the others implicit in calculation? I am new to this topic as well and would like to learn some tips and tricks. 

Thanks again, 
Jaslene


Bob Carpenter

unread,
May 25, 2017, 5:09:00 PM5/25/17
to stan-...@googlegroups.com
Someone else will probably jump in here (if not, try our
new Discourse list at: http://discourse.mc-stan.org

What you want to do in general is look at the operations.
If that Kronecker is multiplied by a vector, then the
result is much smaller and you don't need to expand out
the Kronecker product---just figure out how to get the right
vector as a result.

- Bob

Jaslene Lin

unread,
May 25, 2017, 5:19:34 PM5/25/17
to Stan users mailing list
Thanks Bob, yes perhaps I think look into the Kronecker matrix theory as currently it is quite redundant to construct a full kronecker matrix and then extract from it. 
Reply all
Reply to author
Forward
0 new messages