Geometric Matrix Multiplication

Question

I'm working on a project building a geometric (as opposed to arithmetic) neural net. To construct the transfer function, I would like to use geometric summation instead of arithmetic summation.

To make things clearer, I'm just going to 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, it's a pretty minimal change. The only problem is that, in the same way I would never actually write and use the above arithmetic code (I would use numpy.dot), I would really not like to actually use the geometric code above. Is there any way to leverage numpy's matrix multiplication to achieve the geometric result? I haven't been able to think of one, and I haven't found anything obvious past the solution above, which is far from optimal.

Answer 1

You are all complicating things unnecessarily... Since you have a single operation, multiplication, which is commutative, you can swap the order in which you perform them at will, i.e. you don't need to multiply items of matrix1 with items of matrix2, and once you have computed all of them, multiply them together. Instead, you can first multiply all the relevant items of matrix1 together, then all the relevant items of matrix2 together, and then multiply the two resulting values. So you can write your function as the very simple:

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

It has the additional advantage that, if you are multiplying matrices of shapes (m, k) and (k, n), you would expect to be having to do m*n*2*k multiplications, while this method only requires m*k + n*k + m*n, which is bound to be much smaller than what you presently are doing for almost any array shapes.

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, :]*y.T[..., None], axis=2)
np.allclose(xy, xy2)

## -- End pasted text --
Out[373]: True

This solution is not very stable with regards of what shapes it can handle. So great if it works, but not something you should use for anything long term, if your data size is variable.

Answer 3

Your problem is a quite perfect candidate for numba. The only thing you should have to add is @autojit. Of course you would also optimize the nested calls for the iteration length. range and xrange are treated both as simple for-loops by numba (without the overhead to create large arrays). In the end it could 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

By making use of jit compiling for this function, you should get C-like speed afterwards. 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