Vectorize calculating all unit vectors for a set of points in numpy

Question

I need to calculate all unit vectors between two sets of points.

Currently I have this:

def all_unit_vectors(points_a, points_b):
    results = np.zeros((len(points_a) * len(points_b), 3), dtype=np.float32)
    
    count = 0
    for pt_a in points_a:
        for pt_b in points_b:
            results[count] = (pt_a - pt_b)/np.linalg.norm([pt_a - pt_b])
            count += 1
            
    return results
            
       
in_a = np.array([[51.34, 63.68, 7.98], 
                 [53.16, 63.23, 7.19],
                 [77.50, 62.55, 4.23],
                 [79.54, 62.73, 3.61]])

in_b = np.array([[105.58, 61.09, 5.50],
                 [107.37, 60.66, 6.50],
                 [130.73, 58.30, 12.33],
                 [132.32, 58.48, 13.38]])


results = all_unit_vectors(in_a, in_b)

print(results)

which (correctly) outputs:

[[-0.368511    0.01759667  0.01684932]
 [-0.3777128   0.02035861  0.00997707]
 [-0.47964868  0.03250422 -0.02628129]
 [-0.4851439   0.03115273 -0.03235091]
 [-0.3551545   0.01449887  0.01145004]
 [-0.3644423   0.01727756  0.00463872]
 [-0.46762046  0.02971985 -0.03098581]
 [-0.4732132   0.02839518 -0.03700341]
 [-0.17814296  0.00926242 -0.00805704]
 [-0.18821244  0.01190899 -0.01430339]
 [-0.3044056   0.02430441 -0.04632135]
 [-0.31113514  0.0230996  -0.05193153]
 [-0.16408844  0.0103343  -0.01190965]
 [-0.1741932   0.01295652 -0.01808905]
 [-0.29113463  0.02519489 -0.04959355]
 [-0.29793915  0.02399093 -0.05515092]]

Can the loops in all_unit_vectors() be vectorized?

*"I need to calculate all unit vectors between two sets of points."* The rows in your sample output are not unit vectors. Why do you divide `pt_a - pt_b` by `np.linalg.norm([pt_a, pt_b])`? — Warren Weckesser, Jul 06 '20 at 04:44
@WarrenWeckesser isn't it? I got that from https://stackoverflow.com/a/52781450/280795 (although I was unclear what `v` looks like in that answer) — Nick Lothian, Jul 06 '20 at 04:49
`np.linalg.norm([-0.368511, 0.01759667, 0.01684932]) ` is 0.369..., not 1. To get a unit vector, you would have to divide by `np.linalg.norm(pt_a - pt_b)`. The difference `pt_a - pt_b` corresponds to `v` in the answer that you linked to, so you have to replace both occurrences of `v` with `pt_a - pt_b`. — Warren Weckesser, Jul 06 '20 at 05:02
@WarrenWeckesser oh of course. Thanks for picking that up. I've edited the question so no one else makes the same mistake. — Nick Lothian, Jul 06 '20 at 08:22

score 2 · Accepted Answer · answered Jul 06 '20 at 03:54

norm is calculated as root of sum squared, you can implement your own norm calculation as follows, and then vectorize your solution with broadcasting:

diff = (in_a[:, None] - in_b).reshape(-1, 3)
norm = ((in_a[:, None] ** 2 + in_b ** 2).sum(2) ** 0.5).reshape(-1, 1)

diff / norm

gives:

[[-0.36851098  0.01759667  0.01684932]
 [-0.3777128   0.02035861  0.00997706]
 [-0.47964868  0.03250422 -0.02628129]
 [-0.4851439   0.03115273 -0.03235091]
 [-0.35515452  0.01449887  0.01145004]
 [-0.36444229  0.01727756  0.00463872]
 [-0.46762047  0.02971985 -0.03098581]
 [-0.4732132   0.02839518 -0.03700341]
 [-0.17814297  0.00926242 -0.00805704]
 [-0.18821243  0.01190899 -0.01430339]
 [-0.30440561  0.02430441 -0.04632135]
 [-0.31113513  0.0230996  -0.05193153]
 [-0.16408845  0.0103343  -0.01190965]
 [-0.1741932   0.01295652 -0.01808905]
 [-0.29113461  0.02519489 -0.04959355]
 [-0.29793917  0.02399093 -0.05515092]]

Play.

Vectorize calculating all unit vectors for a set of points in numpy

1 Answers1