Multiplicar una matriz 3D con una matriz 2D

StackOverflow https://stackoverflow.com/questions/1745299

Pregunta

Supongamos que tengo un AxBxC X matriz y un BxD Y matriz.

¿Hay un método no bucle por la que puedo multiplicar cada uno de los C AxB matrices con Y?

¿Fue útil?

Solución

Usted puede hacer esto en una línea usando las funciones NUM2CELL para romper el X matriz en una matriz celular y CELLFUN para operar a través de las células: 

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

El Z resultado es un 1-por-C serie de células, donde cada célula contiene un matriz A-por-D. Si quieres Z a ser un A-por-D-by-C matriz, se puede utilizar el función CAT : 

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


Nota: Mi vieja solución utiliza MAT2CELL en lugar de NUM2CELL , que no era lo más sucinta: 

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

Otros consejos

Como una preferencia personal, me gusta mi código sea lo más breve y fácil de leer como sea posible.

Esto es lo que habría hecho, a pesar de que no cumple con sus necesidades 'no-loops': 

for m = 1:C

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

end

Esto resulta en un A x D x C matriz Z .

Y, por supuesto, siempre se puede pre-asignar Z para acelerar las cosas mediante el uso de Z = zeros(A,D,C);.

Aquí es una solución de una línea (dos si desea dividir en tercera dimensión): 

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

Por lo tanto ahora: Z(:,:,i) contiene el resultado de X(:,:,i) * Y

Explicación:

Lo anterior puede parecer confuso, pero la idea es simple. Primero comienzo de tomar la tercera dimensión de la X y hacer una concatenación vertical a lo largo del primer tenue: 

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

... la dificultad era que C es una variable, por lo tanto, no se puede generalizar que el uso de la expresión cat o vertcat . A continuación multiplicamos esto por Y: 

ZZ = XX * Y;

Finalmente he dividido de nuevo en la tercera dimensión: 

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

Así se puede ver que sólo requiere una multiplicación de matrices, pero hay que reshape la matriz antes y después.

Me estoy acercando el mismo problema exacto, con un ojo para el método más eficiente. Hay aproximadamente tres enfoques que ven a su alrededor, a falta de uso de bibliotecas externas (es decir, mtimesx ):

  1. Loop a través de las rebanadas de la matriz 3D
  2. repmat-y-permute magia
  3. cellfun multiplicación

Recientemente he comparado los tres métodos para ver cuál fue el más rápido. Mi intuición era que (2) sería el ganador. Aquí está el código: 

% 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

Los tres enfoques producen la misma salida (¡uf!), Pero, sorprendentemente, el bucle fue el más rápido: 

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

Tenga en cuenta que los tiempos pueden variar mucho de un ensayo a otro, ya veces (2) sale el más lento. Estas diferencias se hacen más dramático con los datos más grandes. Pero con mucho de datos más grandes, (3) tiempos (2). El método de bucle es todavía mejor. 

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

Pero el método de bucle puede ser más lenta que (2), si la dimensión de bucle es mucho más grande que los otros. 

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.

Así que (2) gana por un factor importante, en este caso (tal extremo). Puede que no haya un enfoque que sea óptima en todos los casos, pero el bucle sigue siendo bastante bueno, y lo mejor en muchos casos. También es mejor en términos de facilidad de lectura. Bucle de distancia!

Nop. Hay varias maneras, pero siempre sale en un bucle, directo o indirecto.

Sólo para complacer mi curiosidad, ¿por qué quieres que de todos modos?

Para responder a la pregunta, y para facilitar la lectura, por favor ver:

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

entrada

  • 2 arrays
  • dim

Ejemplo 

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

Fuente

Fuente original. Añadí comentarios en línea. 

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

Le recomiendo que utilice el MMX caja de herramientas de MATLAB. Se puede multiplicar matrices n-dimensionales tan rápido como sea posible.

Las ventajas de MMX son:

  1. fácil a utilizar.
  2. Multiplicar n-dimensional matrices (en realidad puede multiplicar matrices de matrices 2-D)
  3. Se lleva a cabo otras operaciones matriz (transponer, Quadratic Multiply, la descomposición Chol y más)
  4. Utiliza C compilador y multi-hilo cálculo para la velocidad hacia arriba.

Para este problema, sólo tiene que escribir este comando: 

C=mmx('mul',X,Y);

aquí es un punto de referencia para todos los métodos posibles. Para más detalles consulte esta pregunta . 

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

Me gustaría pensar recursividad, pero ese es el único otro método no se puede hacer bucle

Se puede "desenrollar" del bucle, es decir, escribir toda la secuencialmente multiplicaciones que se producirían en el bucle

Me gustaría compartir mi respuesta a los problemas de:

1) haciendo que el producto tensorial de dos tensores (de cualquier valencia);

2) hacer que la contracción de dos tensores a lo largo de cualquier dimensión.

Aquí están mis subrutinas para el primer y segundo tareas:

1) producto tensorial: 

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

2) contracción: Aquí A y B son los tensores para ser contratados a lo largo de los dimesions I y J respectivamente. Las longitudes de estas dimensiones deben ser iguales, por supuesto. No hay verificación para esta (esto oscurecería el código), pero aparte de esto funciona bien. 

   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

Cualquier críticos?

Licenciado bajo: CC-BY-SA con atribución
No afiliado a StackOverflow
scroll top