Multiplique uma matriz 3D com uma matriz 2D

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

Pergunta

Suponha que eu tenho um AXBXC matriz X e a Bxd matriz Y.

Existe um método não loop pelo qual eu posso multiplicar cada um dos C AXB matrizes com Y?

Foi útil?

Solução

Você pode fazer isso em uma linha usando as funções Num2Cell Para quebrar a matriz X em uma matriz celular e Cellfun Para operar através das células:

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

O resultado Z é um 1 por C. matriz celular onde cada célula contém um A-BY-D matriz. Se você quiser Z ser um A-BY-D-BY-C matriz, você pode usar o GATO função:

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


NOTA: Minha solução antiga usada Mat2Cell ao invés de Num2Cell, que não foi tão sucinto:

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

Outras dicas

Como preferência pessoal, gosto que meu código seja o mais sucinto e legível possível.

Aqui está o que eu teria feito, embora não atenda ao seu requisito de 'loops':

for m = 1:C

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

end

Isso resulta em um A X D X C matriz Z.

E, claro, você sempre pode pré-alocar z para acelerar as coisas usando Z = zeros(A,D,C);.

Aqui está uma solução de uma linha (dois se você quiser se dividir na 3ª dimensão):

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

Daí agora: Z(:,:,i) contém o resultado de X(:,:,i) * Y

Explicação:

O exposto acima pode parecer confuso, mas a ideia é simples. Primeiro eu começo por tomar a terceira dimensão de X e faça uma concatenação vertical ao longo da primeira dim:

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

... a dificuldade era que C é uma variável, portanto, você não pode generalizar essa expressão usando gato ou Vertcat. Em seguida, multiplicamos isso por Y:

ZZ = XX * Y;

Finalmente, divida -o de volta à terceira dimensão:

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

Então você pode ver isso apenas requer uma multiplicação da matriz, mas você precisa remodelar a matriz antes e depois.

Estou abordando exatamente o mesmo problema, com um olho no método mais eficiente. Existem aproximadamente três abordagens que eu vejo, sem usar bibliotecas externas (ou seja, mtimesx):

  1. Loop através de fatias da matriz 3D
  2. Magia Repmat-and-Permute
  3. Multiplicação de Cellfun

Recentemente, comparei todos os três métodos para ver qual foi o mais rápido. Minha intuição era que (2) seria o vencedor. Aqui está o 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

Todas as três abordagens produziram a mesma saída (ufa!), Mas, surpreendentemente, o loop foi o mais rápido:

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

Observe que os tempos podem variar bastante de um teste para outro e, às vezes (2), sai mais lento. Essas diferenças se tornam mais dramáticas com dados maiores. Mas com Muito de dados maiores, (3) batidas (2). O método de loop ainda é o melhor.

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

Mas o método de loop posso Seja mais lento que (2), se a dimensão em loop for muito maior que os outros.

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.

Então (2) vence por um grande fator, neste caso (talvez extremo). Pode não haver uma abordagem ideal em todos os casos, mas o loop ainda é muito bom e melhor em muitos casos. Também é melhor em termos de legibilidade. Faça uma garra!

Não. Existem várias maneiras, mas sempre sai em um loop, direto ou indireto.

Só para agradar minha curiosidade, por que você quer isso?

Para responder à pergunta, e Para legibilidade, consulte:

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

Entrada

  • 2 matrizes
  • escurecido

Exemplo

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

Fonte

Fonte original. Eu adicionei comentários embutidos.

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

Eu recomendo que você use o Caixa de ferramentas MMX de Matlab. Pode multiplicar as matrizes N-dimensionais o mais rápido possível.

As vantagens de Mmx são:

  1. Isso é fácil usar.
  2. Multiplicar matrizes n-dimensionais (Na verdade, ele pode multiplicar matrizes de matrizes 2-D)
  3. Ele executa outros operações da matriz (Transpor, multiplique quadrático, decomposição de CHOL e muito mais)
  4. Ele usa C compilador e Multi fio Computação para acelerar.

Para esse problema, você só precisa escrever este comando:

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

Aqui está uma referência para todos os métodos possíveis. Para mais detalhes, consulte isso pergunta.

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

Eu pensaria recursão, mas esse é o único outro método que não é loop que você pode fazer

Você poderia "desmoltar" o loop, ou seja, escrever todas as multiplicações sequencialmente que ocorreriam no loop

Eu gostaria de compartilhar minha resposta para os problemas de:

1) fazer o produto tensor de dois tensores (de qualquer valência);

2) Fazendo a contração de dois tensores ao longo de qualquer dimensão.

Aqui estão minhas sub -rotinas para a primeira e a segunda tarefas:

1) Produto Tensor:

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

2) Contração: Aqui A e B são os tensores a serem contratados ao longo das dializações I e J, respectivamente. Os comprimentos dessas dimensões devem ser iguais, é claro. Não há verificação para isso (isso obscureceria o código), mas, além disso, funciona bem.

   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

Algum crítico?

Licenciado em: CC-BY-SA com atribuição
Não afiliado a StackOverflow
scroll top