Vectorization of matrices power
1 view (last 30 days)
Show older comments
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.
1 Comment
Matt J
on 12 Nov 2013
It is doubtful that a for-loop is to be feared here. Surely you can't be doing this for N very large?
Accepted Answer
Matt J
on 12 Nov 2013
>> [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
0 Comments
More Answers (2)
Sean de Wolski
on 12 Nov 2013
This should give you the tools you need:
x = (1:3).^(1:3)
xm = tril(toeplitz(x))
0 Comments
Ubaldo Tiberi
on 13 Nov 2013
1 Comment
Matt J
on 13 Nov 2013
Edited: Matt J
on 13 Nov 2013
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?
See Also
Categories
Find more on Matrices and Arrays in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!