English 中文(简体)
Multiply a 3D matrix with a 2D matrix
原标题:

Suppose I have an AxBxC matrix X and a BxD matrix Y.

Is there a non-loop method by which I can multiply each of the C AxB matrices with Y?

最佳回答

You can do this in one line using the functions NUM2CELL to break the matrix X into a cell array and CELLFUN to operate across the cells:

Z = cellfun(@(x) x*Y,num2cell(X,[1 2]), UniformOutput ,false);

The result Z is a 1-by-C cell array where each cell contains an A-by-D matrix. If you want Z to be an A-by-D-by-C matrix, you can use the CAT function:

Z = cat(3,Z{:});



NOTE: My old solution used MAT2CELL instead of NUM2CELL, which wasn t as succinct:

[A,B,C] = size(X);
Z = cellfun(@(x) x*Y,mat2cell(X,A,B,ones(1,C)), UniformOutput ,false);
问题回答

As a personal preference, I like my code to be as succinct and readable as possible.

Here s what I would have done, though it doesn t meet your no-loops requirement:

for m = 1:C

    Z(:,:,m) = X(:,:,m)*Y;

end

This results in an A x D x C matrix Z.

And of course, you can always pre-allocate Z to speed things up by using Z = zeros(A,D,C);.

Here s a one-line solution (two if you want to split into 3rd dimension):

A = 2;
B = 3;
C = 4;
D = 5;

X = rand(A,B,C);
Y = rand(B,D);

%# calculate result in one big matrix
Z = reshape(reshape(permute(X, [2 1 3]), [A B*C]), [B A*C])  * Y;

% # split into third dimension
Z = permute(reshape(Z ,[D A C]),[2 1 3]);

Hence now: Z(:,:,i) contains the result of X(:,:,i) * Y


Explanation:

The above may look confusing, but the idea is simple. First I start by take the third dimension of X and do a vertical concatenation along the first dim:

XX = cat(1, X(:,:,1), X(:,:,2), ..., X(:,:,C))

... the difficulty was that C is a variable, hence you can t generalize that expression using cat or vertcat. Next we multiply this by Y:

ZZ = XX * Y;

Finally I split it back into the third dimension:

Z(:,:,1) = ZZ(1:2, :);
Z(:,:,2) = ZZ(3:4, :);
Z(:,:,3) = ZZ(5:6, :);
Z(:,:,4) = ZZ(7:8, :);

So you can see it only requires one matrix multiplication, but you have to reshape the matrix before and after.

I m approaching the exact same issue, with an eye for the most efficient method. There are roughly three approaches that i see around, short of using outside libraries (i.e., mtimesx):

  1. Loop through slices of the 3D matrix
  2. repmat-and-permute wizardry
  3. cellfun multiplication

I recently compared all three methods to see which was quickest. My intuition was that (2) would be the winner. Here s the code:

% generate data
A = 20;
B = 30;
C = 40;
D = 50;

X = rand(A,B,C);
Y = rand(B,D);

% ------ Approach 1: Loop (via @Zaid)
tic
Z1 = zeros(A,D,C);
for m = 1:C
    Z1(:,:,m) = X(:,:,m)*Y;
end
toc

% ------ Approach 2: Reshape+Permute (via @Amro)
tic
Z2 = reshape(reshape(permute(X, [2 1 3]), [A B*C]), [B A*C])  * Y;
Z2 = permute(reshape(Z2 ,[D A C]),[2 1 3]);
toc


% ------ Approach 3: cellfun (via @gnovice)
tic
Z3 = cellfun(@(x) x*Y,num2cell(X,[1 2]), UniformOutput ,false);
Z3 = cat(3,Z3{:});
toc

All three approaches produced the same output (phew!), but, surprisingly, the loop was the fastest:

Elapsed time is 0.000418 seconds.
Elapsed time is 0.000887 seconds.
Elapsed time is 0.001841 seconds.

Note that the times can vary quite a lot from one trial to another, and sometimes (2) comes out the slowest. These differences become more dramatic with larger data. But with much bigger data, (3) beats (2). The loop method is still best.

% pretty big data...
A = 200;
B = 300;
C = 400;
D = 500;
Elapsed time is 0.373831 seconds.
Elapsed time is 0.638041 seconds.
Elapsed time is 0.724581 seconds.

% even bigger....
A = 200;
B = 200;
C = 400;
D = 5000;
Elapsed time is 4.314076 seconds.
Elapsed time is 11.553289 seconds.
Elapsed time is 5.233725 seconds.

But the loop method can be slower than (2), if the looped dimension is much larger than the others.

A = 2;
B = 3;
C = 400000;
D = 5;
Elapsed time is 0.780933 seconds.
Elapsed time is 0.073189 seconds.
Elapsed time is 2.590697 seconds.

So (2) wins by a big factor, in this (maybe extreme) case. There may not be an approach that is optimal in all cases, but the loop is still pretty good, and best in many cases. It is also best in terms of readability. Loop away!

Nope. There are several ways, but it always comes out in a loop, direct or indirect.

Just to please my curiosity, why would you want that anyway ?

To answer the question, and for readability, please see:

  • ndmult, by ajuanpi (Juan Pablo Carbajal), 2013, GNU GPL

Input

  • 2 arrays
  • dim

