Powers of a matrix

Question

I have a square matrix A (nxn). I would like to create a series of k powers of this matrix into an nxnxk multidimensional matrix (Not element-wise but actual powers of the matrix), i.e.getting [A^0 A^1 A^2..A^k]. It's sort of a varied vandermonde for matrix case.

I am able to do it with loops but it is annoying and slow. I tried using bsxfun but no luck since I am probably missing something here.

Here is a simple loop that I did:

for j=1:1:100 
    final(:,:,j)=A^(j-1); 
end

Answer 1

You are trying to perform cummulative version of mpower with a vector of k values.

Sadly, bsxfun hasn't evolved yet to handle such a case. So, the best I could suggest at this point would be having a running storage that accumulates the matrix-product at each iteration to be used at the next one.

Your original loop code looked something like this -

final = zeros([size(A),100]);
for j=1:1:100 
    final(:,:,j)=A^(j-1); 
end

So, with the suggestion, the modified loopy code would be -

final = zeros([size(A),100]);
matprod = A^0;
final(:,:,1) = matprod;
for j=2:1:100 
    matprod = A*matprod;
    final(:,:,j)= matprod;
end

Benchmarking -

%// Input
A = randi(9,200,200);

disp('---------- Original loop code -----------------')
tic
final = zeros([size(A),100]);
for j=1:1:100 
    final(:,:,j)=A^(j-1); 
end
toc

disp('---------- Modified loop code -----------------')
tic
final2 = zeros([size(A),100]);
matprod = A^0;
final2(:,:,1) = matprod;
for j=2:1:100 
    matprod = A*matprod;
    final2(:,:,j)= matprod;
end
toc

Runtimes -

---------- Original loop code -----------------
Elapsed time is 1.255266 seconds.
---------- Modified loop code -----------------
Elapsed time is 0.205227 seconds.