how to get an incremental power matrix in matlab

Question

I wanted to compute the following matrix in Matlab:

    g=[I
       A 
       .
       .
       .
      A^N]

I used the following program in Matlab:

A=[2 3;4 1];
s=A;
for n=1:1:50
   s(n)=A.^n;
end
g=[eye(1,1),s];

I am getting the following error:

In an assignment A(I) = B, the number of elements in B and I must be the same.
Error in s_x_calcu_v1 (line 5)
s(n)=A.^n;

you are mixing two operations in your question: `A^N` is **not** `A.^n` - do you want to compute matrix product or element-wise products? — Shai, May 11 '15 at 09:07
I wanted to do here A to the power N where n starts from 1 to 50,with 1 increment each time. — tina066, May 12 '15 at 09:44

Forss · Accepted Answer · 2015-05-11T08:04:42.877

The problem is that you are trying to assign a matrix to a single element. In matlab calling s(n) mean you get the nth element of s, regardless of the dimensions of s. You can use a three dimensional matrix

N = 50;
A=[2 3;4 1];

[nx,ny] = size(A);
s(nx,ny,N) = 0; %makes s a nx x ny x N matrix
for n=1:1:N
    s(:,:,n)=A.^n; %Colon to select all elements of that dimension
end
g=cat(3, eye(size(A)) ,s); %Add the I matrix of same size as A

Or a vectorized version

s = bsxfun(@power, A(:), 1:N);
s = reshape(s,2,2,N);
g = cat(3, eye(size(A)) ,s);

And a third solution using cumprod

s = repmat(A(:), [1 N]);
s = cumprod(s,2);
s = reshape(s,2,2,N);
g = cat(3, eye(size(A)) ,s);

1

`cumprod` saves computations! – Shai May 11 '15 at 09:10

score 3 · Answer 2 · answered May 11 '15 at 07:31

3

Your s array is a 2-by-2 array, you cannot index it to store the result of your compuation at each step of your loop.

For this, the simpler is probably to define s as a cell:

% --- Definitions
A = [2 3;4 1];
N = 50;

% --- Preparation
s = cell(N,1);

% --- Computation
for n=1:N
    s{n} = A.^n;
end

Best,

answered May 11 '15 at 07:31

Ratbert

5,463
2
18
37

I wanted to do here A to the power N where n starts from 1 to 50,with 1 increment each time. – tina066 May 12 '15 at 09:44

Shai · Answer 3 · 2015-05-12T09:49:16.070

When you loop from 1 to N computing each time A.^n you are doing LOTS of redundant computations! Note that

A.^n = (A.^(n-1)).*A; %//element-wise power
A^n = (A^n) * A; %// matrix power

Therefore,

A = [2 3;4 1];
N = 50;
s = cell(N+1,1);
s{1} = eye(size(A,1));
for ii=1:N
    s{ii+1} = s{ii}.*A; %// no powers, just product!
end
g = vertcat( s{:} );

BTW, the same holds if you want to compute matrix power (instead of element-wise powers), all you need is changing to s{ii+1} = s{ii}*A;

how to get an incremental power matrix in matlab

3 Answers3

Linked