Is there a way to do mpower on a matrix repmat array?

Question

I wonder if there is a way raise power of a matrix A as an array?

Assume that we have this matrix

Then we repeat its shape.

>> repmat(A, 5, 1)
ans =

   5   4
   3   6
   5   4
   3   6
   5   4
   3   6
   5   4
   3   6
   5   4
   3   6

Now I want to rase the power so the long repeated matrix looks like this:

>> [A^1; A^2; A^3; A^4; A^5]
ans =

       5       4
       3       6
      37      44
      33      48
     317     412
     309     420
    2821    3740
    2805    3756
   25325   33724
   25293   33756

Is it possible to do that without a for loop in MATLAB/Octave?

MATLAB or Octave? Please pick one. Answers might not be the same. Notice the tag usage guideline: "Don’t use both the [matlab] and [octave] tags, unless the question is explicitly about the similarities or differences between the two." — Cris Luengo, Mar 25 '19 at 22:29
Also, why do you want to avoid a loop? Is this part of your code a bottleneck? — Cris Luengo, Mar 25 '19 at 22:30
@CrisLuengo It does not matter. Matlab and Octave have the same syntax and Octave can do MATLAB code. — euraad, Mar 25 '19 at 22:35
For loops are not at all slow any more in MATLAB, and many vectorization efforts actually slow down execution. Don't vectorize until you've programmed the simple solution and seen it is too slow for your application. Otherwise you're just wasting time. And you need to be able to compare the timing of your vectorized code, to make sure you are actually picking the fastest code. — Cris Luengo, Mar 25 '19 at 23:04
Octave has similar syntax to MATLAB, but not identical. Octave also lacks a lot of the functions that are available in MATLAB, and it has some functions that do not exist in MATLAB. Octave cannot do a lot of MATLAB code I have written, and MATLAB cannot do a lot of Octave code people have posted in their questions here. — Cris Luengo, Mar 25 '19 at 23:05
@CrisLuengo Not so much. Please update https://www.gnu.org/software/octave/NEWS-5.1.html — euraad, Mar 26 '19 at 07:35

score 3 · Answer 1 · edited Jun 20 '20 at 09:12

EDIT

To also mention this in my answer: Recursion is not vectorization as of Matlab/Octave users usually might think of. I just had the idea of a recursive, anonymous function in my mind, and found the given task a nice small example to test the referenced solution on.

I was looking for recursive, anonymous functions and found this great answer. I incorporated the ideas from there to meet the expectations stated in the question, and came to this short code snippet.

% Input.
A = [5 4; 3 6]

% Set up recursive anonymous function.
iif = @(varargin) varargin{2*find([varargin{1:2:end}], 1, 'first')}();
recPower = @(A, B, n, f) iif(n > 1, @() [B; f(A, A * B, n - 1, f)], true, @() B);
nPower = @(A, n) recPower(A, A, n, recPower);

% Calculate for arbitrary n.
nPower(A, 5)

For explanations, please have a look at the linked answer.

Output:

A =

   5   4
   3   6

ans =

       5       4
       3       6
      37      44
      33      48
     317     412
     309     420
    2821    3740
    2805    3756
   25325   33724
   25293   33756

@SardarUsama I know, and actually I wanted to mention it in my answer. Generally, I don't think there is a "plain" vectorization possibility in terms of Matlab/Octave users usually think of. Nevertheless, you're right. — HansHirse, Mar 26 '19 at 10:10

obchardon · Accepted Answer · 2019-03-26T09:57:00.703

Another option using arrayfun

B = cell2mat(arrayfun(@(x)A^x,1:5,'UniformOutput',0).')

Result:

B =
   5       4
   3       6
  37      44
  33      48
 317     412
 309     420
2821    3740
2805    3756
25325   33724
25293   33756

But a basic for-loop would probably be the fastest option in this case.

Benchmarking with octave:

tic
iif = @(varargin) varargin{2*find([varargin{1:2:end}], 1, 'first')}();
recPower = @(A, B, n, f) iif(n > 1, @() [B; f(A, A * B, n - 1, f)], true, @() B);
nPower = @(A, n) recPower(A, A, n, recPower);
for ii = 1:1000
% Calculate for arbitrary n.
    nPower(A, 5);
end
toc

Elapsed time is 1.023 seconds.

tic
for ii = 1:1000
    B = cell2mat(arrayfun(@(x)A^x,1:5,'UniformOutput',0).');
end
toc

Elapsed time is 4.8684 seconds.

tic
for ii = 1:1000
    B=[];
    for jj = 1:5
    B = [B;A^jj];
    end
end
toc

Elapsed time is 0.039371 seconds

`arrayfun` is a loop. The tags indicate that vectorization is required — Sardar Usama, Mar 26 '19 at 09:56
Using `arrayfun` was also my first idea. Principally, I agree on using loops for that task, as also initially stated by Cris Luengo in the above comments. — HansHirse, Mar 26 '19 at 10:15

obchardon · Answer 3 · 2019-03-26T13:10:58.883

If you really want to use vectorization (which is IMO overkill in this case), you could also use the property:

A^n = P*D^n*P^-1 %A SHOULD BE a diagonalizable matrix

Where

[P,D] = eig(A) %the eigenvectors and eigenvalue

So

A = [5 4; 3 6]
n = 5;
% get the eigenvalue/eigenvector
[P,D]=eig(A);
% create the intermediate matrix
MD = diag(D).^[1:n];
MD = diag(MD(:));
% get the result
SP = kron(eye(n,n),P)*MD*kron(eye(n,n),P^-1);

With:

SP =

      5      4      0      0      0      0      0      0      0      0
      3      6      0      0      0      0      0      0      0      0
      0      0     37     44      0      0      0      0      0      0
      0      0     33     48      0      0      0      0      0      0
      0      0      0      0    317    412      0      0      0      0
      0      0      0      0    309    420      0      0      0      0
      0      0      0      0      0      0   2821   3740      0      0
      0      0      0      0      0      0   2805   3756      0      0
      0      0      0      0      0      0      0      0  25325  33724
      0      0      0      0      0      0      0      0  25293  33756

You still just need to extract those values. It could be interesting to use sparse matrix in this case to reduce memory usage.

Something like this:

SP = sparse(kron(eye(n,n),P))*sparse(MD)*sparse(kron(eye(n,n),P^-1));

You should add this to your timing in the other answer. – Cris Luengo Mar 26 '19 at 12:27 — Cris Luengo, Mar 26 '19 at 12:27

Is there a way to do mpower on a matrix repmat array?

3 Answers3

EDIT