Vectorization of matrices power

Hi all,

Let A a matrix n*n and let N an integer. I wish to create a matrix of power of the form AN = [I A A^2,...,A^N], where I=eye(n), without resorting to a "for" cycle. I have tried to take a look at several commands, including cumprod, kron, etc, and trying to combine them, but I failed. I made it only for the scalar case, i.e. n=1.

After that, I wish to create a matrix

AA = [I 0 0 0;A I 0 0;A^2 A I 0;A^3 A^2 A I]

without using any "for" cycle. I have noticed that the first column of AA is equal to AN' (if it may help). Any hints? Thanks.

[Matt J]

It is doubtful that a for-loop is to be feared here. Surely you can't be doing this for N very large?

Answer 1

[Matt J]

>> [AN,T]=matpowers(diag([2,3]),3),

AN =

       1     0     2     0     4     0     8     0
       0     1     0     3     0     9     0    27

T =

       1     0     0     0     0     0     0     0
       0     1     0     0     0     0     0     0
       2     0     1     0     0     0     0     0
       0     3     0     1     0     0     0     0
       4     0     2     0     1     0     0     0
       0     9     0     3     0     1     0     0
       8     0     4     0     2     0     1     0
       0    27     0     9     0     3     0     1

function  [AN,T]=matpowers(A,c)

   persistent pcell N

   if isempty(pcell), pcell={eye(size(A))}; end
   if nargin>1, N=c; end

   if ~isempty(N) && N>0

       pcell=[pcell,{pcell{end}*A}];

       N=N-1;

      AN=matpowers(A); 
   else

       AN=pcell; pcell=[]; N=[];

   end

   if nargout>1
    z=[AN,{zeros(size(A))}];
    T=toeplitz(1:c+1,[1,ones(1,c)*(c+2)]);
    T=cell2mat(z(T)) ;
    AN=cell2mat(AN);

   end

Answer 2

[Sean de Wolski]

This should give you the tools you need:

x = (1:3).^(1:3)
xm = tril(toeplitz(x))

Answer 3

[Ubaldo Tiberi]

Thanks for the reply.

Sean, unfortunately it seems that your solution works only with scalar values and not for matrices.

Matt, thanks for the impressive answer! However, I was wondering if I can get the same result without using a function, but just by combining elementary Matlab operations such as blkdiag, kron, cumprod, etc. But I do like this function! However, how would you write a "for" cycle that solves my problem?

Thanks again!

[Matt J]

but just by combining elementary Matlab operations such as blkdiag, kron, cumprod

If A is symmetric you might be able to. But I think a for-loop will be the most efficient, regardless. Are you sure you need to build these matrices explicitly? They contain a lot of redundant data. What are you planning to use them for?