1

Basically, I have two matrices A and B, and I want C (dimensions marked by the side of the matrices), with computation like this:

The formula below is what I do now. I take advantage of some broadcasting, but I am still left with a loop. I am novel to Python so maybe I am wrong, but I just have a hunch that this loop can be eliminated. Can anyone share some ideas?

EDIT: 2018-04-27 09:48:28 as requested, an example:

In [5]: A
Out[5]: 
array([[0, 1],
       [2, 3],
       [4, 5],
       [6, 7],
       [8, 9]])

In [6]: B
Out[6]: 
array([[0, 1],
       [2, 3],
       [4, 5],
       [6, 7]])

In [7]: C = np.zeros ((B.shape[0], A.shape[0]))

In [8]: for m in range (B.shape[0]):
   ...:     C[m] = np.sum (np.square (B[m] - A), axis=1).flatten ()
   ...:     

In [9]: C
Out[9]: 
array([[  0.,   8.,  32.,  72., 128.],
       [  8.,   0.,   8.,  32.,  72.],
       [ 32.,   8.,   0.,   8.,  32.],
       [ 72.,  32.,   8.,   0.,   8.]])
qzhangjhu
  • 131
  • 1
  • 6
  • 1
    It's traditional to give an example of the slower way, e.g. something like `A = np.random.random((D, N)); B = np.random.random((D, M)); C = np.zeros((N, M)); for m in range(M): C[m] = np.sum((B[m]-A)**2, axis=1).flatten() ` (or whatever the actual formula is.) – DSM Apr 27 '18 at 13:33

1 Answers1

1

This appears to work at the cost of some extra memory:

C = ((B[:, :, None] - A.T)**2).sum(axis=1)

Testing:

import numpy
D = 10
N = 20
M = 30

A = numpy.random.rand(N, D)
B = numpy.random.rand(M, D)
C = numpy.empty((M, N))

Timing:

for m in range(M):
    C[m] = numpy.sum((B[m, :] - A)**2, axis=1)

514 µs ± 13.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

C2 = ((B[:, :, None] - A.T)**2).sum(axis=1)

53.6 µs ± 529 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

chthonicdaemon
  • 19,180
  • 2
  • 52
  • 66