Ponga las diferencias por pares de las filas de la matriz en una matriz tridimensional

Tengo una matriz Y de forma (n, d). Ya calculé las diferencias entre filas por pares de la siguiente manera:

I, J = np.triu_indices(Y.shape[0], 0) rowDiffs = (Y[I, :] - Y[J, :]) 

No, quiero crear una matriz 3D que contenga las diferencias de las filas i y j de Y en la posición (i, j, :). ¿Como lo harias?

El objective de esto es reemplazar este bucle ineficiente:

  for i in range(Y.shape[0]): for j in range(Y.shape[0]): C[i,:] = C[i,:] + W[i, j] * (Y[i, :]-Y[j, :]) 

He encontrado cierto éxito con esto:

 row_diffs = Y[:, np.newaxis] - Y 

Y[:, np.newaxis] crea una versión de Y con dimensiones (n, 1, 3). Luego, la resta usa la transmisión para hacer lo que quieras.

Desafortunadamente, he encontrado que este enfoque es relativamente lento y todavía no he encontrado una mejor manera.

Ejemplo completo:

 >>> x = np.random.randint(10, size=(4, 3)) >>> x array([[4, 0, 8], [8, 5, 3], [4, 1, 6], [2, 2, 4]]) >>> x[:, np.newaxis] - x array([[[ 0, 0, 0], [-4, -5, 5], [ 0, -1, 2], [ 2, -2, 4]], [[ 4, 5, -5], [ 0, 0, 0], [ 4, 4, -3], [ 6, 3, -1]], [[ 0, 1, -2], [-4, -4, 3], [ 0, 0, 0], [ 2, -1, 2]], [[-2, 2, -4], [-6, -3, 1], [-2, 1, -2], [ 0, 0, 0]]]) 

Aquí hay un enfoque vectorizado que usa broadcasting y np.einsum

 np.einsum('ij,ijk->ik',W,Y[:,None] - Y) 

Prueba de tiempo de ejecución –

 In [29]: def original_app(Y,W): ...: m = Y.shape[0] ...: C = np.zeros((m,m)) ...: for i in range(Y.shape[0]): ...: for j in range(Y.shape[0]): ...: C[i,:] = C[i,:] + W[i, j] * (Y[i, :]-Y[j, :]) ...: return C ...: In [30]: # Inputs ...: Y = np.random.rand(100,100) ...: W = np.random.rand(100,100) ...: In [31]: out = original_app(Y,W) In [32]: np.allclose(out, np.einsum('ij,ijk->ik',W,Y[:,None] - Y)) Out[32]: True In [33]: %timeit original_app(Y,W) 10 loops, best of 3: 70.8 ms per loop In [34]: %timeit np.einsum('ij,ijk->ik',W,Y[:,None] - Y) 100 loops, best of 3: 4.01 ms per loop