power of a matrix

[pcdjo...@gmail.com]

Hi,

Can I raise a square matrix to a power? Not elementwise, but where the result of Q^n is is Q matrix-multiplied by itself n-1 times. I understand that pow() doesn't do this. I could use eigenvectors * eigenvalues^n * eigenvectors^-1, but eigenvalues_sym and eigenvectors_sym require a symmetric matrix and Q is not symmetric.

Thanks for any help,
Paul

[Ben Goodrich]

On Thursday, October 2, 2014 7:00:02 PM UTC-4, pcdjo...@gmail.com wrote:

Can I raise a square matrix to a power? Not elementwise, but where the result of Q^n is is Q matrix-multiplied by itself n-1 times. I understand that pow() doesn't do this. I could use eigenvectors * eigenvalues^n * eigenvectors^-1, but eigenvalues_sym and eigenvectors_sym require a symmetric matrix and Q is not symmetric.

No, but you can just loop it

for(i in 1:(n-1)) Q <- Q * Q;

 Ben

[Ben Goodrich]

I meant to say

for(i in 1:(n-1)) Q <- Q * original_Q;

Ben
 

[pcdjo...@gmail.com]

Thanks very much, that worked:

Q <- original_Q;
for(i in 1:(n-1)) Q <- Q * original_Q;

I hadn't thought of this before because it isn't possible in JAGS -- still learning the differences. My main motivation for changing from JAGS to Stan was that Stan appears to allow more matrix operations.

Paul

[Ben Goodrich]

On Thursday, October 9, 2014 9:04:01 AM UTC-4, pcdjo...@gmail.com wrote:

Thanks very much, that worked:

Q <- original_Q;
for(i in 1:(n-1)) Q <- Q * original_Q;

I hadn't thought of this before because it isn't possible in JAGS -- still learning the differences. My main motivation for changing from JAGS to Stan was that Stan appears to allow more matrix operations.

That is true, and as time goes on Stan will probably support even more matrix operations. However, the fundamental difference is HMC vs Gibbs-with-fallbacks and the best reason to switch is that Stan should have less dependence among its MCMC draws.

Ben

[Bob Carpenter]

JAGS could do this with an array to store intermediate
values:

Q[1] <- original_Q %*% original_Q;
for (n in 2:N)
Q[n] <- original_Q %*% Q[n-1];

Then you'd get your result in Q[N] = original_Q^N.
And a whole lot of intermediate outputs you would not want
to monitor.

We spend a fair bit of time suggesting that people don't
literally translate JAGS code like the above to Stan but
to use reassignment.

We could certainly add a matrix power operator. I added
a feature request, but did not assign it to anyone or make
it an immediate priority:

https://github.com/stan-dev/stan/issues/1073

I would suggest putting the loop in a function for the time
being.

functions {
matrix matrix_pow(matrix a, int n) {
if (n == 0)
return diag_matrix(rep_vector(0, rows(a)));
else if (n == 1)
return a;
else
return a * matrix_pow(a, n - 1);
}
}


- 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.
> For more options, visit https://groups.google.com/d/optout.

[pcdjo...@gmail.com]

I wasn't aware that I could use an array to do this in JAGS (I'm almost as new to JAGS as I am to Stan). Here is the JAGS code that worked, with the help of a local user (QW is the original 4x4 matrix, QWt[t, 1:4, 1:4]  is QW^t):

          QWt[1, 1:4, 1:4] <- QW
          for (gen in 2:t) { QWt[gen, 1:4, 1:4] <- QW %*% QWt[gen - 1, 1:4, 1:4] }

Many thanks for your help,
Paul