matlab - Producto triple de Matrix/Tensor?

numpy matrix-multiplication

Question

Un algoritmo en el que estoy trabajando requiere computar, en un par de lugares, un tipo de producto triple matriz.

La operación toma tres matrices cuadradas con dimensiones idénticas, y produce un tensor de 3 índices. Etiquetando los operandos A , B y C , el elemento (i,j,k) -th del resultado es

X[i,j,k] = /sum_a A[i,a] B[a,j] C[k,a]

En numpy, puedes calcular esto con einsum(''ia,aj,ka->ijk'', A, B, C) .

Preguntas:

• ¿Esta operación tiene un nombre estándar?
• ¿Puedo calcular esto con una sola llamada BLAS?
• ¿Hay alguna otra biblioteca C / Fortran numérica fuertemente optimizada que pueda calcular expresiones de este tipo?

Answer 1

Introducción y código de solución

np.einsum , es realmente difícil de superar, pero en casos raros, todavía puede vencerlo, si puede traer matrix-multiplication en los cálculos. Después de algunos ensayos, parece que puede traer matrix-multiplication with np.dot para superar el rendimiento con np.einsum(''ia,aj,ka->ijk'', A, B, C) .

La idea básica es dividir la operación "all einsum" en una combinación de np.einsum y np.dot como se np.dot continuación:

• Las sumas para A:[i,a] y B:[a,j] se hacen con np.einsum para obtener una 3D array:[i,j,a] .
• Esta matriz 3D se vuelve a formar en una 2D array:[i*j,a] y la tercera matriz, C[k,a] se transpone a [a,k] , con la intención de realizar matrix-multiplication entre estos dos, dándonos [i*j,k] como el producto de la matriz, ya que perdemos el índice [a] allí.
• El producto se reconfigura en una 3D array:[i,j,k] para la salida final.

Aquí está la implementación de la primera versión discutida hasta ahora:

import numpy as np

def tensor_prod_v1(A,B,C):   # First version of proposed method
    # Shape parameters
    m,d = A.shape
    n = B.shape[1]
    p = C.shape[0]

    # Calculate /sum_a A[i,a] B[a,j] to get a 3D array with indices as (i,j,a)
    AB = np.einsum(''ia,aj->ija'', A, B)

    # Calculate entire summation losing a-ith index & reshaping to desired shape
    return np.dot(AB.reshape(m*n,d),C.T).reshape(m,n,p)

Como estamos sumando el índice a a-th en las tres matrices de entrada, uno puede tener tres métodos diferentes para sumar a lo largo del índice a-ésimo. El código enumerado anteriormente fue para (A,B) . Por lo tanto, también podemos tener (A,C) y (B,C) dándonos dos variaciones más, como se detalla a continuación:

def tensor_prod_v2(A,B,C):
    # Shape parameters
    m,d = A.shape
    n = B.shape[1]
    p = C.shape[0]

    # Calculate /sum_a A[i,a] C[k,a] to get a 3D array with indices as (i,k,a)
    AC = np.einsum(''ia,ja->ija'', A, C)

    # Calculate entire summation losing a-ith index & reshaping to desired shape
    return np.dot(AC.reshape(m*p,d),B).reshape(m,p,n).transpose(0,2,1)

def tensor_prod_v3(A,B,C):
    # Shape parameters
    m,d = A.shape
    n = B.shape[1]
    p = C.shape[0]

    # Calculate /sum_a B[a,j] C[k,a] to get a 3D array with indices as (a,j,k)
    BC = np.einsum(''ai,ja->aij'', B, C)

    # Calculate entire summation losing a-ith index & reshaping to desired shape
    return np.dot(A,BC.reshape(d,n*p)).reshape(m,n,p)

Dependiendo de las formas de las matrices de entrada, los diferentes enfoques producirían diferentes aceleraciones entre sí, pero tenemos la esperanza de que todo sea mejor que el enfoque de all-einsum . Los números de rendimiento se enumeran en la siguiente sección.

Pruebas de tiempo de ejecución

Esta es probablemente la sección más importante, ya que tratamos de ver las cifras de aceleración con las tres variaciones del enfoque propuesto sobre el enfoque de all-einsum como se propuso originalmente en la pregunta.

Conjunto de datos n. ° 1 (matrices de forma igual):

In [494]: L1 = 200
     ...: L2 = 200
     ...: L3 = 200
     ...: al = 200
     ...: 
     ...: A = np.random.rand(L1,al)
     ...: B = np.random.rand(al,L2)
     ...: C = np.random.rand(L3,al)
     ...: 

In [495]: %timeit tensor_prod_v1(A,B,C)
     ...: %timeit tensor_prod_v2(A,B,C)
     ...: %timeit tensor_prod_v3(A,B,C)
     ...: %timeit np.einsum(''ia,aj,ka->ijk'', A, B, C)
     ...: 
1 loops, best of 3: 470 ms per loop
1 loops, best of 3: 391 ms per loop
1 loops, best of 3: 446 ms per loop
1 loops, best of 3: 3.59 s per loop

Dataset # 2 (Bigger A):

In [497]: L1 = 1000
     ...: L2 = 100
     ...: L3 = 100
     ...: al = 100
     ...: 
     ...: A = np.random.rand(L1,al)
     ...: B = np.random.rand(al,L2)
     ...: C = np.random.rand(L3,al)
     ...: 

In [498]: %timeit tensor_prod_v1(A,B,C)
     ...: %timeit tensor_prod_v2(A,B,C)
     ...: %timeit tensor_prod_v3(A,B,C)
     ...: %timeit np.einsum(''ia,aj,ka->ijk'', A, B, C)
     ...: 