Example

 nT = 100;
 t = 2*pi*linspace (0,1,nT)’;

 # 2 experiments measuring 3 signals at nT timestamps
 signals = zeros(nT,3,2);
 signals(:,:,1) = [sin(2*t) cos(2*t) sin(4*t).^2];
 signals(:,:,2) = [sin(2*t+pi/4) cos(2*t+pi/4) sin(4*t+pi/6).^2];

 sT(:,:,1) = signals(:,:,1)’;
 sT(:,:,2) = signals(:,:,2)’;
   G = ndmult (signals,sT,[1 2]);

Source

Original source. I added inline comments.

function M = ndmult (A,B,dim)
  dA = dim(1);
  dB = dim(2);

  # reshape A into 2d
  sA = size (A);
  nA = length (sA);
  perA = [1:(dA-1) (dA+1):(nA-1) nA dA](1:nA);
  Ap = permute (A, perA);
  Ap = reshape (Ap, prod (sA(perA(1:end-1))), sA(perA(end)));

  # reshape B into 2d
  sB = size (B);
  nB = length (sB);
  perB = [dB 1:(dB-1) (dB+1):(nB-1) nB](1:nB);
  Bp = permute (B, perB);
  Bp = reshape (Bp, sB(perB(1)), prod (sB(perB(2:end))));

  # multiply
  M = Ap * Bp;

  # reshape back to original format
  s = [sA(perA(1:end-1)) sB(perB(2:end))];
  M = squeeze (reshape (M, s));
endfunction

I highly recommend you use the MMX toolbox of matlab. It can multiply n-dimensional matrices as fast as possible.

The advantages of MMX are:

  1. It is easy to use.
  2. Multiply n-dimensional matrices (actually it can multiply arrays of 2-D matrices)
  3. It performs other matrix operations (transpose, Quadratic Multiply, Chol decomposition and more)
  4. It uses C compiler and multi-thread computation for speed up.

For this problem, you just need to write this command:

C=mmx( mul ,X,Y);

here is a benchmark for all possible methods. For more detail refer to this question.

    1.6571 # FOR-loop
    4.3110 # ARRAYFUN
    3.3731 # NUM2CELL/FOR-loop/CELL2MAT
    2.9820 # NUM2CELL/CELLFUN/CELL2MAT
    0.0244 # Loop Unrolling
    0.0221 # MMX toolbox  <===================

I would like to share my answer to the problems of:

1) making the tensor product of two tensors (of any valence);

2) making the contraction of two tensors along any dimension.

Here are my subroutines for the first and second tasks:

1) tensor product:

function [C] = tensor(A,B)
   C = squeeze( reshape( repmat(A(:), 1, numel(B)).*B(:).  , [size(A),size(B)] ) );
end

2) contraction: Here A and B are the tensors to be contracted along the dimesions i and j respectively. The lengths of these dimensions should be equal, of course. There s no check for this (this would obscure the code) but apart from this it works well.

   function [C] = tensorcontraction(A,B, i,j)
      sa = size(A);
      La = length(sa);
      ia = 1:La;
      ia(i) = [];
      ia = [ia i];

      sb = size(B);
      Lb = length(sb);
      ib = 1:Lb;
      ib(j) = [];
      ib = [j ib];

      % making the i-th dimension the last in A
      A1 = permute(A, ia);
      % making the j-th dimension the first in B
      B1 = permute(B, ib);

      % making both A and B 2D-matrices to make use of the
      % matrix multiplication along the second dimension of A
      % and the first dimension of B
      A2 = reshape(A1, [],sa(i));
      B2 = reshape(B1, sb(j),[]);

      % here s the implicit implication that sa(i) == sb(j),
      % otherwise - crash
      C2 = A2*B2;

      % back to the original shape with the exception
      % of dimensions along which we ve just contracted
      sa(i) = [];
      sb(j) = [];
      C = squeeze( reshape( C2, [sa,sb] ) );
   end

Any critics?

I would think recursion, but that s the only other non- loop method you can do

You could "unroll" the loop, ie write out all the multiplications sequentially that would occur in the loop





相关问题
Matrix to Represent a Triangle in Screen Space

So i have a set of four points in 3D Space. P1 [0, 0, 0] P2 [128, 0, 0] P3 [0, 128, 0] P4 [128, 128, 0] Which I m then projecting orthographically to the screen effectively giving me two ...

Multiply a 3D matrix with a 2D matrix

Suppose I have an AxBxC matrix X and a BxD matrix Y. Is there a non-loop method by which I can multiply each of the C AxB matrices with Y?

matrix and vector template classes in c++

#include <array> template <typename T> class Vector4<T> { std::array<T, 4> _a; // or T _a[4]; ? }; template <typename T> class Matrix4<T> { std::array<...

Linear Independence Matrix

Suppose we have a m by n matrix A with rank m and a set K⊆{1..n} such that the columns of A indexed by K are linearly independent. Now we want to extend K and find a set L so that k⊆L and columns ...

Difference between MATLAB s matrix notations

How do you read the following MATLAB codes? #1 K>> [p,d]=eig(A) // Not sure about the syntax. p = 0.5257 -0.8507 -0.8507 -0.5257 d = ...

Strange error when using sparse matrices and glmnet

I m getting a weird error when training a glmnet regression. invalid class "dgCMatrix" object: length(Dimnames[[2]]) must match Dim[2] It only happens occasionally, and perhaps only under larger ...

Java large datastructure for storing a matrix

I need to store a 2d matrix containing zip codes and the distance in km between each one of them. My client has an application that calculates the distances which are then stored in an Excel file. ...

Checking row and column for a word in python

I am trying to create a checking program to see if the word is in a matrix horizontally or vertically. I have the code for checking the row, but would checking the column be similar to the row code? ...

热门标签