construct block matrix consisting of different powers of basis matrices

user_12158757 0
Viewed 128

Assuming that we have a basis matrix A

>> A = [1,2;-1,3]
A =

   1   2
  -1   3

I am wondering if there is an easy command that can construct a block matrix consisting of A, A^2, A^3 ... A^n with any given positive integer n, like

B = [A; A^2; A^3; ...; A^n]

Evidently, if A is a scalar, we can do like A.^(1:n). I was wondering if there is a similar manipulation or some command for matrix A. That's my motivation.

I can make it with loops (for loop or arrayfun) but I am looking for some easy and elegant command for it.

octave matrix matlab
5 Answers

Super verbose way of doing it in one line:

B= cell2mat(cellfun(@mpower,mat2cell(repmat(A,n,1),[ones(n,1)*size(A,1)],size(A,2)),mat2cell([1:n]',ones(n,1),1),'UniformOutput',false))

However, I think a for loop (with preallocation) is likely faster and much, much clearer. A for loop is not "inellegant" if its clear. Clarity should be your goal, unless its speed is super critical (and again, this one-liner still does loops under the hood).

I'm sure someone will come up with a clearer one liner, yet my comments will still apply.

PD: The magician @LuisMendo proposes: B=cell2mat(arrayfun(@(k){A^k},1:n).') as a much cleaner option. (Still loops inside)

user_14858720

Not necessarily an answer either, since it also relies on cellfun, and similar in spirit to other answers here, but I think this is slightly more elegant, and makes use of a built-in rather than an anonymous function (which is guaranteed to be faster, if that matters at all to you).

Also, this version is generalised for a matrix of exponents of any shape, producing the appropriate "block" output.

A = [ 1,2
     -1,3 ];

N = [ 1,2,3
      2,3,4
      3,4,5 ];

expblk = @(A,N) cell2mat(
             cellfun( @mpower, {A}, num2cell(N), 'UniformOutput', false )
         );

expblk(A,N)

% ans =
%     1    2   -1    8   -9   22
%    -1    3   -4    7  -11   13
%    -1    8   -9   22  -31   48
%    -4    7  -11   13  -24   17
%    -9   22  -31   48  -79   82
%   -11   13  -24   17  -41    3
user_41831910

In the case the matrix A is diagonalizable, you get an easy solution (not a one-liner, but still!) :

[V,D] = eig(A);
k = size(A,1);
B = kron(eye(n),V)*(reshape(D.^reshape(1:n,[1 1 n]),[k k*n]).')*inv(V);

Taken from the basic theory. Beside this solution in specific cases, I don't think that the general form involving Jordan's form is easy to implement in Matlab.

user_15070140

eval can be helpful:

A = [1,2;-1,3];
n = 3;
B = eval(['[A' sprintf(';A^%d',2:n) ']']);

However using a loop to multiply the matrix by itself is efficient specially when the matrix is large.

C = cell (n,1);
D = 1;
for k = 1:n
    D *= A; % In MATLAB use D = D * A;
    C{k} = D;
end
B = vertcat(C{:});

And its eval version in Octave syntax (for fun and not recommended):

B = eval(['[D=A' repmat(';D*=A',1,n-1) ']']);
user_65797440

If you don't mind using a recursive method, maybe the following could be an option

function y = f(A,n)
  if n==1 
    y = A;
    return;
  end 
  y = [A;f(A,n-1)*kron(eye(n-1),A)];
end

then

B =

    1    2
   -1    3
   -1    8
   -4    7
   -9   22
  -11   13
  -31   48
  -24   17
  -79   82
  -41    3
Admin1
Related