Fill a matrix with matrix powers

7 views (last 30 days)
Luigi Emanuel di Grazia on 22 Jul 2021
Answered: Steven Lord on 22 Jul 2021
Hi to everyone,
I was wondering if anyone knows the fastest way to achieve the following:
Given A [n x n], fill a matrix B such as:
B= [A 0n ... 0n;
0n A^2 ... 0n;
.... ;
0n 0n ... A^n]
where 0n=zeros(n).

Rik on 22 Jul 2021
n=3;
A=rand(n,n);
Zero=zeros(size(A));
C=repmat({Zero},n,n);
C(logical(eye(n)))=arrayfun(@(n)A^n,1:n,'uni',false);
C=cell2mat(C)
C = 9×9
0.6504 0.2955 0.1950 0 0 0 0 0 0 0.4001 0.8820 0.5801 0 0 0 0 0 0 0.7784 0.0231 0.7369 0 0 0 0 0 0 0 0 0 0.6931 0.4574 0.4420 0 0 0 0 0 0 1.0647 0.9096 1.0171 0 0 0 0 0 0 1.0892 0.2675 0.7082 0 0 0 0 0 0 0 0 0 0.9779 0.6185 0.7262 0 0 0 0 0 0 1.8482 1.1404 1.4847 0 0 0 0 0 0 1.3668 0.5742 0.8895

Steven Lord on 22 Jul 2021
A = magic(3);
AM = {A^0, A^1, A^2};
celldisp(AM)
AM{1} = 1 0 0 0 1 0 0 0 1 AM{2} = 8 1 6 3 5 7 4 9 2 AM{3} = 91 67 67 67 91 67 67 67 91
B = blkdiag(AM{:})
B = 9×9
1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 8 1 6 0 0 0 0 0 0 3 5 7 0 0 0 0 0 0 4 9 2 0 0 0 0 0 0 0 0 0 91 67 67 0 0 0 0 0 0 67 91 67 0 0 0 0 0 0 67 67 91
You could create AM automatically rather than hard-coding it if you wanted a larger B.
AM2 = arrayfun(@(x) A^x, 0:2, 'UniformOutput', false);
check = isequal(AM, AM2)
check = logical
1
C = blkdiag(AM2{:});
isequal(B, C)
ans = logical
1