Есть ли способ сделать mpower на матрицу repmat массив?

Question

Интересно, есть ли способ повысить мощность матрицы A как массива?

Предположим, что у нас есть эта матрица

A =

   5   4
   3   6

Затем мы повторимего форма.

>> 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

Теперь я хочу увеличить мощность, чтобы длинная повторяющаяся матрица выглядела так:

>> [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

Возможно ли это сделать без цикла for в MATLAB /Октав

Answer 1

Другой вариант с использованием arrayfun

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

Результат:

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

Но основной цикл for, вероятно, будет самым быстрым вариантом в этом случае.

Бенчмаркинг с октавой:

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

Истекшее время составляет 1,023 секунды.

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

Истекшее время составляет 4,8684 секунды.

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

Истекшее время составляет 0,039371 секунды

Answer 2

EDIT

Также упомянуть об этом в моем ответе: рекурсия - это не векторизация, как обычно думают пользователи Matlab / Octave. У меня просто возникла мысль о рекурсивной анонимной функции, и я обнаружил, что данная задача - хороший небольшой пример для тестирования ссылочного решения.

---

Я искал рекурсивные анонимные функции и нашел отличный ответ . Я включил идеи оттуда, чтобы оправдать ожидания, изложенные в вопросе, и пришел к этому короткому фрагменту кода.

% 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)

Для объяснения, пожалуйста, посмотрите на связанный ответ.

Выход:

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

Answer 3

Если вы действительно хотите использовать векторизацию (в данном случае это IMO-избыточность), вы также можете использовать свойство:

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

Где

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

Итак

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);

С:

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

Вам все еще нужно извлечь эти значения. В этом случае было бы интересно использовать разреженную матрицу для уменьшения использования памяти.

Примерно так:

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