一种从一系列产品中创建矩阵的快速有效的方法

Posted 2023-03-11

【中文标题】一种从一系列产品中创建矩阵的快速有效的方法

【英文标题】：A Fast and Efficient way to create a matrix from a series of product

【发布时间】：2011-09-08 14:35:49

【问题描述】：

Ax、Ay、Az：[N×N]

B=AA（二元乘积）

意思是：

B(i,j)= [Ax(i,j);Ay(i,j);Az(i,j)]*[Ax(i,j) Ay(i,j) Az(i,j)]

B(i,j) ：一个 3x3 矩阵。 构造 B 的一种方法是：

N=2;
Ax=rand(N); Ay=rand(N); Az=rand(N);     %# [N-by-N]
t=1;
F=zeros(3,3,N^2);
for  i=1:N
    for j=1:N
        F(:,:,t)= [Ax(i,j);Ay(i,j);Az(i,j)]*[Ax(i,j) Ay(i,j) Az(i,j)];
        t=t+1;                          %# t is just a counter
    end
end
%# then we can write
B = mat2cell(F,3,3,ones(N^2,1));
B = reshape(B,N,N)'; 
B = cell2mat(B);

当N很大时有没有更快的方法。

编辑：

感谢您的回答。 （更快） 让我们说： N=2； 斧头=[1 2;3 4]; Ay=[5 6;7 8];阿兹=[9 10;11 12];

B =

 1     5     9     4    12    20
 5    25    45    12    36    60
 9    45    81    20    60   100
 9    21    33    16    32    48
21    49    77    32    64    96
33    77   121    48    96   144

运行： ？？？使用 ==> 时出错 内部矩阵尺寸必须一致。

如果我写：P = Ai*Aj; 那么

B  =

 7    19    31    15    43    71
23    67   111    31    91   151
39   115   191    47   139   231
10    22    34    22    50    78
34    78   122    46   106   166
58   134   210    70   162   254

这与上面不同 A(:,:,1) 不同于 [Ax(1,1) Ay(1,1) Az(1,1)]

编辑：

N=100;
Me      :Elapsed time is 1.614244 seconds.
gnovice :Elapsed time is 0.056575 seconds.
N=200;
Me      :Elapsed time is 6.044628 seconds.
gnovice :Elapsed time is 0.182455 seconds.
N=400;
Me      :Elapsed time is 23.775540 seconds.
gnovice :Elapsed time is 0.756682 seconds.
Fast!

rwong: B was not the same.

编辑：

对我的应用程序进行一些修改后： 通过新手代码

 1st code :  19.303310 seconds
 2nd code:  23.128920  seconds
 3rd code:  13.363585  seconds

似乎任何像 ceil,ind2sub ... 这样的函数调用都会使 thw 循环变慢，如果可能的话应该避免。

symIndex 很有趣！谢谢。

【参考方案1】：

这是一个相当简单和通用的实现，它使用单个 for 循环来执行 linear indexing 并避免处理 3 维变量或重塑：

%# General solution:
%# ----------------
B = cell(N);
for index = 1:N^2
  A = [Ax(index) Ay(index) Az(index)];
  Bindex = A(:)*A;
end
B = cell2mat(B);

编辑#1：

针对如何减少对对称矩阵B 执行的计算次数的附加问题，您可以使用上述代码的以下修改版本：

%# Symmetric solution #1:
%# ---------------------
B = cell(N);
for index = find(tril(ones(N))).'  %'# Loop over the lower triangular part of B
  A = [Ax(index) Ay(index) Az(index)];
  Bindex = A(:)*A;
  symIndex = N*rem(index-1,N)+ceil(index/N);  %# Find the linear index for the
                                              %#   symmetric element
  if symIndex ~= index       %# If we're not on the main diagonal
    BsymIndex = Bindex;  %#   then copy the symmetric element
  end
end
B = cell2mat(B);

但是，在这种情况下，您可以通过前面的线性索引并使用两个带有下标索引的 for 循环来获得更好的性能（或者至少看起来更简单的代码），如下所示：

%# Symmetric solution #2:
%# ---------------------
B = cell(N);
for c = 1:N    %# Loop over the columns
  for r = c:N  %# Loop over a subset of the rows
    A = [Ax(r,c) Ay(r,c) Az(r,c)];
    Br,c = A(:)*A;
    if r ~= c           %# If we're not on the main diagonal
      Bc,r = Br,c;  %#   then copy the symmetric element
    end
  end
end
B = cell2mat(B);

编辑 #2：

第二个对称解决方案可以通过将对角线计算移到内部循环之外（消除对条件语句的需要）并用结果A(:)*A 覆盖A 来进一步加快速度，这样我们就可以更新对称单元格条目Bc,r 使用A 而不是Br,c（即用另一个单元格的内容更新一个单元格似乎会带来一些额外的开销）：

%# Symmetric solution #3:
%# ---------------------
B = cell(N);
for c = 1:N    %# Loop over the columns
  A = [Ax(c,c) Ay(c,c) Az(c,c)];
  Bc,c = A(:)*A;
  for r = c+1:N  %# Loop over a subset of the rows
    A = [Ax(r,c) Ay(r,c) Az(r,c)];
    A = A(:)*A;
    Br,c = A;
    Bc,r = A;
  end
end
B = cell2mat(B);

以下是使用以下示例对称矩阵 Ax、Ay 和 Az 的 3 个对称解决方案的一些时序结果：

N = 400;
Ax = tril(rand(N));     %# Lower triangular matrix
Ax = Ax+triu(Ax.',1);  %'# Add transpose to fill upper triangle
Ay = tril(rand(N));     %# Lower triangular matrix
Ay = Ay+triu(Ay.',1);  %'# Add transpose to fill upper triangle
Az = tril(rand(N));     %# Lower triangular matrix
Az = Az+triu(Az.',1);  %'# Add transpose to fill upper triangle

%# Timing results:
%# --------------
%# Solution #1 = 0.779415 seconds
%# Solution #2 = 0.704830 seconds
%# Solution #3 = 0.325920 seconds

解决方案 3 的大幅加速主要是因为使用局部变量 A 更新 B 的单元格内容，而不是使用另一个单元格的内容更新一个单元格。换句话说，使用Bc,r = Br,c; 之类的操作复制单元格内容的开销超出了我的预期。

【讨论】：

非常感谢 gnovice。这太棒了！是否有可能在线性索引中我们对算法说当 B 对称时只计算矩阵的一半？

【参考方案2】：

A = cat(3, Ax, Ay, Az);   % [N-by-N-by-3]
F = zeros(3, 3, N^2);
for i = 1:3, 
  for j = 1:3,
    Ai = A(:,:,i);
    Aj = A(:,:,j);
    P = Ai(:) .* Aj(:);
    F(i,j,:) = reshape(P, [1, 1, N^2]);
  end
end

%# then we can write
B = mat2cell(F,3,3,ones(N^2,1));
B = reshape(B,N,N)'; 
B = cell2mat(B);

【讨论】：

P = Ai(:) * Aj(:); --> P = Ai * Aj; ?

@user784433：已修复。它应该按元素相乘，得到一个[N^2, 1, 1] 向量。谢谢。