Geometric Matrix Multiplication

Question

I am working on a project creating a geometric (as opposed to arithmetic) neural network. To construct the transfer function, I would like to use geometric summation instead of arithmetic summation.

To make things clearer, I will simply describe in code:

def arithmetic_matrix_multiplication(matrix1, matrix2): new_matrix = np.zeros(len(matrix1),len(matrix2[0])) for i in range(len(matrix1)): for j in range(len(matrix2[0])): for k in range(len(matrix2)): new_matrix[i][j] += matrix1[i][k]*matrix2[k][j] return new_matrix def geometric_matrix_multiplication(matrix1, matrix2): new_matrix = np.ones(len(matrix1),len(matrix2[0])) for i in range(len(matrix1)): for j in range(len(matrix2[0])): for k in range(len(matrix2)): new_matrix[i][j] *= matrix1[i][k]*matrix2[k][j] return new_matrix 

As you can see, this is a pretty minimal change. The only problem is that in the same way I would never write and use the above arithmetic code (I would use numpy.dot ), I really would not want to use the geometric code above. Is there a way to use matrix matrix multiplication to achieve a geometric result? I could not think of one thing, and I did not find anything obvious in the past above, which is far from optimal.

Answer 1

You all complicate things unnecessarily ... Since you have one operation, multiplication, which is commutative, you can change the order in which you perform them at will, i.e. you do not need to multiply matrix1 elements with matrix1 elements, and once you have calculated all of them, multiply them together. Instead, you can first multiply all the corresponding elements of matrix1 together, then all the corresponding elements of matrix2 together, and then multiply the two resulting values. This way you can write your function as very simple:

 def fast_geometric_matrix_multiplication(matrix1, matrix2): return np.prod(matrix1, axis=1)[:, None] * np.prod(matrix2, axis=0) 

It has an additional advantage: if you multiply the matrices of the forms (m, k) and (k, n) , you expect that you will have to do the multiplications m*n*2*k , whereas this method requires only m*k + n*k + m*n , which is bound to be much smaller than what you are currently doing for almost any array form.

And of course:

 In [24]: a = np.random.rand(100, 200) ...: b = np.random.rand(200, 50) ...: In [25]: np.allclose(geometric_matrix_multiplication(a, b), ...: fast_geometric_matrix_multiplication(a, b)) Out[25]: True In [26]: %timeit geometric_matrix_multiplication(a, b) 1 loops, best of 3: 1.39 s per loop In [27]: %timeit fast_geometric_matrix_multiplication(a, b) 10000 loops, best of 3: 74.2 us per loop In [28]: %timeit np.prod(a[:, None, :]*b[..., None].T, axis=2) 100 loops, best of 3: 5.97 ms per loop 

Answer 2

 In [373]: %paste x = np.arange(5*5, dtype=np.long).reshape(5,5) y = np.arange(5*5, 5*5*2, dtype=np.long).reshape(5,5) xy = geometric_matrix_multiplication(x,y) xy2 = np.prod(x[:, None, :]*yT[..., None], axis=2) np.allclose(xy, xy2) ## -- End pasted text -- Out[373]: True 

This solution is not very stable as to what forms it can handle. So great if it works, but not something you should use for anything in the long run if your data size is variable.

Answer 3

Your problem is the perfect candidate for numba . The only thing you need to add: @autojit . Of course, you would also optimize nested calls for the iteration length. range and xrange treated as simple for-loops numba (no overhead for creating large arrays). In the end, it might look like this:

 from numba import autojit @autojit def geometric_matrix_multiplication(matrix1, matrix2): new_matrix = np.ones((len(matrix1),len(matrix2[0]))) jlen = len(matrix2[0]) klen = len(matrix2) for i in range(len(matrix1)): for j in range(jlen): for k in range(klen): new_matrix[i,j] *= matrix1[i,k]*matrix2[k,j] return new_matrix 

Using jit compilation for this function, you should get a C-like speed later. To give some idea of ​​the speed gain:

 a = np.random.rand(100*100).reshape((100,100)) b = np.random.rand(100*100).reshape((100,100)) %timeit geometric_matrix_multiplication_org(a,b) 1 loops, best of 3: 4.1 s per loop %timeit geometric_matrix_multiplication(a,b) 100 loops, best of 3: 5 ms per loop