有没有办法在矩阵 repmat 阵列上做 mpower？

Question

我想知道是否有办法提高矩阵 A 作为数组的幂？

假设我们有这个矩阵

然后我们重复它的形状。

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

是否可以在 MATLAB/Octave 中不使用 for 循环来做到这一点？

Answer 1

编辑

在我的回答中也提到这一点：递归不是 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 2

另一个选项使用 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 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));

有没有办法在矩阵 repmat 阵列上做 mpower？

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

matlab

vectorization

编辑