A Fast and Efficient way to create a matrix from a series of product

Ax, Ay, Az: [N-by-N]

B=AA (a dyadic product)

It means :

B(i,j)= [Ax开发者_开发问答(i,j);Ay(i,j);Az(i,j)]*[Ax(i,j) Ay(i,j) Az(i,j)]

B(i,j) : a 3x3 matrix. One way to construct B is:

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

Is there a faster way for when N is large.

Edit:

Thanks for your answer. (It's faster) Let's put: N=2; Ax=[1 2;3 4]; Ay=[5 6;7 8]; Az=[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

Run:

??? Error using ==> mtimes Inner matrix dimensions must agree.

If I write :P = Ai*Aj; then

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

That is defferent from above A(:,:,1) deffer from [Ax(1,1) Ay(1,1) Az(1,1)]

Edit:

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.

Edit:

After some modification for my application : by gnovice codes

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

It seems that any function calling like ceil,ind2sub ... make thw loops slow and shoud avoid if possible.

symIndex was interesting! Thank you.

---

Here's a fairly simple and general implementation that uses a single for loop to perform linear indexing and avoids dealing with 3-dimensional variables or reshaping:

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

EDIT #1:

In response to the additional question of how to reduce the number of calculations performed for a symmetric matrix B, you could use the following modified version of the above code:

%# 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)];
  B{index} = 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
    B{symIndex} = B{index};  %#   then copy the symmetric element
  end
end
B = cell2mat(B);

However, in such a case you may get better performance (or at least simpler looking code) by foregoing the linear indexing and using two for loops with subscripted indexing like so:

%# 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)];
    B{r,c} = A(:)*A;
    if r ~= c           %# If we're not on the main diagonal
      B{c,r} = B{r,c};  %#   then copy the symmetric element
    end
  end
end
B = cell2mat(B);

EDIT #2:

The second symmetric solution can be further sped up by moving the diagonal calculation outside of the inner loop (removing the need for a conditional statement) and overwriting A with the result A(:)*A so that we can update the symmetric cell entry B{c,r} using A instead of B{r,c} (i.e. updating one cell with the contents of another appears to carry some extra overhead):

%# Symmetric solution #3:
%# ---------------------
B = cell(N);
for c = 1:N    %# Loop over the columns
  A = [Ax(c,c) Ay(c,c) Az(c,c)];
  B{c,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;
    B{r,c} = A;
    B{c,r} = A;
  end
end
B = cell2mat(B);

And here are some timing results for the 3 symmetric solutions using the following sample symmetric matrices Ax, Ay, and Az:

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

The big speed-up for solution 3 results primarily from updating the cell contents of B with the local variable A instead of updating one cell with the contents of another. In other words, copying cell contents with operations like B{c,r} = B{r,c}; carries more overhead than I expected.

---

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