1 loops, best of 3: 442 ms per loop
1 loops, best of 3: 355 ms per loop
1 loops, best of 3: 303 ms per loop
1 loops, best of 3: 2.42 s per loop

Dataset # 3 (Bigger B):

In [500]: L1 = 100
     ...: L2 = 1000
     ...: L3 = 100
     ...: al = 100
     ...: 
     ...: A = np.random.rand(L1,al)
     ...: B = np.random.rand(al,L2)
     ...: C = np.random.rand(L3,al)
     ...: 

In [501]: %timeit tensor_prod_v1(A,B,C)
     ...: %timeit tensor_prod_v2(A,B,C)
     ...: %timeit tensor_prod_v3(A,B,C)
     ...: %timeit np.einsum(''ia,aj,ka->ijk'', A, B, C)
     ...: 
1 loops, best of 3: 474 ms per loop
1 loops, best of 3: 247 ms per loop
1 loops, best of 3: 439 ms per loop
1 loops, best of 3: 2.26 s per loop

Dataset # 4 (Bigger C):

In [503]: L1 = 100
     ...: L2 = 100
     ...: L3 = 1000
     ...: al = 100
     ...: 
     ...: A = np.random.rand(L1,al)
     ...: B = np.random.rand(al,L2)
     ...: C = np.random.rand(L3,al)

In [504]: %timeit tensor_prod_v1(A,B,C)
     ...: %timeit tensor_prod_v2(A,B,C)
     ...: %timeit tensor_prod_v3(A,B,C)
     ...: %timeit np.einsum(''ia,aj,ka->ijk'', A, B, C)
     ...: 
1 loops, best of 3: 250 ms per loop
1 loops, best of 3: 358 ms per loop
1 loops, best of 3: 362 ms per loop
1 loops, best of 3: 2.46 s per loop

Conjunto de datos n. ° 5 (longitud de la dimensión a-ésima más grande):

In [506]: L1 = 100
     ...: L2 = 100
     ...: L3 = 100
     ...: al = 1000
     ...: 
     ...: A = np.random.rand(L1,al)
     ...: B = np.random.rand(al,L2)
     ...: C = np.random.rand(L3,al)
     ...: 

In [507]: %timeit tensor_prod_v1(A,B,C)
     ...: %timeit tensor_prod_v2(A,B,C)
     ...: %timeit tensor_prod_v3(A,B,C)
     ...: %timeit np.einsum(''ia,aj,ka->ijk'', A, B, C)
     ...: 
1 loops, best of 3: 373 ms per loop
1 loops, best of 3: 269 ms per loop
1 loops, best of 3: 299 ms per loop
1 loops, best of 3: 2.38 s per loop

Conclusiones: estamos viendo una aceleración de 8x-10x con las variaciones del enfoque propuesto sobre el enfoque de all-einsum se menciona en la pregunta.

Answer 2

Deje n x n ser los tamaños de la matriz. En Matlab, puedes

1. Grupo A y C en una matriz n^2 x n AC , de manera que las filas de AC correspondan a todas las combinaciones de filas de A y C
2. Publicar-multiplicar AC por B Eso da el resultado deseado, solo que en una forma diferente.
3. Reforme y permute las dimensiones para obtener el resultado en la forma deseada.

Código:

AC = reshape(bsxfun(@times, permute(A, [1 3 2]), permute(C, [3 1 2])), n^2, n); % // 1
X = permute(reshape((AC*B).'', n, n, n), [2 1 3]);                               %''// 2, 3

Verifique con un enfoque literal basado en bucle:

%// Example data:
n = 3;
A = rand(n,n);
B = rand(n,n);
C = rand(n,n);

%// Proposed approach:
AC = reshape(bsxfun(@times, permute(A, [1 3 2]), permute(C, [3 1 2])), n^2, n);
X = permute(reshape((AC*B).'', n, n, n), [2 1 3]); %''

%// Loop-based approach:
Xloop = NaN(n,n,n); %// initiallize
for ii = 1:n
    for jj = 1:n
        for kk = 1:n
            Xloop(ii,jj,kk) = sum(A(ii,:).*B(:,jj).''.*C(kk,:)); %''
        end
    end
end

%// Compute maximum relative difference:
max(max(max(abs(X./Xloop-1))))

ans =
    2.2204e-16

La diferencia relativa máxima es del orden de eps , por lo que el resultado es correcto dentro de la precisión numérica.

Answer 3

Sé que esto es un poco viejo ahora, pero este tema surge mucho. En Matlab es difícil vencer a tprod, un archivo MEX escrito por Jason Farquhar disponible aquí

https://www.mathworks.com/matlabcentral/fileexchange/16275-tprod-arbitary-tensor-products-between-and-arrays

tprod funciona muy parecido al einsum aunque está limitado a una operación binaria (2 tensores). Probablemente esto no sea realmente una limitación porque sospecho que einsum simplemente realiza una serie de operaciones binarias. El orden de estas operaciones hace una gran diferencia y tengo entendido que einsum simplemente las ejecuta en el orden en que se transfieren las matrices y no permite productos intermedios múltiples.

tprod también está limitado a matrices densas (completas). El Tensor Toolbox de Kolda (mencionado en una publicación anterior) admite tensores dispersos, pero tiene una funcionalidad más limitada que tprod (no permite índices repetidos en la salida). Estoy trabajando en llenar estos vacíos, pero ¿no sería bueno si Mathworks lo hiciera?